Impact of Changing Stellar and Planetary Magnetic Fields on (Exo)planetary Environments and Atmospheric Mass Loss

A star interacts with the planets it hosts via magnetized stellar wind, which directly impacts their magnetospheric configuration and atmospheric loss. Gravitational (tidal) interactions can be neglected if the planet lies far away (as in the case of the Sun–Earth system), and in such a scenario the stellar wind independently affects the planetary environment without the possibility of the coronal structure of the star getting modified. The evolution of planetary atmospheres and the nature of its interaction with the stellar wind are profoundly affected by changes in their respective magnetic field strengths (Cohen et al. 2015; Basak & Nandy 2021). The discovery of several exoplanets (Mayor & Queloz 1995; Lunine et al. 2009; Pepe et al. 2014) has sparked interest in exploring the signatures of life-sustaining conditions and understanding habitability from the perspective of the origin and evolution of planetary atmospheres (Pollack & Yung 1980; Lammer et al. 2009; Lammer 2013; Cridland et al. 2017). Habitability—in the astrophysical context—depends on the ability of the stellar wind–forced planet to hold on to its atmosphere (Alvarado-Gómez et al. 2020) and some other aspects; this work focuses on stellar wind planetary interactions.

The age of the star determines its activity, which dictates the properties of the outflowing wind. Observations confirm the varying magnetic activity of the Sun during its evolutionary phases (Nandy & Martens 2007; Vidotto 2021), as well as variations in the field strength of solar system planets and stars outside our solar system (Stevenson 1983; do Nascimento et al. 2016; Kiefer et al. 2017). The exploration of magnetic field effects on exoplanets is a research pursuit of high topical interest (Grießmeier 2015; Oklopčić et al. 2020). With the detection of diverse magnetic activity in a number of star–planet systems, the question that naturally arises is—what are the consequences of magnetic interactions of (exo)planets with their host star? This is a critical question because the answer to this determines whether a planet can retain its atmosphere during its long-term interaction with the harboring star and, consequently, its habitability timescale. We explore how the nature of this interaction is determined by the relative strengths of the stellar wind magnetic field and planetary magnetosphere. Several efforts have been made to understand how stellar winds influence the planetary magnetosphere as well as its atmosphere (Lammer et al. 2012; See et al. 2014; Vidotto & Cleary 2020; Harbach et al. 2021; Nandy et al. 2021). Variation in stellar magnetic activity (Nandy 2004; Nandy & Martens 2007; Brun et al. 2014; Nandy et al. 2021; Tripathi et al. 2021; Vidotto 2021) introduces variations in stellar radiation (Spina et al. 2020), stellar wind speed (Finley et al. 2018), the magnetic field strength of plasma winds (Vidotto et al. 2015), and may lead to magnetic storms as well. As a consequence of these interconnected phenomena, the planetary dipolar field is deformed, resulting in a magnetospheric structure that may vary from what we often observe on the Earth (Gallet et al. 2017).

In this paper, we perform a thorough parameter space study using numerical modeling for understanding interactions in star–planet systems and the effect of stellar activity evolution on planets with different magnetospheric strengths. The magnetized stellar wind and the planet's inherent magnetic field are projected to have a significant impact on planetary habitability and atmospheric mass loss (Nandy 2004; Khodachenko et al. 2008; Nandy et al. 2017). We carry out three-dimensional global magnetohydrodynamic (MHD) simulations using a broad range of stellar and planetary magnetic field values for providing a comprehensive picture of the interaction process. Most extrasolar giant planets that have been discovered so far are anticipated to have significant ionospheres. To simplify the global simulation, and based on prior studies (Koskinen et al. 2010; Gronoff et al. 2011; Yan et al. 2019), we consider the planet to be surrounded by a perfectly conducting plasma atmosphere. We study the impact of an increasing stellar wind magnetic field on the magnetopause standoff distance and explore the conditions that lead to the magnetic pileup in the dayside region. Another significant occurrence we examine is the formation of Alfvén wings, where the magnetotail opens up when the upstream wind Alfvénic Mach number is low and the current sheet length shortens and bifurcates in the nightside region of the planet. We find a criterion for current sheet bifurcation that is consistent with all cases of our study and establish an analytical relation for the variation of the mass-loss rate with changing stellar and planetary magnetic fields.

