Classification of magnetized star–planet interactions:
bow shocks, tails, and inspiraling flows
Key Words.:
Planetstar interactions  Stars: winds, outflows  Magnetohydrodynamics (MHD)  Methods: numericalAbstract
Context: Closein exoplanets interact with their host stars gravitationally as well as via their magnetized plasma outflows. The rich dynamics that arises may result in distinct observable features.
Aims: Our objective is to study and classify the morphology of the different types of interaction that can take place between a giant closein planet (a Hot Jupiter) and its host star, based on the physical parameters that characterize the system.
Methods: We perform 3D magnetohydrodynamic numerical simulations to model the star–planet interaction, incorporating a star, a Hot Jupiter, and realistic stellar and planetary outflows. We explore a wide range of parameters and analyze the flow structures and magnetic topologies that develop.
Results: Our study suggests the classification of star–planet interactions into four general types, based on the relative magnitudes of three characteristic length scales that quantify the effects of the planetary magnetic field, the planetary outflow, and the stellar gravitational field in the interaction region. We describe the dynamics of these interactions and the flow structures that they give rise to, which include bow shocks, cometarytype tails, and inspiraling accretion streams. We point out the distinguishing features of each of the classified cases and discuss some of their observationally relevant properties.
Conclusions: The magnetized interactions of star–planet systems can be categorized, and their general morphologies predicted, based on a set of basic stellar, planetary, and orbital parameters.
1 Introduction
The exoplanet class “Hot Jupiters” consists of massive gaseous planets (of mass –, where the subscript “J” denotes the planet Jupiter) that orbit their host stars on closein trajectories (semimajor axis of , corresponding to an orbital period of a few days). Among other exoplanets, the size and location of Hot Jupiters favors both their discovery and the collection of detailed data. The two leading and complementary detection techniques, radial velocity and transit observations, have resulted in a wealth of information regarding their orbital characteristics, such as obliquity and eccentricity (e.g., Fabrycky & Winn 2009), as well as their physical properties, such as surface temperature and atmospheric composition (e.g., Bean et al. 2013). The study of the phenomenology of Hot Jupiters also provides valuable input to research on the physical conditions and processes in smaller, harder to probe, Neptune and Earthlike bodies (see, e.g., Seager 2011 for a recent reference on exoplanets).
Hot Jupiter atmospheres are believed to be heated by photoionizing extremeultraviolet (EUV) stellar irradiation (e.g., Burrows et al. 2000, but see also Owen & Jackson 2012 as well as Buzasi 2013 and Lanza 2013, who proposed, respectively, Xrays for very young stars and magnetic reconnection for very closein planets). The energy input could be large enough to induce an evaporative mass loss at a rate of –. Such outflows have been inferred in a couple of systems (HD 189733b and HD 209458b) in which detailed observations have been carried out (e.g., VidalMadjar et al. 2004; Ehrenreich et al. 2008; Lecavelier des Etangs et al. 2010, 2012; Bourrier et al. 2013, although see BenJaffel 2007, 2008 for a different interpretation of the data). The theory of evaporative planetary outflows has been developed in a number of papers, taking into account the thermal and ionization structure of the irradiated planetary atmosphere as well as tidal and magnetic effects (e.g., Lecavelier des Etangs et al. 2004; Yelle 2004, 2006; Tian et al. 2005; García Muñoz 2007; MurrayClay et al. 2009; Adams 2011; Trammell et al. 2011, 2014; Koskinen et al. 2013; Owen & Adams 2014, and references therein).
The host stars of closein planets typically drive magnetized stellar winds that are accelerated to speeds of a few hundred . The properties of these winds change as the star evolves, from strong and often collimated outflows in young stars (e.g., Matt & Pudritz 2008; Matt et al. 2012b) to weak and more isotropic winds during the mainsequence phase (e.g., Matt et al. 2012a; Cohen & Drake 2014, and references therein). The magnetized plasma that is driven out in these winds constitutes the medium into which the evaporative outflows from closein expoplanets expand.
A closein planet and its host stars could thus interact through their respective outflows and/or magnetic fields. In either case, the interaction could have a potentially observable signature. For example, it has been suggested that there might be a detectable bow shock formed by the supersonic stellar wind in front of the planet (e.g., Vidotto et al. 2010, 2011a, 2011b, 2014; Llama et al. 2011, 2013; BenJaffel & Ballester 2013), or an observable cometarytype tail associated with the sweptup planetary wind that trails the planet in its orbit (e.g., Schneider et al. 1998; Mura et al. 2011; Kulow et al. 2014). According to these proposals, a bow shock could potentially contribute to the absorption of stellar radiation before an ingress of a transiting Hot Jupiter, whereas a tail might contribute after an egress. It was furthermore proposed that the interaction with the stellar magnetosphere might lead to phased hot spots on the stellar surface or to flaring activity (Shkolnik et al. 2003, 2004, 2005, 2008; Walker et al. 2008; Lanza 2008, 2009; Cohen et al. 2009; Pillitteri et al. 2010, 2011, 2014), and that the planetary magnetic field could be affected in both its topology and its strength (e.g., Laine et al. 2008; Cohen et al. 2011b, 2014; Strugarek et al. 2014). It was also suggested (Cohen et al. 2010) that the interaction with a magnetized Hot Jupiter could lead to a reduction in the mass loss rate, and even more strongly in the angular momentum loss rate, from the host star.
In this paper we perform 3D magnetohydrodynamic (MHD) simulations that aim to address the following questions:

How can one classify the different types of star–planet interactions that arise when a Hot Jupiter orbits its host star?

What are the dynamical features that develop in each case, and what are their expected observational signatures?