This study is relevant for any star–planet system and those systems in which the magnetic fields of either or both entities vary considerably with time. Our findings are significant for a better understanding of (exo)planetary atmospheres and determining the impact of varying magnetic field strength on habitability. This paper begins with the model description in Section 2, in which we give an overview of the theory and numerical setup employed in this study. The detailed findings are presented in Section 3, followed by the conclusion in Section 4.

We adapt the Star–Planet Interaction Module (CESSI-SPIM) developed by Das et al. (2019) for simulating different configurations of stellar wind and planetary magnetospheres. The governing set of resistive MHD equations are given by:

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{v}}+({\boldsymbol{v}}\cdot {\rm{\nabla }}){\boldsymbol{v}}+\displaystyle \frac{1}{4\pi \rho }{\boldsymbol{B}}\times ({\rm{\nabla }}\times {\boldsymbol{B}})+\displaystyle \frac{1}{\rho }{\rm{\nabla }}P={\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\partial }_{t}E+{\rm{\nabla }}\cdot [(E+P){\boldsymbol{v}}-{\boldsymbol{B}}({\boldsymbol{v}}\cdot {\boldsymbol{B}})+(\eta \cdot {\boldsymbol{J}})\times {\boldsymbol{B}}]=\rho {\boldsymbol{v}}\cdot {\boldsymbol{g}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}+{\rm{\nabla }}\times ({\boldsymbol{B}}\times {\boldsymbol{v}})+{\rm{\nabla }}\times (\eta \cdot {\boldsymbol{J}})=0,\end{eqnarray} \tag{ 4 }$$

where the symbols ρ, v, B , P, E, and g denote density, velocity, magnetic field, pressure, total energy density, and gravitational acceleration due to the planet, respectively. J is the current density given by ∇ × B , ignoring the displacement current. The expression for the total energy density is given by

$$\begin{eqnarray}&&E=\displaystyle \frac{P}{\gamma -1}+\displaystyle \frac{\rho {v}^{2}}{2}+\displaystyle \frac{{B}^{2}}{8\pi },\end{eqnarray} \tag{ 5 }$$

for an ideal gas equation of state.

The computational domain extends from −80 R_p to 200 R_p in the x-direction, −45 R_p to 45 R_p in the y-direction, and −200 R_p to 200 R_p in the z-direction, where R_p is the radius of the planet, which is located at the origin of the Cartesian box. The region extending from −2 R_p to 2 R_p is resolved by 12 grids in all three directions, i.e., 1 R_p is resolved using three grids in the planetary vicinity. Keeping in mind our modest computational facility, a combination of stretched and uniform grid types is used in the x- and z-directions. In the x-direction, the regions extending from −80 R_p to −2 R_p and 2 R_p to 10 R_p are resolved using 156 and 16 grids, respectively, i.e., 1 R_p is resolved using two grids. Grids with a stretching ratio of 1.005175 are used from 10 R_p to 200 R_p . In the y-direction, the regions from −45 R_p to −2 R_p and from 2 R_p to 45 R_p are resolved by 86 grids each, i.e., 1 R_p is resolved using two grids. In the z-direction, the regions from −22 R_p to −2 R_p and 2 R_p to 22 R_p are resolved by 40 grids each, i.e., 1 R_p is resolved using two grids, whereas grids with a stretching ratio of 1.004958 are used for the regions extending from −200 R_p to −22 R_p and 22 R_p to 200 R_p .

In this study, an Earth-like planet is considered, with similar mass, radius, and tilt, while the stellar wind is assumed to have a speed of 270 km s⁻¹, corresponding to quiet times. We perform a set of simulations by varying the magnetic field strength of both the planetary magnetosphere (0.1B_e , 0.5B_e , 1B_e , 2B_e , and 5B_e ) and stellar wind (−0.1 nT, −0.5 nT, −1 nT, −2 nT, −5 nT, −10 nT, −30 nT, −50 nT, and −75 nT) with all possible combinations. Here, B_e = 3.1 × 10⁴ nT denotes the surface equatorial magnetic field strength of the Earth's dipolar field.

In order to study the effect of stellar activity evolution on planetary atmospheric loss, which directly impacts habitability, we initialize a conducting plasma atmosphere surrounding the planet. The density profile in the vicinity of the planet is defined by

$$\begin{eqnarray*}&&{\rho }_{\mathrm{pl}}={10}^{6}{\rho }_{\mathrm{atm}}\qquad r\leqslant {R}_{p}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{atm}}(r) & ={\rho }_{\mathrm{pl}}+\displaystyle \frac{({\rho }_{\mathrm{amb}}-{\rho }_{\mathrm{pl}})}{2} & \left[\tanh \left\{9\left(\displaystyle \frac{r}{{R}_{P}}-2\right)\right\}+1\right]\\ & & {R}_{P}\leqslant r\leqslant 3{R}_{p},\end{array}\end{eqnarray} \tag{ 6 }$$

where ρ_pl and ρ_amb are densities of the planet and ambient medium, respectively, and r is the radial distance from the origin. For more details on the justifications for the choice of the above atmospheric profile, interested readers may refer to Section 2.1 of Basak & Nandy (2021). The pressure distribution in the atmosphere is evaluated by numerical integration of the equation

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dr}}=-{\rho }_{\mathrm{atm}}(r)g(r).\end{eqnarray} \tag{ 7 }$$

Here, the gravitational field is given by $g(r)=-\tfrac{{{GM}}_{\mathrm{pl}}}{{r}^{2}}$, where M_pl is the planetary mass. The pressure inside the planet is evaluated by extrapolating the value of the pressure at the planet–ionosphere boundary. The density and pressure in the region r > 3R_pl are initialized to be equal to those in the ambient medium (Table 1).

The input parameters at the stellar wind injection boundary (at x = –80 R_p in the yz plane) are obtained by solving the Rankine–Hugoniot magnetized jump conditions for the given shock velocity. A southward-oriented (southward interplanetary magnetic field or SIMF) stellar wind is injected with density ρ_sw = 4 ρ_amb. For all other boundary faces of the Cartesian box, a force-free outflow boundary condition is implemented. The magnetic fields of the wind and (or) planet are varied in each run and the steady-state configuration of the planetary magnetosphere and atmospheric loss are analyzed.

When the stellar wind impinges upon the planetary magnetosphere, a dynamic pressure balance is achieved at the dayside of the planet, resulting in the formation of a magnetopause that inhibits the penetration of stellar plasma and protects the inner atmosphere from erosion. The balance between thermal, magnetic, and dynamic pressures determines the shape and location of the magnetopause. An earlier study by Basak & Nandy (2021) considers the magnetopause standoff distance (R_mp ) to be the location of balance between the dynamic P_d (incoming wind) and magnetic pressure P_m (near the planet) curves. This approach works well if the stellar wind magnetic field (B_sw) is weak, since the magnetic pressure of the incoming wind is very small and may be neglected in the calculations. However, when the stellar magnetic field is relatively stronger, it accumulates outside the magnetopause and the magnetic pressure in the magnetosheath region cannot be ignored anymore. Our study shows that the point of intersection between P_m and P_d migrates outward nearer to the bowshock as ∣B_sw∣ is increased (Figure 1). For a strongly magnetized wind, it is found that the dynamic pressure P_d is sufficiently high, even after this intersection point, which indicates that stellar wind plasma is able to penetrate beyond this location, and therefore it cannot be regarded as the magnetopause standoff distance. Thus, we need a different approach for evaluating the magnetopause standoff distance for the case of strongly magnetized winds. It has been shown in previous studies (Schield 1969; Shue & Chao 2013; Lu et al. 2015) that at the magnetopause the sum of the magnetic and thermal pressures (denoted by P_th+m henceforth) on the planetary side balances the incoming stellar wind (see also Shue & Chao 2013). These considerations result in the relation