What is the impact of this interaction on the planetary and stellar outflows?
We classify such interactions by studying the evolution of a grid of numerical models. In particular, we investigate whether it is the planetary magnetosphere or the planetary outflow that intercepts the stellar plasma, whether the stellar outflow terminates in a bow shock, the circumstances for forming a stable tail along the planet’s orbit, and the conditions under which accretion streams of planetary material reach the stellar surface. We carry out a parametric study and present the criteria that can be used to classify a system’s morphology and dynamics based on its physical parameters. We also consider some general properties of the interaction regions that bear on their potential detectability.
Previous work by Cohen et al. (2009, 2011a) has identified and described some of the physical processes that may take place in a system consisting of a magnetized star and a magnetized Hot Jupiter. These studies were, however, focused on the interpretation of two specific systems. In addition, the planetary outflows were not explicitly included; they developed based on the boundary conditions applied on the surface of the planet. In contrast, we perform a systematic study to classify the types of star–planet interactions, incorporating consistent planetary winds that are based on detailed models.
The paper is structured as follows. Sect. 2 explains the methodology and presents the numerical models. Sect. 3 reports the results of the simulations and describes the dynamics. Sect. 4 classifies the types of interaction and considers some of their potential observational manifestations. Finally, Sect. 5 summarizes the conclusions of this work.
2 Numerical modeling
As discussed in MurrayClay et al. (2009, see also ), evaporative Hot Jupiter outflows need to be modeled as a fluid (using a hydrodynamic or an MHD formulation) rather than as a collection of individual particles that undergo Jeans escape, because the plasma typically remains collisional beyond the flow’s critical point. MurrayClay et al. (2009) considered an unmagnetized outflow from a Hot Jupiter exposed to a UV flux near the low end of the range of values that we model (see Sect. 2.6), and evaluated the collisional mean free path to proton–proton collisions in a hydrogen plasma of density and temperature K (for which the relevant cross section is ). They inferred a value (cm) that is smaller than the density scale height at the sonic point (which lies at a distance of planetary radii from the center of the planet). The large value of this ratio indicates that the collisional approximation should be applicable to this problem over a wide range of circumstances. Similarly, our simple stellarwind model (Sect. 2.8) is also adequately modeled within the hydrodynamic framework (e.g., Parker 1960). The fact that both the planetary and the stellar outflows are likely magnetized sharpens these conclusions: for typical parameters, the proton Larmor radius (where is the speed of sound, is the proton mass, is the speed of light, is the electron charge, and is the magnetic field amplitude) is cm — much smaller than . Given that the region where the planetary and stellar outflows interact is located at a distance from the source of each of these winds that is comparable to the distance of its respective critical surface, an MHD treatment of this interaction should be amply justified.
A numerical simulation that encompasses both the star and the planet needs to adequately resolve the sizes of both objects. For the case of a Hot Jupiter orbiting a solartype host, the characteristic length scales of the planet and the star differ by roughly an order of magnitude. Since 3D numerical simulations are computationally demanding, the best strategy to model the star–planet interaction is to adopt multiple levels of refinement. One approach is to configure an adaptive mesh that follows the orbit of the planet and provides high resolution for its interior and environment. Another is to set up a static, multiply refined grid and combine it properly with a corotating frame of reference, so that the location of the planet stays fixed in the highly resolved region. Here we adopt the latter approach, which also simplifies the treatment of the interior of the planet. In general, a static, multiply refined grid is also appropriate for inclined orbits, as long as they are circular.
In the following, we denote the star with “” and the planet with “”. The simulations are performed in Cartesian coordinates , but for convenience we also use the notation of spherical coordinates . The planet is assumed to orbit in the  plane, or, equivalently, at . The center of the star is placed at the origin, , whereas the center of the planet at . For the initialization we use an extra set of coordinate systems denoted with primes, and , with their origins placed at (i.e., centered on the planet). At each spatial point we first calculate and then perform a coordinate transformation to obtain . In order to keep the center of the planet fixed, we adopt a corotating frame with , where is the angular velocity of the frame and is the Keplerian angular speed of the planet (with being the mass of the star and the gravitational constant). Since (where is the mass of the planet), we assume that and hence that the origin of the coordinate system is effectively located at the center of mass.
2.1 The MHD equations
The ideal MHD equations in a frame that corotates with the planet are:
(1) 
(2) 
(3) 
(4) 
The variables , , , and denote the density, pressure, velocity, and magnetic field, respectively, and is the polytropic index. The magnetic field is required to obey the solenoidal condition . Since we consider the nonrelativistic regime, , , and have the same values in the corotating frame and in the lab (star’s center of mass) frame. The velocity in the lab frame is obtained from the expression . The vector is the gravitational acceleration, and the inertial force that appears in the corotating frame, , has centrifugal
(5) 
and Coriolis
(6) 
components.
The temperature at the base of an EUVinduced Hot Jupiter outflow is found to be (e.g., MurrayClay et al. 2009): this value reflects the balance between photoionization heating and collisionally excited radiative cooling (e.g., Spitzer 1978). Furthermore, in cases where the stellar UV flux is sufficiently high, radiative cooling continues to roughly balance the photoionization heating as the wind moves away from the surface, with the result that the temperature changes little along the flow (e.g., MurrayClay et al. 2009). The isothermal approximation was also found to apply to the inner regions of the solar wind (e.g., Cranmer et al. 2007), and has thus been adopted in previous simulations of both stellar (e.g., Matt & Pudritz 2008) and planetary (e.g., Trammell et al. 2014) outflows. Based on these considerations, we model both the stellar and the planetary outflows in our simulations as being nearly isothermal throughout, characterizing them by a polytropic index as in the fiducial model of Matt & Pudritz (2008). This approach enables us to simulate the dynamics of these winds without having to deal with the intricacies of a detailed calculation of the thermal structure of the outflow. As we show in Sect. 2.8, the adopted framework can be made to yield density and velocity profiles for the planetary wind that are compatible with the results obtained by employing the more elaborate (albeit nonmagnetic) 1D model of MurrayClay et al. (2009), which incorporates photoionization heating and radiative as well as adiabatic cooling. Koskinen et al. (2013) — who adopted a more detailed modeling approach that included, among other factors, a consideration of the full spectrum of ionizing photons (rather than just a single representative frequency) — found quantitative differences from the MurrayClay et al. (2009) findings (such as a slower increase of the ionization fraction with distance from the planet’s surface; see also Trammell et al. 2011), but their results remain qualitatively similar to those of MurrayClay et al. (2009). We therefore consider the adopted approximation to be justified for the present study, especially in view of the fact that the classification criteria presented in Sect. 4 do not depend on the details of the wind models.
We further assume, for simplicity, that the plasma consists of pure atomic hydrogen, and that it is fully ionized, so that the temperature is given by (where the factor represents the mean molecular weight, and is the Boltzmann constant). The strong ionization and nearisothermality assumptions for the planetary wind are consistent with the “high UV flux” case (corresponding to a young host star; e.g., Ribas et al. 2005) presented in MurrayClay et al. (2009). By contrast, in the “low UV flux” case considered by these authors (corresponding to a solaranalog host), the wind temperature drops with distance from the stellar surface on account of adiabatic cooling. To mimic this case, we also consider models with a characteristic temperature that is lower than .
At locations where the magnetized stellar wind interacts with the planetary outflow or magnetosphere, fieldline reconnection may potentially take place. This process can be realized in the simulations — even though resistivity is not explicitly included in the formulation — because numerical diffusion allows the magnetic field to reconnect in the presence of strong current sheets. We note, however, that the model setup adopted in this paper represents the star and the planet as aligned magnetic dipoles, a configuration that does not give rise to Xpoint structures and is therefore not prone to reconnection.
2.2 Stellar and planetary parameters
Parameter  Symbol  Stellar value  Planetary value  Units 

Radius  ,  –  
Mass  ,  –  
Temperature at the base of the outflow  ,  –  K  
Density at the base of the outflow  ,  –  g cm  
Equatorial surface magnetic field  ,  –  G  
Escape speed at the base of the outflow  ,  –  km s  
Sound speed at the base of the outflow  ,  –  km s  
Escape velocity parameter  ,  –  –  
Plasma at the base of the outflow  ,  –  –  
Rotation period  ,  –  days  
Rotation speed parameter  ,  –  –  
Lagrange point 1  ,  –  –  
Stellar UV flux  –  –  erg cm s  
Orbital radius  –  –  AU  
Orbital period  –  –  days  
Orbital speed  –  –  km s 
Table 1 lists the basic model parameters that characterize the host star and the Jupiterlike closein planet. The masses (, ) enter only in the gravitational field at the exterior of each body, whereas their radii (, ) define the spheres within which an internal boundary is applied. The temperatures correspond to the values at the base of the respective outflows (the stellar corona and the UV absorption layer in the outer planetary atmosphere). The value of the density of the stellar outflow, , is chosen so that the massloss rate be comparable to that of the solar wind (a few times ), an approach adopted from Matt & Pudritz (2008). For the planet, the base density value is determined from the requirement that the simulated wind match highresolution 1D numerical simulations that include the relevant heating and cooling processes (see Sect. 2.8). Thus, should not be regarded as the actual density at optical depth (the physical base of the outflow) but rather as a numerically motivated value that gives the appropriate wind profile. The initial magnetic fields of the star and the planet are assumed to be dipolar, with equatorial surface values equal to and , respectively.
The escape speed from the stellar surface is given by , whereas the sound speed at the base of the stellar wind is . A useful model parameter for determining the properties of the wind is . Another important nondimensional variable is the ratio of the thermal and magnetic pressures at the stellar surface, (the stellar plasma beta), which characterizes the dynamical significance of the magnetic field.^{1}^{1}1For plasma values less than 1, the magnetic field dominates the dynamics and constrains the plasma to move along flux tubes. Conversely, when the value of this parameter exceeds 1, the fluid dominates the dynamics and the magnetic field is dragged by the flow. The third relevant parameter for characterizing the outflow is , which measures the effect of surface rotation. Corresponding expressions for , , , , and are used for the planet.
The stellar UV flux, , is a critical quantity because it determines the energy input and hence the strength of the planetary wind. Finally, Table 1 gives the rotational periods of the star () and the planet () as well as the orbital characteristics, , , and . The host star is taken to rotate at of the breakup angular speed (as compared to of for the present sun), whereas the planet is assumed to be tidally locked (i.e., ), which is a reasonable approximation for a closein planet that interacts tidally with its host star (e.g., Matsumura et al. 2010). Our model planets thus rotate at – of breakup.
2.3 Stellar and planetary winds
We base the numerical treatment of the two outflows on the simplified approach of Matt & Balick (2004). In brief, we initialize dipolar magnetospheres (forcefree by definition) and keep fixed at each surface the physical quantities that drive the flow. The temporal evolution of such configurations opens up the magnetospheres and leads to steadystate winds (see Matt & Pudritz 2008). Although this method does not require the winds to be specified explicitly, we have opted to initialize a simple form of outflow to ensure that the starting wind profiles are in agreement with detailed, selfconsistent models in the literature. We emphasize that, because of the transonic nature of these flows, the final steady state does not depend on the initial configurations, only on the boundary conditions that are imposed at the surfaces of the two objects. However, in lowresolution simulations that do not include all the relevant physics, simple assumptions about the density and temperature at the boundary might not give correct results. Therefore, we use the profiles of detailed outflow models as a guide in order to set up the appropriate boundary conditions. This procedure is particularly useful for capturing the correct mass loss rate of the planet, as described in more detail in Sect. 2.8.
We initialize the velocity field of each outflow at time as an isotropic, isothermal Parker wind (Parker 1958). In particular, the velocity profile for the stellar outflow, , is obtained by solving numerically the equation:
(7) 
where , , and is defined in Sect. 2.2. In this equation we take the speed of sound to be strictly isothermal, i.e., (corresponding to a polytropic index that is exactly 1). Since the star is assumed to be rotating, we make the approximation that the initial wind is rotating rigidly with it. This is a minor contribution to the outflow velocity because is an order of magnitude smaller than (and, in any case, the flow selfconsistently assumes the correct values of as the simulation evolves). The Cartesian components of the initial stellar outflow velocity are then given by the sum of the wind, rotation, and frame speeds:
(8) 
(9) 
(10) 
To derive the corresponding expressions for the planetary outflow, we first calculate from Eq. (7), using , , and . We then add the planet’s orbital velocity, , as well as the frame and planetaryrotation velocities (using the rigidbody assumption), to obtain:
(11) 
(12) 
(13) 
where , , and are functions of , as explained at the beginning of this section. For simplicity, we have assumed that the initial planetary wind orbits and spins rigidly with the planet, an approximation that is appropriate only close to the surface. Note that, for a tidally locked planet, the orbital velocity cancels out with the frame and planetary rotation velocities, and the planet stays fixed at the location . For example, the components of these three terms at a distance along the line that connects the centers of the star and the planet are . Under this assumption, given that we are working in the corotating frame, we could have just as well initialized the velocity by using only the thermal (Parker) wind, with no extra terms.
To obtain the pressure profile of the stellar wind, we solve analytically the hydrodynamic radial momentum equation,
(14) 
which gives:
(15) 
The density is then found using the isothermalplasma assumption:
(16) 
In a similar vein, we get for the planet
(17) 
and
(18) 
We set up the stellar wind everywhere in the computational box, apart from a sphere of radius , centered on the planet. Within this volume we initialize the planetary wind. Both outflows are initially supersonic and superAlfvénic (with the Alfvén speed given by ) at the interface that separates them.^{2}^{2}2Recall that in MHD there are three types of waves (listed in order of increasing propagation speed): slowmagnetosonic (SMS), Alfvén, and fastmagnetosonic (FMS). Along field lines, an FMS wave propagates at and an SMS one at . Across field lines, the SMS and Alfvén wave propagation speeds are zero, and an FMS wave propagates at . The FMS critical surface is the separatrix beyond which no perturbation can propagate oppositely to the flow (analogous to the sonic critical surface in hydrodynamics); we henceforth refer to this surface as the “critical surface”. We also follow a common practice in the literature and plot the Alfvén and sonic separatrices rather than the SMS, Alfvén, and FMS ones.
It is worth noting that the Parker wind (Eq. 7) is nondimensional (with the radial distance expressed in units of the radius of the object and the velocity in units of the sound speed at the base) and that its solution depends only on the parameter , which can be written in the form
(19) 
Interestingly, for a solar analog (, , ) and a Hot Jupiter (, , ) we obtain . As a result, the two outflows reach their corresponding sonic speeds at the same scaled distance. For the adopted fiducial parameters, , which implies that the stellar and planetary outflows become supersonic at a radius that is slightly smaller than and , respectively.
In a real system, the stellar flux that heats the planetary surface varies with longitude (in particular, there is no direct irradiation of the planet’s night side) as well as latitude (in particular, the planetary poles also receive no direct heating). This implies that even an unmagnetized outflow would be anisotropic, although the quantitative effect of the uneven irradiation on the evaporative mass outflow rate need not be large (e.g., MurrayClay et al. 2009). In this work we simplify the treatment by assuming that the base temperature and density are uniform over the entire planetary surface.
The initial magnetic field configuration in the computational domain is taken to be a sum of two magnetic dipoles, , which are given by
(20) 
and
(21) 
for the star and planet, respectively. Both dipoles are oriented along the axis and are aligned with each other. We also write explicitly the components of the total gravitational field outside the two bodies:
(22) 
(23) 
2.4 Wind boundary conditions
In order to continuously accelerate the outflows and guarantee a constant supply of mass, we keep fixed the density, pressure, and velocity at their bases. Specifically, during the temporal evolution, we impose the initial values of these quantities at the regions defined by and for the stellar and planetary winds, respectively. However, the magnetic field is free to evolve within these shells, readjusting selfconsistently from its initial dipolar topology. This is a simplified treatment as compared with the formulation of Matt & Pudritz (2008), who modeled 2.5D axisymmetric stellar winds using a thinner, fourlayer shell above the surface of the star.^{3}^{3}3A 2.5D MHD simulation refers to the evolution of all three components of the velocity and magnetic field in a 2D computational domain. Within that region, they progressively relaxed the timedependency constraint on the physical quantities. Here we employ a less detailed setup on account of the lower resolution imposed by our 3D modeling. Nevertheless, the steadystate winds that we obtain compare well with the more refined 2.5D simulations (see Sect. 2.8).
The above implementation of stellar winds differs from the approach adopted by Cohen et al. (2011a, b, 2014) and Vidotto et al. (2014) (and references therein), which includes nonaxisymmetric and/or temporal variabilities. We neglect such effects in this work in order to investigate the star–planet interaction in a general manner that does not depend on temporary local fluctuations.
Our wind boundary conditions do not incorporate the effect of tidal forces on the planetary outflow. As discussed by Trammell et al. (2011), these forces could suppress the wind from the polar regions of a tidally locked planet with a dipolar field geometry. In particular, these authors showed that the wind would not undergo a sonic transition in this case if (where the parameters and are defined in Sect. 2.2). Our neglect of this suppression in the adopted boundary conditions is not a serious omission in view of the fact that our simulations incorporate the effect of these forces for and that the sonic surface typically lies at a distance of a few from the planet’s center.^{4}^{4}4 We note in this connection that the only two models in Table 2 that develop supercritical outflows, FvRb and FvRB, do not satisfy the inequality . This is consistent with the fact that the sonic surface plotted in Fig. 9 is represented by a closed curve. Furthermore, for the star–planet configurations considered in this paper, which are characterized by parallel orbital and spin axes, the bulk of the interaction between the stellar and the planetary plasmas occurs in the equatorial plane, with the details of the planetary outflow from the polar regions having little effect on the results.
2.5 Stellar and planetary interiors
The interiors of both the star and the planet are treated as an internal boundary. In particular, the values of all physical quantities are kept fixed by having them overwritten at each time step. To avoid spurious effects at the surfaces of the objects, or extreme dynamics that could affect the global time step of the simulation, we prescribe an approximate equilibrium within their volumes.
Since the physical conditions in the interior of the star or the planet are not important for the simulation, the simplest approach is to consider them as uniformdensity spheres. The interior gravitational accelerations are then given by
(24) 
and
(25) 
from which the pressures can be inferred using the hydrostatic equilibrium condition:
(26) 
(27) 
We keep the dipolar expression for the magnetic field, but we smooth it out at the center where it becomes singular. Specifically, we assume that the central regions of the star and the planet, and , are uniformly magnetized spheres, with magnetic fields and that connect smoothly with the dipolar configurations at larger radii. Finally, we approximate the star and the planet as rotating rigid bodies and specify their angular velocities by and , respectively.
2.6 The models
Model  Irradiation  