$$\begin{eqnarray}&&{\left[{P}_{{th}}+{P}_{m}\right]}_{\mathrm{planet}-\mathrm{side}}={\left[{P}_{{th}}+{P}_{m}\right]}_{\mathrm{star}-\mathrm{side}.}\end{eqnarray} \tag{ 8 }$$

The distance from the planet where the P_th+m curve peaks and its gradient becomes zero along the subsolar line is then considered to be the magnetopause standoff distance. The intersection of the P_th and P_d curves gives the location of the bowshock standoff distance (Liu et al. 1999; Ma et al. 2004; Holmberg et al. 2019). We cross-validate the estimates of the bowshock standoff distance by studying the spatial variations of the Alfvén Mach number; these are in good agreement with the results obtained from our pressure balance analysis.

Zoom In Zoom Out Reset image size

Figure 1 shows the pressure balance for obtaining the bowshock distance (R_bs ) and the magnetopause standoff distance (R_mp ) for the case of planetary magnetosphere B_p = B_e (B_e = 3.1 × 10⁴ nT) with varying stellar wind magnetic field. Note that only the SIMF is considered for this study. The total thermal pressure (P_th) is obtained as an output of the PLUTO code (Mignone et al. 2007), while the dynamic pressure of stellar wind (P_d = $\tfrac{{\rho }_{\mathrm{sw}}{v}^{2}}{2}$) and magnetic pressure (P_m = $\tfrac{{B}^{2}}{8\pi }$) are evaluated from evolving simulation variables. As the stellar wind magnetic field strength is increased, the magnetopause moves closer to the planet and the bowshock moves farther away. For a weakly magnetized stellar wind, as shown in Figures 1(a)–(c), the stellar field accumulation is not significant and the intersection of P_d and P_m gives the R_mp location as depicted in Basak & Nandy (2021). As ∣B_sw∣ increases further, it is evident that the stellar wind magnetic pressure builds up significantly in front of the planet (Figures 1(d)–(f)). With the increase in ∣B_sw∣, both the thermal and magnetic pressure increase in the magnetosheath region initially up to some distance. However, the enhancement in the magnetic pressure gets prominent gradually, due to the accumulation of the stellar wind magnetic field outside the magnetopause.

Just after the magnetopause, magnetic pressure drops to nearly zero, representing the dayside reconnection X-point. Thermal pressure rapidly increases near the magnetopause due to energy conversion in the magnetic reconnection region (see, e.g., Lu et al. 2015).

The pressure balance for different planetary magnetic field strengths (B_p ), keeping the stellar wind magnetic field constant, reveals that with decreasing B_p both the bowshock and magnetopause move toward the planet (Basak & Nandy 2021). This implies that for lower values of B_p , the stellar wind is able to penetrate closer to the planet, resulting in greater atmospheric loss. We find significant increments of the magnetic pressure in the magnetosheath region for the case of weaker planetary magnetic fields as compared to stronger-field cases. Thus, either strengthening the stellar wind magnetic field or weakening the planetary magnetic field leads to magnetic field accumulation at the dayside of the planet. Hence, the ratio of the intrinsic planetary magnetic field and stellar wind magnetic fields is an essential factor that determines the dynamics of the star–planet interaction. This is in keeping with expectations.