FvRB  5  High UV flux  
FvRb  High UV flux  
FvrB  High UV flux  
fvRB  Low UV flux  
fvRb  Low UV flux  
fvrB  Low UV flux  
FVRB  High UV flux  
FVRb  High UV flux  
FVrB  High UV flux  
fVRB  Low UV flux  
fVRb  Low UV flux  
fVrB  Low UV flux 
Table 2 lists the models that we study numerically. We aim to explore the observationally relevant regions of parameter space and attempt to deduce, based on the behaviors exhibited by the simulated systems, a general — and largely model independent — classification scheme for magnetized star–planet interactions. We usually adopt two representative values for each of the physical variables that need to be specified in the model. Thus, we consider two combinations of planetary parameters, the first having Jupiter’s mass and radius and the other corresponding to a less massive () and largerradius () planet. These choices are intended to sample the typical range of values for the escape speed from a hot Jupiter, which is relevant to the mass evaporation rate in the low limit (where ; e.g., MurrayClay et al. 2009).
We similarly consider two magnitudes for the UV flux, one low () and one high (), adopting numerical values similar to those given for these two limits in MurrayClay et al. (2009). These two values correspond to different evaporation regimes. In the highflux case, the excess photoionization energy is largely lost by radiative cooling (radiation/recombinationlimited evaporation), whereas in the lowflux case it is mostly absorbed and drives the outflow (energylimited evaporation). The outflow from a strongly irradiated planet is denser, hotter, and faster than that from a weakly irradiated one (e.g., MurrayClay et al. 2009). Since in this paper we do not explicitly model the radiative processes in the planetary atmosphere, we specify the wind properties through the values of the base temperature and density, which are adjusted to match the results of more detailed calculations (see Sect. 2.8).
We also vary the distance of the Hot Jupiter from the host star. Smaller orbital radii imply a stronger tidal force as well as a locally weaker, and possibly not fully developed, stellar wind.^{5}^{5}5We do not take into account changes to the incident UV flux that result from variations in the orbital radius since their magnitudes remain small in comparison with the difference between our chosen representative values. The value of the orbital radius also determines whether the planet lies within the critical surface of the stellar wind or whether it is impacted by a supercritical flow.
The two values that we adopt for the magnetic field amplitude at the planet’s surface differ by an order of magnitude and are intended to represent a “strong” and a “weak” field. Stronger fields correspond to higher magnetic pressure and tension and exhibit a greater “rigidity.” In this case the field resists being opened up by the stellar outflow over a larger region near the equator, and as a result the planetary magnetosphere presents a larger obstacle in its interaction with the stellar wind.
Although we vary the stellar UV flux, our simulations employ only a single set of stellar wind parameters. Our coverage of the relevant parameter phase space is, however, sufficiently broad to yield general results, which in principle apply also to systems with either a weaker or a stronger stellar wind.
2.7 Numerical setup
The simulations are performed with PLUTO (version 4.0.1), a code for computational astrophysics (Mignone et al. 2007, 2012).^{6}^{6}6PLUTO is freely available at http://plutocode.ph.unito.it The computational domain consists of a cube, with , which is resolved by a static, multiply refined grid of . We have also carried out one simulation (model FvRb) in higher resolution, , in order to validate our results. The region around the star, , is resolved with ( cells) in most cases, and with ( cells) in the highresolution simulation. However, for the Hot Jupiter and its surroundings, , we use a resolution that is higher by one order of magnitude, i.e., ( cells) and ( cells) for the standard and highly resolved cases, respectively. In the region between the star and the planet we apply the resolution .
The total time of the simulation is typically days, only a fraction of which is required for the system to attain a quasisteady state. For comparison, a stellar wind propagating at travels times the distance from the star to the outer boundary over this time interval. We adopt highly accurate numerical schemes to compensate for the resolution limits imposed by the threedimensional character of the simulations. In particular, we use the HLLD Riemmann solver and apply thirdorder accuracy in space (piecewise parabolic interpolation) and time ( order RungeKutta). The condition is implemented using the divergencecleaning method, an approach based on the generalized Lagrange multiplier (GLM) formulation (see Mignone et al. 2012 and references therein).
2.8 Verification and calibration of 3D wind models
The acceleration and final steadystate properties of a simulated wind depend on the included physics as well as on the resolution of the grid. 3D simulations are computationally expensive and therefore cannot incorporate the same level of detail that is possible in 1D and 2D models. In this work we have implemented a simplified procedure for treating the stellar and planetary outflows, and we now check on the validity of the adopted approximations. We first compare our 3D stellarwind model with a more detailed (and higherresolution) 2.5D model and verify its consistency. We then demonstrate our method of ensuring the consistency of 3D planetary wind models with 1D hydrodynamic calculations that incorporate the relevant physics of evaporative outflows.
Figure 1 compares the 3D stellar wind model adopted in this work (left panel), with a 2.5D simulation with the same parameters (right panel), but with a times higher resolution along each direction (modeled as in Matt & Pudritz 2008). The density, velocity, and magnetic field have very similar profiles, despite the fact that every cell on the left would be resolved by a thousand cells on the right if it were a 3D setup. Any discrepancies can be neglected since the outflow serves its purpose of providing a generic stellar wind.
Figure 2 shows the profiles of the density (left panel) and radial velocity (right panel) for the steady state reached by our stellar wind model. The radial velocity is zero along the equator up to , corresponding to the extent of the stellar dead zone (see Fig. 1). More generally, the finding that the wind speed remains lower along than along for any given value of can be attributed to the fact that a smaller fraction of the field lines open up near the equator than near the poles. This behavior is also exhibited by more detailed models of winds from solartype stars (e.g., Cohen & Drake 2014). In the case of faster rotators, in which magnetocentrifugal effects play a more prominent role in the acceleration of the wind, the outflow speed along the equator can become higher than along the spin axis (e.g., Matt & Pudritz 2008). This regime is pertinent to our models of tidally locked planets.
To ensure that the planetary winds are properly initialized, we specify the boundary conditions so that the resulting profiles match the results of a more elaborate (and highly resolved) 1D model that incorporates relevant additional physics. Specifically, we simulate the outflow along the line joining the centers of the planet and the star, incorporating explicit energy and ionization balance equations following the formulation of MurrayClay et al. (2009). The energy equation includes advective and work terms as well as source terms in the form of photoionization heating (assuming a single ionizing photon energy of eV and perfect efficiency) and Ly cooling (under the assumption of collisional excitation by free electrons). This equation is integrated using a modified version of the Simplified NonEquilibrium Cooling (SNEq) module of PLUTO (Teşileanu et al. 2008). The hydrogen ionization balance equation accounts for photoionization, Case B radiative recombination, and ion advection. Using the results of these calculations, we finetuned the boundary conditions of the simplified 3D models until the latter were found to effectively capture the general features of the 1D models for our representative planets in both the high and lowflux limits. The steadystate planetary wind profiles obtained in this way are shown in Figs. 3 and 4 together with the corresponding 1D results.
As can be seen from Figs. 3 and 4, the 1D outflow models predict a very steep density drop near the surface of the planet. This implies that, if one were to choose the value of the surface density for the 3D simulation simply on the basis of the location of the nominal base of the flow (where ), this would likely not lead to a good match with the 1D density profile. On the other hand, if one were to choose even a slightly different reference radius for fixing the surface density, this might result in a significant under or overestimate of the mass outflow rate. These considerations strongly suggest that the correct procedure for modeling the underlying planetary outflow is to choose the boundary values of density and temperature on the surface of the planet so that they lead to a good match with the reference profiles. The small discrepancy seen in the velocity profile of the lowflux case in Fig. 4 can be ignored because, as we show below, the interaction with the stellar plasma chokes the outflow from a Jupiterlike planet for this value of the flux.
3 Numerical results
3.1 Unmagnetized star–planet interactions
Before turning to the MHD simulations, we make a few qualitative remarks on the expected behavior in the absence of magnetic fields. Assuming that there is no planetary outflow, the stellar wind will impact the planet directly. For a supersonic stellar flow, the height above the planet’s surface where the plasma will get shocked is determined by the location where the atmospheric thermal pressure, , is equal to the local value of the wind ram pressure, . For a subsonic flow, there will be no shock, but rather a smooth transition of the physical quantities between the two bodies: the density and pressure will have high values at the surface of each object and will be lower in between. The converse will hold for the velocity: zero values at the surfaces, and a smooth acceleration — and subsequent deceleration — on moving from the host to the Hot Jupiter.
If a planetary outflow is present, its strength will determine whether it will be suppressed by the stellar wind (weak outflow) or else become supersonic and then collide with the stellar plasma well away from the planet’s surface (strong outflow). When both flows are supersonic, they interact at the location where the ram pressure of the planetary wind, , is equal to that of the stellar gas (as measured in the corotating frame).
3.2 General behavior
We start our presentation of the MHD simulations by describing the general structure of the interacting outflows; in the ensuing subsections we present the results for each of the simulated models in greater detail.
We use model FVRB for this illustration, and show its initial configuration in Fig. 5; the initial structure of the other models is similar except for the planetary wind profiles, which depend on the choice of surface parameters (see Figs. 3 and 4). The star is located at the center of the left panel and of the top right panel, which respectively provide faceon and edgeon views of the orbital plane. The stellar wind is depicted with arrows, and its initially dipolar magnetic field with solid lines. Next to the star, at (in units of ), is the spherical volume within which the Hot Jupiter and its outflow are located (left and top right panels). A closeup of the region near the planet is shown in the bottom right panel. Note that both the stellar and the plentary winds are superAlfvénic and supersonic at the initial interface between them (a sphere of radius , centered on the planet).
At the beginning of the simulation, the plasma escaping from the surface of each object forces the magnetosphere to open. In general, depending on the strength of the magnetic field and the physical conditions at the base of the wind, this might happen over the entire surface (for a weaker field and/or a stronger outflow) or just at the polar regions (for a stronger field and/or a weaker outflow). For the adopted value of the stellar surface magnetic field (), the star forms a dead zone for . In this region, the hot plasma cannot overcome the opposing forces of stellar gravity and magnetic tension, both of which are almost perpendicular to the flow. The behavior of the planetary magnetosphere will be discussed on a case by case basis: its structure depends both on and the planetary wind properties and on its environment.
3.3 Models exhibiting a planetary tail
Figure 6 displays the final steady state of model fvRb. This case corresponds to a low UV flux, and the planetary outflow is weak and cannot overcome the ram pressure of the stellar wind (projected along the orbit). The stellar wind is stopped by the magnetic pressure of the planet’s closed field lines (which exceeds the thermal pressure above the planet’s surface), terminating in a bow shock that stands off a few planetary radii from the planet’s surface. The magnetosphere is dragged backward by the stellar wind, forming a cometarytype tail. This process is akin to the formation of the geomagnetic tail in the interaction of the solar wind with Earth’s magnetosphere (e.g., Ness 1965). The tail does not follow the trajectory of the planet since it is continuously pushed away radially by the stellar wind, resulting in a tilted structure.
Figure 7 illustrates how the structure of the tail is modified when the model parameters are changed. The left panel shows the results for model fvRB, which has the same parameters as model fvRb (Fig. 6) except that the planetary magnetic field is a factor of stronger, corresponding to an increase by a factor of in the magnetic pressure at a given distance from the planet. As expected, this results in the bow shock being located farther away from the planetary surface, which, in turn, gives rise to a wider tail. The right panel presents the final state for model FVRb, which corresponds to a higher UV flux (which promotes an outflow) but also a larger escape speed (which hinders a wind). Overall, its behavior is very similar to that of model fvRb, although the tail in this case is noticeably denser. This can be understood from the fact that, in the highflux limit, the base density scales roughy as (MurrayClay et al. 2009).
The snapshot shown in the right panel of Fig. 7 highlights a few interesting dynamical features of cometarytype tails. The plasma that trails the planet has a velocity comparable to the orbital speed, whereas the component of the stellar wind along the orbit is negligible at that location. This results in the development of strong shear, which may trigger a KelvinHelmholtz instability. Furthermore, the stellar plasma that pushes the tail outward is of lower density than the trailing stream, which can induce a RayleighTaylor instability. In this particular simulation, their effect on the tail structure remains comparatively mild, likely on account of the relatively high density of the trailing stream. We have, however, found that the induced distortions can be more persistent in other cases.
When the UV flux remains low but the escape speed is high, the planetary outflow weakens and the density in the tail becomes measurably lower than in the models considered so far. This result follows directly from the expression for the mass evaporation rate in the lowflux (energylimited) limit, which implies that (see Eq. 9 in MurrayClay et al. 2009). In this case the tail becomes thin and barely noticeable. This is demonstrated in Fig. 8 for model fVRb — the outflow is even weaker when the magnetic field amplitude is increased (model fVRB, not plotted).
3.4 Models exhibiting both a tail and an accretion stream
Figures 9 and 10 display the final quasisteady state of model FvRb. This model differs from the reference case of Sect. 3.3 (model fvRb, shown in Fig. 6) in having a high (rather than low) incident UV flux. This results in the formation of a strong outflow, which opens up the magnetosphere and becomes supersonic (bottom right panel).^{7}^{7}7The comparatively weak magnetic field in this case implies that the Alfvén critical surface is located closer to the planet than the sonic surface.
The planetary outflow propagates unperturbed over a few planetary radii, giving rise to the eyeshaped region around the planet at (left and top right panels of Fig. 9). Beyond that region, the planetary outflow collides with the stellar wind and forms a shock. This happens on both the day and night sides of the Hot Jupiter on account of the high orbital speed (which shifts the apex of the bow shock away from the substellar point) and the axisymmetry assumption for the planetary wind.
A noteworthy feature of the interaction in this case is the accretion of part of the shocked planetary outflow onto the host star. The infalling material spirals for a quarter of an orbit and then impacts the stellar surface. The fact that this gas falls in, rather than form a toruslike structure, can be attributed to the action of the shock and of the stellar wind, which slow down and disrupt the flow, and to the increase in the gravitational pull of the star as the gas spirals in. The rest of the outflowing planetary gas is, however, pushed back and forms a tail, as in the cases considered in the preceding subsection. Both the inspiraling and the trailing streams are affected by the velocity shear with the lowdensity stellar gas and by the outward acceleration that this component induces, which trigger KelvinHelmholtz and RayleighTaylor instabilities. These effects lead to the destruction of the tail in this case, with its fragments being pushed outward by the ram pressure of the stellar wind (left panel of Fig. 9). A very similar outcome is produced for model FvRB (not plotted), in which, however, the stronger planetary magnetic field results in the Alfvén surface lying farther away from the planet than the sonic one.
Accretion streams onto the stellar surface can form also in cases where the planetary outflow remains weak and the stellar wind is stopped by the planetary magnetic field. This is illustrated in Fig. 11, which shows the final configuration of model FVRB. This case is similar to model FVRb, shown in the right panel of Fig. 7, in that a combination of high UV flux and large escape speed generate a dense discharge from the planetary surface but only a weak outflow. In this case, however, the stronger magnetic field causes the interface between the planetary magnetosphere and the shocked stellar wind to lie at a larger distance from the planet’s surface; in particular, it now lies beyond the point. At that location, the gravitational pull from the star and the thermal pressure of the dead zone combine to deform the planetary magnetic field lines in such a way that accretion streams are formed. These streams, however, do not flow smoothly: as is seen in the left panel of Fig. 11, their interaction with the stellar wind and the magnetic field of the star causes their trajectories to split several times before they reach the stellar surface.
A similar situation characterizes our highflux “near” models, as illustrated in Fig. 12 with model FvrB. In this case, the planetary outflow is massive enough and the point is located close enough to the planet’s surface that both a planetary tail and an accretion stream are formed, resembling the behavior of model FVRB. However, the different position of the planet relative to the star modifies the evolution in this case. First, the fragmented tail is not pushed away by the stellar outflow, which is still weak at this location. Instead, the stellar outflow mixes with these fragments and slows them down, leading to their accretion by the star. Second, the denser ambient stellar gas in this model also enhances the mixing with the accretion stream, slowing the latter down and causing it to hit the stellar surface close to the orbital location of the planet. Finally, as the orbital motion of the planet is now faster than the stellar rotation, the stream does not attain a steadystate configuration. Instead, the magnetic field topology undergoes a continuous readjustment, with new magnetic accretion channels forming periodically. Model FVrB (not plotted) exhibits a similar morphology to that of model FvrB, with accretion flows onto the star originating from both the upstream and downstream regions of the planet, although the larger escape speed results in less mass leaving the planet in this case. We however find that there is no transfer of mass to the star when the UV flux is low, irrespective of whether the escape speed is large or small (models fVrB and fvrB, respectively; these are also not plotted).
4 Classification of star–planet interactions
Our simulations of the dynamical interaction between a magnetized, winddriving host star and a magnetized hot Jupiter that loses mass to photoevaporation were described in Sect. 3 in terms of the resulting morphological structures. We now attempt to distill from this phenomenology a general classification scheme that is based on the underlying dynamical processes operating in such systems. We quantify the relative influence of these processes with the help of characteristic length scales, and use the latter to identify four basic types of interaction. We then apply this scheme to the categorization of the models listed in Table 2.
4.1 Characteristic length scales
The formation of a quasistationary morphological structure in the interaction between a star and a closein planet entails the establishment of pressure equilibrium between the stellar and planetary plasmas at the interface separating these two media.^{8}^{8}8This interface formally constitutes either a contact or a tangential discontinuity. As confirmed by our simulations, such discontinuities may be subject to instabilities, which, among other effects, can induce mixing of the two plasmas. The relevant frame of reference for considering this equilibrium is that of the planet, and we label the pressure on the stellar side of the boundary in this frame by . The characteristic distance of this surface from the center of the planet, measured along (or close to) the line to the stellar center, is determined by whether it is the magnetic pressure, outflow ram pressure, or thermal pressure that dominates on the planet’s side of the boundary.
If the pressure of a closed planetary magnetosphere dominates, we can estimate the relevant scale () by assuming a dipolar field, , and equating its pressure to . This gives
(28) 
If, however, the planetary outflow is strong enough to become supercritical by the time it is stopped (in a shock) by colliding with the stellar gas, the relevant scale () is given, instead, by equating to the ram pressure . Expressing the density in terms of the planetary mass loss rate , we get
(29) 
Which of these pressure components (magnetic or ram) dominates the interaction with the stellar plasma is reflected in the relative magnitudes of the associated characteristic radii.
One could in principle also consider the contribution of the thermal pressure of the planetary atmosphere, but in practice its role is negligible because of its expected rapid decline with radius. In particular, under the isothermal approximation, (see Eq. 17), which represents a much faster drop than that of the magnetic pressure (). In fact, the magnetic pressure likely already dominates the thermal pressure at the (subsonic) base of the planetary outflow, where, using Jupiter’s mass, radius, and magnetic field (), and the base density () and temperature () inferred from our high–UVflux model (see Fig. 4), we estimate
(30) 
Note, however, that when the stellar wind is stopped by a supercritical planetary wind, the plasma on the planet’s side of the contact discontinuity is dominated by its thermal pressure (having passed through a shock).
The ambient pressure can be approximated by the sum of four contributions:
(31) 
which represent, respectively, the stellar plasma’s thermal, magnetic, and “intrinsic” ram pressure components, and the ram pressure associated with the relative motion between the planet and the ambient (stellar) gas. In this approximation, we take the direction of the stellar wind to be purely radial and the direction of the orbital motion to be purely azimuthal in the lab frame. If the stellar wind is supercritical, a bow shock will form at the interaction site. In the limit , the apex of this shock will lie along the line connecting the centers of the star and the planet, and the shock surface will be symmetric about this line. However, given that the orbital speeds of hot Jupiters are typically and are thus highly supersonic with respect to the ambient gas, such a shock will form even when the stellar wind is still subcritical when it is stopped. In this case, however, the shock will be oriented at a finite angle to the aforementioned line (see Eq. 37 below), and its apex will be displaced away from this line (in the direction of the planet’s motion).^{9}^{9}9Note that the speed of the orbital motion for these systems is of the order of the terminal speed of the stellar wind, so that, when the interaction occurs in the supercritical region of the stellar wind, the resulting bow shock will be significantly inclined with respect to the instantaneous “line of centers.” For typical parameters, we expect that the last two terms in Eq. (31) constitute the dominant contributions to .
One other characteristic length scale for this problem is the tidal (or Hill) radius,
(32) 
which gives the distance of the Lagrange point from the planet’s center. This distance is relevant to the question of whether the morphology of the planetary gas that interacts with the stellar plasma will be shaped only by the relative motion between these two media (through the effect of the ram pressure terms in Eq. 31) or whether the stellar gravitational field will also play a role. It can be expected that at least some of the planet’s plasma will be accreted onto the stellar surface if , but that essentially all of the planet’s gas will remain confined to the vicinity of its orbital radius if this inequality is not satisfied.
Before concluding this subsection, we note that the components of can be expressed in terms of the stellar parameters just as was done above for the planet’s pressure components. In particular, using , , and the isothermalgas assumption, we can express the terms appearing on the righthand side of Eq. (31) in the form
(33)  
(34)  
(35)  
(36) 
Thus, the values of the three characteristic length scales (, , and ) can be estimated from a given set of planetary (, , , , ), orbital (, ), and stellar (, , , , , ) parameters.
4.2 Types of interaction
We base our classification of the flow structures that arise from the hydromagnetic interaction between a closein giant planet and its host star on the relative magnitudes of the three characteristic length scales considered in Sect. 4.1: , , and . The four basic types of interaction that we identify in this way are sketched in Fig. 13 and discussed below. We note, however, that any given system may not fall exclusively under a single category. This could be due to the planet having an eccentric orbit, to the stellar wind being nonaxisymmetric, or to the presence of multipolar components in the stellar magnetic field, which could each lead to variations in the values of and on a timescale of days. A stellar magnetic cycle would be likely to induce variations on timescales of years in the properties of the stellar magnetic field, the stellar wind, and the stellar UV flux. Changes on even longer timescales could result from the effects of stellar evolution.
Type I:
Type I interactions occur when the planetary outflow is weak (corresponding to either a low irradiating flux or a large escape speed), so that the stellar plasma is intercepted by the planetary magnetic field (i.e., ). As we noted in Sect. 4.1, the relative motion between the planet and the ambient gas is typically high enough to lead to the formation of a bow shock upstream of the planet. The shocked stellar flow sweeps back the planetary magnetic field and diverts the planetary outflow in the downstream direction, leading to the formation of a thin planetary tail. The stronger the planetary field, the larger its pressure — hence the farther away from the planetary surface does the shock appear ( is larger), and the wider the tail that is produced. The tail gas remains confined to the vicinity of the planet’s orbit because .
The orientation of the tail is determined by the direction of the incident stellar wind as seen in the corotating frame at the apex of the bow shock. Approximating the wind velocity to be nearly radial in the lab frame and of magnitude (see Fig. 5), the angle that the tail makes with respect to the tangent to the planet’s orbit (indicated on the top left drawing in Fig. 13) can be approximated by
(37) 
(see also Vidotto et al. 2010). In the limit of an extended and static stellar corona (), the tail will trail the trajectory, forming a torus around the star. In the other limit, when a strong wind hits a comparatively distant planet, the tail will be almost perpendicular to the orbit, and so will lie nearly along the line of sight during transits. Note that we did not include the effect of radiation pressure, which is the dominant force in the formation of comettype tails in closein rocky planets (e.g., Rappaport et al. 2012, 2014). In the case of the gaseous giant planets considered in this paper, radiation pressure on Ly