The interaction of plasma winds and the planetary magnetosphere results in a certain deformation in the planetary dipolar field. This deformation leads to a teardrop-shaped magnetospheric structure that we usually observe for the Earth. A series of simulation results of interactions of the planetary magnetosphere with the stellar wind magnetic field is shown in Figure 2. In each subplot, the magnetic field lines (white streamlines) are traced over the current density (colormap) in the y = 0 plane, with the stellar wind plasma flowing from left to right. In plots (a)–(c), the simulations show the formation of extended magnetotail structures for varying planetary magnetic fields B_p = 0.1 B_e , 1 B_e , and 5 B_e , respectively. The stellar wind magnetic field strength is fixed at ∣B_sw∣ = 2 nT for all cases.

Zoom In Zoom Out Reset image size

As ∣B_sw∣ is increased to 10 nT (panels (d)–(f)) for the same set of planetary magnetic field values, the magnetotail begins to open up, and on the further increase of ∣B_sw∣ to 50 nT (panels (g)–(i)), the lobes of the tail open up even further (Turnpenney et al. 2020). The Alfvén Mach number is an important control parameter that governs the behavior of plasma embedded in the magnetic field. With increasing stellar wind magnetic field strength, the upstream Alfvén Mach number decreases, which leads to the appearance of Alfvén wings across the planet (Ridley 2007). Figure 3 illustrates the sub-Alfvénic plasma interaction with the magnetized planet and the formation of three-dimensional Alfvén wings in the planetary magnetotail, as shown by the region enclosed by the Alfvénic Mach number M_A isosurface (M_A ≤ 0.3). When the stellar wind magnetic field is strong, field lines bend but do not strongly drape the planet. Chané et al. (2012) have shown observational evidence of the formation of Alfvén wings on the Earth during 2002 May 24 and 25. Several moons in our solar system, including the Lunar wake (Zhang et al. 2016), Ganymede, Io, Europa, and Enceladus, have been observed to form Alfvén wings. Saur et al. (2013) have found signatures of sub-Alfvénic plasma interactions in several star–(exo)planet systems. With a decreasing upstream Alfvén Mach number, the wings form at larger angles from the equatorial plane and the length of the current sheet becomes shorter, i.e., reconnection takes place closer to the planet in the magnetotail region, as shown in the column-wise panels of Figure 2. The row-wise panels of Figure 2 show that for a given stellar wind magnetic field strength, Alfvén wings are more prevalent in the case of weaker planetary magnetic fields. Here, we present only a few simulation results to show how the magnetotail opens up. A similar pattern is found for other simulations that involve all possible combinations of parameters (described in Table 1). Vernisse et al. (2017) have shown that interactions that involve either a weakly magnetized obstacle or sub-Alfvénic upstream plasma velocities or both can lead to the dominance of an Alfvén wing structure—and this is corroborated by our simulations.

To analyze the opening up of the magnetotail in greater detail, we locate the distance from the planet where the current sheet bifurcates in the nightside region. Figure 4(a) shows the current density plot in the xz plane (left) and the corresponding current density magnitude (J_mag) variation along the subsolar line (right). Since the current sheet has a finite thickness along the y- and z-directions, we consider a hypothetical rectangular box outside the planet's atmosphere encapsulating the current sheet and extending from −0.5 R_p to 0.5 R_p along both the y- and z-directions. Thereafter, the spatially averaged J_mag is calculated on yz slices of the cube and its variation along the x-direction is plotted. We find that the current sheet starts to bifurcate, typically when the current density magnitude drops off to about 45%–50% of its peak value. These are represented by line B and line A, respectively, in Figure 4(a). This result is robust, considering all different sets (Table 1) of our simulations. The current sheet, extending from the planet to its bifurcation point on the nightside region, is hereby referred to as the current sheet length. The variation of the current sheet length with planetary and stellar wind magnetic field is shown in Figure 4(b). As the planetary magnetic field B_p increases, the bifurcation point recedes away from the planet. For B_p = 2B_e , we do not see any bifurcation within the simulation domain until ∣B_sw∣ = 30 nT. For a stronger planetary magnetic field of B_p = 5B_e , bifurcation starts around ∣B_sw∣ = 50 nT onward.

Zoom In Zoom Out Reset image size

It is evident from Figure 4(b) that the current sheet bifurcation point shifts nearer to the planet when either the stellar wind magnetic field increases or the planetary magnetic field weakens.

We also investigate how the variations in stellar wind and planetary magnetic field strengths impact the mass-loss rate of planetary atmospheres. To compute the total mass-loss rate, we consider a cube of edge length 6.6 R_p (extending from –3.3 R_p to 3.3 R_p in all three directions), with its origin coinciding with the center of the planet, such that the whole planet including its atmosphere is enclosed within the cube. Mass-loss rates are computed across all six faces of the cube, which contribute to the total mass loss. The rate of atmospheric mass loss depends on the extent of the stellar wind penetration into the planetary atmosphere. This is essentially denoted by a point at the dayside—the reconnection X-point—beyond which planetary magnetic field lines merge with the stellar field. If the magnetic field strength of the stellar wind is increased gradually, keeping that of the planet fixed, the X-point is expected to shift closer to the planet, inducing greater atmospheric loss. However, the distance of the X-point from the planet does not decrease continually with increasing stellar field. The reconnection neutral point saturates at a certain level, and then onward the stellar field starts accumulating at the dayside (Figures 1(d)–(f)), similar to that of an imposed magnetosphere; at this stage, the mass-loss rate attains a limiting value. Figure 5(a) confirms the increase in atmospheric mass-loss rate with increasing stellar wind magnetic field strength. Figure 5(b) shows the plot of the X-point distance from the planetary surface with increasing stellar wind magnetic field for different values of the planetary field. This signifies that for a weak planetary field, the X-point distance saturates to a steady value at lower stellar magnetic fields, while for a stronger planetary field, the saturation occurs for higher values of the stellar wind field. Figure 5(a) shows that for a weaker planetary field, the loss rate increases initially and then saturates at lower values of the stellar wind field, while for the stronger planetary field, the loss rate requires higher values of the stellar wind field to attain saturation. The atmospheric loss rate is found to have a strong negative correlation with the X-point distance from the planet. The Pearson correlation coefficient values are −0.9800, −0.9683, −0.9753, and −0.9659 for planetary magnetic fields 0.5B_e , 1B_e , 2B_e , and 5B_e , respectively.

Zoom In Zoom Out Reset image size

Figure 5(a) also indicates that for a conducting plasma atmosphere, as considered in the present study, the mass-loss rate increases with the increasing strength of the planetary magnetosphere. This is because a stronger magnetosphere provides a larger interaction area with the ambient medium and therefore may lead to greater loss when it interacts with the stellar wind. Moreover, as the atmospheric matter in our study is not considered to be charge neutral, it is expected that the atmospheric plasma will take part in magnetic reconnections and be lost in the process. Earlier studies show that the atmospheric matter stretches out and escapes through polar caps and cusps during reconnections (Gunell et al. 2018; Egan et al. 2019; Sakata et al. 2020; Carolan et al. 2021), resulting in a larger atmospheric mass loss for a stronger magnetosphere. Gunell et al. (2018) show that an intrinsic magnetic field does not necessarily protect a planet against atmospheric escape and the escape rate can be higher even for a highly magnetized planet. Therefore, the role of the intrinsic magnetic field in the protection of planetary atmosphere is a complex question.

3.1. Semi-analytical Expression for Planetary Mass-loss Rates

Our study shows that the global dynamics of the planet, current sheet length, and atmospheric mass-loss rates are significantly influenced by the relative strength of the planetary and stellar wind magnetic fields. For different sets (Table 1) of our simulations, Figure 6 shows that the ratio of planetary to stellar wind magnetic fields is an important factor in governing atmospheric mass-loss rates. Based on theoretical considerations, we establish an analytical relationship between these two in order to ascertain the dependence of the mass-loss rate on the variation of planetary and stellar wind magnetic fields.

Zoom In Zoom Out Reset image size

The unperturbed planetary dipolar magnetic field falls as $\tfrac{1}{{r}^{3}}$, thus the pressure balance at the magnetopause can be written as

$$\begin{eqnarray}&&{p}_{\mathrm{dynsw}}+\displaystyle \frac{{B}_{\mathrm{sw}}^{2}}{2{\mu }_{0}}+{p}_{\mathrm{thsw}}=\displaystyle \frac{{k}_{m}^{2}{B}_{p}^{2}{R}_{p}^{6}}{2{\mu }_{0}{R}_{{mp}}^{6}}+{p}_{\mathrm{th}p},\end{eqnarray} \tag{ 9 }$$

where the dynamic pressure due to stellar wind, the thermal pressure due to stellar wind, and the thermal pressure due to planetary magnetosphere at the magnetopause are depicted by p_dynsw, p_thsw, and p_thp , respectively. The symbols B_sw and B_p represent the stellar wind and planetary magnetic field strengths, respectively, while R_mp is the magnetopause standoff distance in terms of the planetary radius R_p . The factor k_m is a measure of the compression of the dipolar magnetic field by the stellar wind at the magnetopause (Mead & Beard 1964; Tsyganenko & Sitnov 2005).

From Equation (9), the magnetopause standoff distance can be expressed as

$$\begin{eqnarray}&&\left(\displaystyle \frac{{R}_{{mp}}}{{R}_{p}}\right)={\left[\displaystyle \frac{{k}_{m}^{2}{\left(\tfrac{{B}_{p}}{{B}_{\mathrm{sw}}}\right)}^{2}}{1+\tfrac{C}{{B}_{p}^{2}}{\left(\tfrac{{B}_{p}}{{B}_{\mathrm{sw}}}\right)}^{2}}\right]}^{1/6},\end{eqnarray} \tag{ 10 }$$

where C = 2μ₀(p_dynsw + p_thsw − p_thp ). As pointed out in Section 3, the atmospheric mass-loss rate (${\dot{M}}$) is linearly anticorrelated with the dayside X-point distance (X_point) and can be expressed as ${\dot{M}}=-{{AX}}_{\mathrm{point}}+B$, where A and B are fitting parameters that need to be determined. Moreover, in the case of a SIMF, the magnetospheric standoff distance R_mp is approximately equal to the dayside reconnection point (X_point). Therefore, Equation (10) can be expressed as

$$\begin{eqnarray}&&\dot{M}\simeq -{{AR}}_{p}{\left[\displaystyle \frac{{\left(\tfrac{{B}_{p}}{{B}_{\mathrm{sw}}}\right)}^{2}}{1+\tfrac{C}{{B}_{p}^{2}}{\left(\tfrac{{B}_{p}}{{B}_{\mathrm{sw}}}\right)}^{2}}\right]}^{1/6}+B.\end{eqnarray} \tag{ 11 }$$

Equation (11) gives an analytical relationship between the mass-loss rate and the ratio of planetary and stellar wind magnetic field strengths.

Figure 6 shows that our modeled output of the atmospheric mass-loss rate fits reasonably well with the analytical expression given by 11. The value of R_square (goodness of fit) for the modeled fit is greater than 0.99 for all instances of planetary magnetic fields.

Despite computational limitations, we have carefully chosen a grid resolution for our simulations that balances accuracy and computational feasibility. While it is important to recognize that higher grid resolutions may result in quantitative variations in the computed atmospheric mass-loss rate and slight deviations in the positions of the magnetopause and magnetotail, the main essence of our study is robust to further increments in resolution. The qualitative global large-scale structures are interdependent on parameter values and the underlying physics remains valid.

In this study, we have used three-dimensional global MHD modeling to explore the effects of varying stellar and (exo)planetary magnetic field strengths on the steady-state magnetospheric configuration and atmospheric loss rates of planets. This study is important for understanding how the magnetic activity evolution of stars affects the environment and habitability of planets with different planetary magnetic field strengths. Our results are therefore applicable to a wide domain of far-out (exo)planetary systems.

Computing the pressure balance at the dayside of the planet illustrates that either strengthening the stellar wind magnetic field or weakening the planetary magnetosphere results in stellar field accumulation in front of the planet, similar to that of an imposed magnetosphere. In such a scenario, the magnetopause standoff distance cannot be obtained using the conventional procedure of balancing the dynamic pressure of the incoming wind and magnetic pressure close to the planet; a more general treatment is necessary, by considering the sum of both the magnetic and thermal pressures. It is found that for a particular strength of planetary magnetopshere, increasing the magnetic field of the stellar wind makes the bowshock shift away from the planet, while the magnetopause is formed closer, implying greater wind penetration.

We find that the relative strengths of the planetary and stellar wind magnetic fields play a critical role in determining the magnetic reconnection in the magnetotail region. The long extended planetary magnetotail, which exists for a moderately magnetized stellar wind or strong intrinsic planetary magnetic field, ceases to exist for cases when the stellar wind magnetic field becomes extremely strong (i.e., the Alfvénic Mach number decreases) or the planetary magnetic field becomes very weak. This leads to the formation of Alfvén wings in the nightside wake region. Our findings suggest that the magnetotail begins to open up and the point of bifurcation shifts closer to the planet when either the magnetic strength of the stellar wind is increased or the intrinsic planetary magnetic field is reduced. We observe that the magnetotail current sheet starts to bifurcate when the current density magnitude drops to about 45%–50% of its maximum value.

We also find that the atmospheric mass-loss rate increases with the increase in the stellar magnetic field, due to greater wind penetration close to the planet. The atmospheric loss rate saturates at lower values of stellar wind magnetic field for a weak planetary field. In the case of a stronger planetary field, however, a stronger stellar wind magnetic field is required for the saturation of the atmospheric loss rate. The variation of the dayside X-point distance from the planet with stellar wind magnetic field provides a thorough explanation for the saturation of atmospheric loss rates, since the dayside X-point distance from the planet and atmospheric loss rate are found to have a strong negative correlation. By modeling the impact of the planetary magnetic field on atmospheric escape processes, we corroborate the study by Gunell et al. (2018), that the escape rate can be higher for strongly magnetized planets.

We establish an analytical relationship between the mass-loss rate and the ratio of (exo)planetary and stellar wind magnetic fields. Our study shows that the global dynamics of the planet, current sheet length, and atmospheric mass-loss rate are significantly influenced by the relative magnetic field strengths of the star and the planet. The results of our numerical simulation indicate that the analytical expression derived in this study succinctly captures the modeled atmospheric mass-loss rates.

The detailed parameter space study presented in this paper will be helpful for identifying stellar and planetary conditions that lead to a particular magnetospheric configuration in solar and exoplanetary systems. The results also illustrate the impact of magnetic field variability on atmospheric loss rates, which is relevant for understanding the habitability of (exo)planets.

D.N. acknowledges useful discussions with team members of Commission E4 (Impact of Magnetic Activity on Solar and Stellar Environments) of the International Astronomical Union. The authors thank Souvik Roy, Soumyaranjan Dash, and Chitradeep Saha for useful discussions. S.G. acknowledges fellowship support (Award No.: 191620059727) from the University Grants Commission, Government of India. The development of the Star–Planet Interaction Module (CESSI-SPIM) and the simulation runs were carried out at the Center of Excellence in Space Sciences India (CESSI), which is funded by IISER Kolkata, Ministry of Education, Government of India.