Astrospheres of Planet-Hosting Cool Stars and Beyond ⋅ When Modeling Meets Observations

Abstract

Thanks to dedicated long-term missions like Voyager and GOES over the past 50 years, much insight has been gained on the activity of our Sun, the solar wind, its interaction with the interstellar medium, and, thus, about the formation, the evolution, and the structure of the heliosphere. Additionally, with the help of multi-wavelength observations by the Hubble Space Telescope, Kepler, and TESS, we not only were able to detect a variety of extrasolar planets and exomoons but also to study the characteristics of their host stars, and thus became aware that other stars drive bow shocks and astrospheres. Although features like, e.g., stellar winds, could not be measured directly, over the past years several techniques have been developed allowing us to indirectly derive properties like stellar mass-loss rates and stellar wind speeds, information that can be used as direct input to existing astrospheric modeling codes. In this review, the astrospheric modeling efforts of various stars will be presented. Starting with the heliosphere as a benchmark of astrospheric studies, investigating the paleo-heliospheric changes and the Balmer $\text{H}\upalpha$ projections to $1~\text{pc}$, we investigate the surroundings of cool and hot stars, but also of more exotic objects like neutron stars. While pulsar wind nebulae (PWNs) might be a source of high-energy galactic cosmic rays (GCRs), the astrospheric environments of cool and hot stars form a natural shield against GCRs. Their modulation within these astrospheres, and the possible impact of turbulence, are also addressed. This review shows that all of the presented modeling efforts are in excellent agreement with currently available observations.

1 Introduction

The space between stars is not empty. On the contrary, charged and neutral particles fill interstellar space. Moreover, complex structures like molecules and dust can be found, forming the interstellar medium (ISM). Due to the existence of the stellar wind (SW), a constant particle stream from the star, the ISM is pushed radially outwards, and a local bubble, the astrosphere, is built up. Stars possess magnetic fields that are transported into the stellar wind and are frozen into the SW (Alfvén 1942). Because of the constant streaming of the SW, the stellar magnetic field is carried outward into the astrosphere building up the astrospheric magnetic field (AMF). Thereby, both the size of the astrosphere and the strength of the AMF depend on the stellar activity. Since the AMF remains rooted at the rotating stellar photosphere, an Archimedean spiral, also known as the Parker spiral (Parker 1958; Owens and Forsyth 2013), is formed, and due to the different polarities of the two stellar hemispheres and the boundary surface between the two, a wavy astrospheric current sheet (ACS) is established.

Figure 1 shows a detailed sketch of a pure hydrodynamic large-scale astrospheric structure (see also Scherer et al. 2015). In the rest system of the star (orange dot), which is moving uniformly through the surrounding ISM, the ISM appears as a uniform flow (here from the left) whose velocity corresponds to the relative motion of the star and is often supersonic. Another supersonic plasma flow, the SW, is propagated radially outward. If both flows are supersonic, there is a shock to subsonic velocities for both the ISM and the SW: a possible Bow Shock (BS) for the ISM and the Termination Shock (TS) in case of the SW (solid red lines), respectively. The area between TS and BS, where both subsonic flows meet, is the astrosheath (AS). Since merging of both flows is not possible, the astropause (AP) is formed: a tangential discontinuity (blue line), at which pressure equilibrium prevails, and stellar and interstellar plasma remain separated from each other. The area between AP and BS is the outer astrosheath (OAS), while the region between AP and TS correspondingly is the inner astrosheath (IAS). When the interstellar flow is subsonic, the BS does not form. The situation, however, gets even more complicated when stellar and interstellar magnetic fields are being included. Then the propagation speed is no longer the sound speed but the fast-magnetosonic speed. Moreover, because the interstellar magnetic field is usually not aligned with the ISM inflow vector, the astrosphere will become asymmetric, and, thus, full 3D magneto-hydrodynamic (MHD) modeling is required. However, in the particular case that the ISM magnetic field is aligned with the inflow vector and the stellar wind magnetic field is negligible (or both fields are zero; HD case), a 2D axis-symmetric model in one half-plane is sufficient (e.g., Baranov and Zaitsev 1995).

Fig. 1

Large-scale structure of an astrosphere (HD). Sketched are the shocks (in red), the tangential discontinuity (in blue), the sonic lines (SL, in light blue), and the tangential contact discontinuity (TD) that starts at the triple point (T). The red shaded area visualises the Mach disk (MD). A reflected shock (RS) emanates from T and is reflected between the AP and TD. In addition, the Mach numbers of the different regions are given. Figure after Scherer et al. (2016b), reproduced with permission © ESO

The relative movement of the star and the ISM distinguishes a preferred direction: the upwind direction, pointing towards the inflowing ISM, and the downwind direction. In the case of HD modeling, the axis of symmetry is also the stagnation line (solid black line), which is the streamline on which the stagnation point lies. However, the BS, AP, and TS only occur in the upwind direction (forming the so-called nose). Thus, as most observations show, astrospheres can be assumed to be bullet-shaped, a result of the conservation of momentum:

$$ \rho _{\mathrm{sw}} v^{2} + P_{\mathrm{sw}} = \rho _{\mathrm{ISM}} v^{2}_{ \mathrm{ISM}} + P_{\mathrm{ISM}}, $$

(1)

where $\rho$, $v$, $P$ represent the density, velocity, and the thermal pressure of the SW and the ISM, respectively. In the downwind direction, the astrotail, another shock of the SW from supersonic to subsonic velocities, is forming the Mach Disk (MD). Starting from the contact point of TS and MD, a tangential discontinuity (TD) enters the tail parallel to the stagnation line. The contact point of the TS, MD, and TD is known as the triple point (T) (Courant and Friedrichs 1948). Due to the entropy condition, a reflected shock (RS) emanates from the triple point and is reflected between the AP and TD. In the upwind direction of the flanks, a sonic line (SL) is formed because the velocity tangential to the shock normal does not change during the shock transition.

The above picture, however, does not hold for MHD simulations. Due to the asymmetry, the stagnation point must not lie on the inflow line, the line parallel to the inflow vector passing through the star, and the above equation needs to be modified (see below). In addition, the shock structures have to be replaced by the MHD shocks. Nevertheless, the TS and AP remain, while the bow shock can become a bow wave, even if the sound speed is smaller than the inflow speed. The AP remains a tangential discontinuity that has to be distinguished from the contact discontinuity, an additional feature of MHD discontinuities (e.g., Scherer et al. 2015).

The existence of a BS in front of the heliosphere has recently been under scientific discussion (see, e.g., Zank et al. 2013; Scherer and Fichtner 2014; Schwadron et al. 2015). The newest study by Kotlarz et al. (2018), utilizing HD and MHD models, shows that the existence of an ISM magnetic field facilitates the existence of a bow shock around every astrophysical object, including our heliosphere.¹ The importance of neutrals and other species like pick-up ions is discussed in Sokół et al. (2022).

2 (Magneto-)Hydrodynamic Modeling

Many problems in space and astrophysics require numerical modeling, in particular environments well described by fluid dynamics. Over the past years, several codes solving such problems with hydrodynamics (HD) or magnetohydrodynamics (MHD) have been developed. To model such problems, numerical codes most often make use of (multi-) fluid (Euler) equations.² Setting the permeability $\mu =1$, the conservative form of the combination of continuity, momentum, and energy equations with the Euler-Maxwell equations is given by:

$$\begin{gathered} \frac{\partial}{\partial t} \begin{bmatrix} \rho _{j} \\ \rho _{j} \mathbf {v}_{j} + P_{1} \mathbf {F}_{{\mathrm{rad}}} \\ E_{j} + P_{2} E_{\mathrm{rad}} \\ \mathbf {B} \end{bmatrix} + \boldsymbol{\nabla} \cdot \begin{bmatrix} \rho _{j} \mathbf {v}_{j} \\ \rho _{j} \mathbf {v}_{j} \otimes \mathbf {v}_{j} + \left (p_{j} + \frac{B^{2}}{2}\right )\widetilde{I} + P_{3}\mathbf {F}_{{\mathrm{rad}}} - \mathbf {B}\otimes \mathbf {B} \\ \left (E_{j}+p_{j}+\frac{B^{2}}{2}\right )\mathbf {v}_{j} + P_{4} \mathbf {F}_{{ \mathrm{rad}}} - \left (\mathbf {B}\cdot \mathbf {v}_{j}\right )\mathbf {B} \\ \mathbf {v}_{j}\otimes \mathbf {B} - \mathbf {B}\otimes \mathbf {v}_{j} \end{bmatrix} = \\ \begin{bmatrix} 0 \\ \rho _{j} \mathbf {F} + \boldsymbol{\nabla}\cdot \hat{\sigma} - \boldsymbol{\nabla} p_{ \mathrm{CR}} \\ \rho _{j} \mathbf {v}_{j}\cdot \mathbf {F} + \boldsymbol{\nabla}\cdot \left (\mathbf {v}_{j} \cdot \hat{\sigma}\right ) - \boldsymbol{\nabla}\cdot \mathbf {Q} - R_{ \mathrm{L}} - \mathbf {v}_{j}\cdot \boldsymbol{\nabla}p_{\mathrm{CR}} \\ 0 \end{bmatrix} + \begin{bmatrix} S_{j}^{\mathrm{c}} \\ \mathbf {S}_{j}^{\mathrm{m}} \\ S_{j}^{\mathrm{e}} \\ \mathbf {A}_{j} \end{bmatrix} \end{gathered}$$

(2)

Here, $j$ gives the particle species, $\mathbf {v}_{j}$ represents its velocity, $\rho _{j}$ its mass density, $p_{j}$ its pressure, and the total energy is

$$\begin{aligned} E_{j} = e_{j} + \frac{1}{2} \rho _{j} v_{j}^{2} + \frac{B^{2}}{2} \end{aligned}$$

(3)

with the inner energy $e_{j}$. Furthermore, $\widetilde{I}$ gives the unit tensor, $\hat{\sigma}$ the stress tensor, $\mathbf {F}$ an external force per unit mass and volume, $\mathbf {Q}$ the heat flow, while $S_{j}^{x}$ gives sources and sinks caused by charge exchange in the continuity ($x = \mathrm{c}$), energy ($x = \mathrm{e}$), and momentum ($x = \mathrm{m}$) equation. ⊗ describes the dyadic/tensorial product. To account for cooling effects $R_{\mathrm{L}}$ represents a cooling function. The subscript $\mathrm{rad}$ accounts for radiation transport of momentum and energy coupling, while $\mathrm{CR}$ accounts for cosmic rays. In addition, $P_{1}$, $P_{2}$, $P_{3}$, and $P_{4}$ are constants, while $\mathbf {A}_{j}$ describes the ambipolar diffusion between neutrals and ions.

If the heat flux is included, the closure for the above set of Euler equations is much more complicated. However, usually the heat flux is assumed to be zero and the ideal gas law is used instead

$$\begin{aligned} p_{j} = (\gamma - 1) e_{j}. \end{aligned}$$

(4)

In the case of ideal HD modeling the following conditions must apply: $\rho \neq 0$, $\mathbf {v}\neq 0$, $E \neq 0$, $P_{1} = P_{2} = P_{3} = P_{4} = 0$, $R_{L} = S_{j}^{x} = \mathbf {F}_{{\mathrm{rad}}}= \mathbf {E}_{{\mathrm{rad}}} = 0$, as well as $\mathbf {B} = 0$. In case of ideal MHD modeling the same criteria apply, however $\mathbf {B} \neq 0$.

Several (M)HD codes exist in the literature (see Kleimann et al. 2022). In the following, the codes utilized to model the astrospheres of planet-hosting cool stars (Sects. 2 and 4), massive runaway stars (Sect. 5), and relativistic objects (Sect. 6) are briefly discussed. In addition, further information on analytic and numerical modeling of the heliosphere can be found in Kleimann et al. (2022).

2.1 cronos and hyperion

cronos is a C++-based semi-discrete finite-volume code (Kissmann et al. 2018) specifically developed to tackle space and astrophysical problems (see, e.g., Röken et al. 2015; Scherer et al. 2015, 2016a,b; Kleimann et al. 2017; Baalmann et al. 2021). The single-fluid MHD computation results presented in this chapter are based on solving the 3D ideal MHD equations given by:

$$\begin{gathered} \frac{\partial}{\partial t} \begin{bmatrix} \rho \\ \rho \mathbf {v} \\ E \\ \mathbf {B} \end{bmatrix} + \boldsymbol{\nabla} \cdot \begin{bmatrix} \rho \mathbf {v} \\ \rho \mathbf {v} \otimes \mathbf {v} + \left (p + \frac{B^{2}}{2}\right ) \widetilde{I} - \mathbf {B} \otimes \mathbf {B} \\ \left (\rho \mathbf {v}\otimes \mathbf {v} + \frac{\gamma}{\gamma -1} p + \frac{B^{2}}{2}\right )\mathbf {v} - \left (\mathbf {B}\cdot \mathbf {v}\right ) \mathbf {B} \\ \mathbf {v}\otimes \mathbf {B} - \mathbf {B}\otimes \mathbf {v} \end{bmatrix} = \begin{bmatrix} 0 \\ \mathbf {0} \\ -R_{\mathrm{L}} \\ 0 \end{bmatrix} \end{gathered}$$

(5)

with an Harten-Lax-van Leer (HLL) Riemann solver and a second-order Runge-Kutta scheme. The cooling and heating functions follow Schure et al. (2009) and Kosiński and Hanasz (2006), respectively, and take into account heat conduction, radiation, and in particular the cooling by inelastic currents and photoionization heating. Modeling based on the latter is performed with the one-fluid module of the code (see also Scherer et al. 2020). The computational grid was arranged in spherical coordinates ($r,\vartheta,\varphi$), producing a spherical polyhedron with equidistant intervals in radius $r$ and equiangular intervals in the polar ($\vartheta$) and azimuthal ($\varphi$) angles. The number of model cells and associated cell sizes were changed according to the model runs performed.

hyperion, on the other hand, is a Fortran-based code that generates synthetic skymaps in different observables from astrosphere models by calculating the respective radiative emission from the model’s cells and summing the corresponding radiances of all model cells that appear within a virtual pixel of the synthetic detector (Baalmann et al. 2020).

Both codes have been utilized to produce results discussed in Sects. 3, 4, and 5.

2.2 3d kinematic mhd model

Based on the pioneering work of Baranov and Malama (1993), the 3d kinematic mhd model, a self-consistent kinetic-gas dynamics model describing the SW/LISM interaction, has been developed by Izmodenov and Alexashov (2015). Here, the partially ionized interstellar plasma is applied as neutral atomic hydrogen and charged plasma consisting of protons, electrons, and helium.

The neutral component is described kinetically with the help of a velocity distribution derived based on the kinetic equation given by

$$ \begin{aligned} &\mathbf {w}_{\mathrm{H}} \cdot \frac{\partial f_{\mathrm{H}}}{\partial \mathbf {r}}+ \frac{\mathbf {F}_{\mathrm{r}}+\mathbf {F}_{\mathrm{g}}}{m_{\mathrm{H}}} \cdot \frac{\partial f_{\mathrm{H}}}{\partial \mathbf {w}_{\mathrm{H}}}=- \nu _{\mathrm{ph}} f_{\mathrm{H}}\left (\mathbf {r}, \mathbf {w}_{\mathrm{H}} \right ) \\ &-f_{\mathrm{H}} \cdot \int \left |\mathbf {w}_{\mathrm{H}}-\mathbf {w}_{ \mathrm{p}}\right | \sigma _{\mathrm{ex}}^{\mathrm{HP}}\left (\left | \mathbf {w}_{\mathrm{H}}-\mathbf {w}_{\mathrm{p}}\right |\right ) f_{ \mathrm{p}}\left (\mathbf {r}, \mathbf {w}_{\mathrm{p}}\right ) \mathrm{d} w_{ \mathrm{p}} \\ &+f_{\mathrm{p}}\left (\mathbf {r}, \mathbf {w}_{\mathrm{H}}\right ) \int \left |\mathbf {w}_{\mathrm{H}}^{*}-\mathbf {w}_{\mathrm{H}}\right | \sigma _{ \mathrm{ex}}^{\mathrm{HP}}\left (\left |\mathbf {w}_{\mathrm{H}}^{*}- \mathbf {w}_{\mathrm{H}}\right |\right ) f_{\mathrm{H}}\left (\mathbf {r}, \mathbf {w}_{\mathrm{H}}^{*}\right ) \mathrm{d} \mathbf {w}_{\mathrm{H}}^{*}, \end{aligned} $$

(6)

with $\mathbf {F}_{\mathrm{r}}$ and $\mathbf {F}_{\mathrm{g}}$ as the forces of solar radiation pressure and gravitation, respectively, $f_{\mathrm{p}}(\mathbf {r}, \mathbf {w}_{\mathrm{H}})$ as the local Maxwellian distribution function of protons, $\sigma _{\mathrm{ex}}^{\mathrm{HP}}(u)$ as the effective charge-exchange cross section with $u$ as the relative atom-proton velocity, and $\nu _{\mathrm{ph}} = 1.67 \times 10^{-7}(R_{\mathrm{E}}/R)^{2}~\text{s}^{-1}$ as the photoionization rate.

The plasma component, on the other hand, is described in the context of an ideal MHD approach (see Eq. (5)) also taking into account the charge exchange with the interstellar hydrogen based on integrals of the H-atom velocity distribution function (see, e.g. Korolkov and Izmodenov 2021).

The 3d kinematic mhd model has been utilized to produce results discussed in Sects. 4, and 5.

2.3 pluto and radmc-3d

pluto is a C-based code particularly suitable for time-dependent, explicit computations of highly supersonic flows in the presence of strong discontinuities that can be employed for different regimes such as classical, relativistic unmagnetized, and magnetized flows (Mignone et al. 2007). The code has been successfully employed in the context of stellar and extra-galactic jets, radiative shocks, accretion disks, and magneto-rotational as well as relativistic Kelvin-Helmholtz instabilities.

Further, to directly compare observations with models, the gas-dynamical outputs are often post-processed with Monte-Carlo radiative transfer tools such as the radmc-3d code (Dullemond 2012), a highly flexible code to compute predictions for observable images and spectra.

Both codes have been utilized to produce the results discussed in Sect. 5.

2.4 3d rmhd pwn

The 3D Reduced MHD (rmhd) model is an incompressible fluid model of plasma behavior. In contrast to a full MHD model the rmhd is simpler, and thus is computationally more efficient than other models (e.g., Oughton et al. 2017). A rmhd pulsar wind nebular (pwn) model (e.g., Olmi and Bucciantini 2019; Barkov et al. 2019; Bucciantini et al. 2020) is the basis of the results discussed in Sect. 6.

3 The Heliosphere as a Test Case for Astrospheric Modeling

Previously cronos, in the context of astrospheric modeling, mainly has been used to model the astrospheres of hot OB-stars (e.g., Scherer et al. 2015, 2020). However, to validate its applicability for the cool-star regime, as a first step the heliosphere has been modeled.

The following model efforts result from a single-fluid 3D model (spherical coordinates) with a radius of $r_{\max} = 1000~\text{AU}$ around the Sun, including values of the number density $n$, the thermal pressure $p_{\mathrm{therm}}$, as well as the vectors of the flow velocity $\mathbf {v}$ and the magnetic flux density $\mathbf {B}$. The values are normalized to $r_{0} = 1~\text{AU}$, $n_{0} = 1~\text{cm}^{-3}$, $v_{0} = v_{\mathrm{A}}(B_{0},n_{0}) = 21.8124~\text{km/s}$, $B_{0} = 1~\text{nT}$, $\rho _{0} = m_{\mathrm{p}}/(1~\text{cm}^{3}) \cdot v_{0} \approx 3.65 \times 10^{-18}~\text{g}/(\text{cm}^{2}\,\text{s})$.

Furthermore, the temperature is given by $T = T_{0} \cdot p_{\mathrm{therm}}/n$, with $T_{0} = v_{0}^{2} \cdot m_{\mathrm{p}}/(2 k_{\mathrm{B}}) \approx 28811~\text{K}$, $k_{\mathrm{B}}$ the Boltzmann constant, and the normalized quantities $p_{\mathrm{therm}}$ and $n$. The total pressure $p_{\mathrm{tot}}$ is composed of the thermal pressure $p_{\mathrm{therm}}$, the dynamic pressure $p_{\mathrm{dyn}}$, and the magnetic pressure $p_{\mathrm{mag}}$. The latter two are given by $p_{\mathrm{mag}} = (B_{0} B)^{2}/(2 \mu _{0} p_{0})$ and $p_{\mathrm{dyn}} =1/(2 p_{0})\cdot m_{\mathrm{p}}\cdot n_{0}\cdot n \cdot (v_{0} v)^{2}$, respectively, where $\mu _{0} = 4 \pi \cdot 10^{-7}~\text{N/A}^{2}$ represents the vacuum permeability. Note that both pressures here are normalized to $p_{0}$, and that the model is set up in a way that the star is in the origin of the coordinates and the stagnation line is the $y$-axis, with positive $y$-values upwind from the star and negative $y$-values in downwind direction.

The left panel of Fig. 2 displays the resulting modeled heliospheric particle densities along the inflow line based on a single-fluid approach. Visible are the TS and the HP, both highlighted as dashed lines. However, a bow shock around the heliosphere is not present in this simulation. The results show that the ISM is not homogeneous but has a weak density gradient towards the HP. Thus, the results underline the findings of the IBEX mission, where the existence of a bow wave rather than a bow shock was proposed (see McComas et al. 2012). The absence of a BS is expected: in the case of MHD modeling, a shock only evolves when the flow velocity $v$ of the plasma exceeds the fast magnetosonic speed, $v_{\mathrm{f}}$. This result, however, strongly depends on the utilized model. According to Scherer and Fichtner (2014), a BS will occur when a multi-fluid simulation is applied.

Fig. 2

Left panel: Modeled heliospheric density distribution shown within the equatorial plane, with the ISM flowing in from the right-hand side. Two distinct discontinuities can be seen, the HP and the TS. Right panel: Corresponding $n$, $T$, $v$, and $B$ profiles

A more detailed analysis of the heliospheric structure is given in the right panels of Fig. 2: the upper panel shows the modeled density ($n$, black) and temperature ($T$, red) profile along the stagnation line, while the lower panel displays the velocity ($v$, black) and magnetic field values ($B$, red) as a function of heliospheric distance. Based on the density and temperature profiles, the modeled locations of the TS and HP can be determined. As can be seen, the simulation predicts the TS and HP to be at $92~\text{AU}$ and $125~\text{AU}$, respectively. Luckily, both Voyager spacecraft (Voyager 1 and 2) launched at the end of 1977 by now have passed the outer boundary of our solar system. Based on observations, Voyager 1 passed the TS at $94~\text{AU}$ and the HP at $122~\text{AU}$, while Voyager 2 passed the TS at $84~\text{AU}$ and the HP at $119~\text{AU}$. Thus, the model results are in good agreement with the observations (see Richardson et al. 2022).

3.1 The “Paleo”-Heliosphere and Its Connection to Nearby Astrospheres

Due to the combination of several observations, it is commonly expected that our solar system is located inside one of the four spiral arms of our Galaxy. Many new stars are being formed inside these spiral arms while other stars end their shell-burning phase with a supernova explosion. Because the position of the heliosphere inside the Milky Way varies, every 70 to 100 million years, it passes through one of the four arms (e.g., Beer et al. 2012) and, thus, passes regions of higher/lower density ISM conditions, leading to changes in the heliospheric shape. Unfortunately, there are no observational data available that would give an insight into the state of the ISM in the past.

Thus, modeling has to rely on sophisticated guesses (Frisch 2006). Borrmann and Fichtner (2005) were the first to study the dynamics of the heliosphere regarding the solar activity on time scales of months up to the Schwabe cycle (11-years). They further investigated heliospheric changes caused by temporal variations of the local interstellar medium. In addition, Scherer et al. (2008) discuss changes in the interstellar proton and hydrogen density as well as the inflow speed of the ISM. HD simulations of an interstellar wave moving over the heliosphere, increasing the three quantities by a factor of four, led to the results shown in Fig. 3. While the upper panel corresponds to the recent heliosphere, the lower panel shows the effect of the perturbed ISM.

Fig. 3

Proton number density of the recent heliosphere (upper panel) in comparison to the “paleo”-heliosphere as discussed in Scherer et al. (2008) (lower panel). Reprinted from Scherer et al. (2008) with permission from Elsevier

For this particular case, the TS is located at 11 AU, while HP and BS can be found at 16 and 22 AU, respectively. Also, the flux of galactic cosmic rays and the hydrogen penetrating the heliosphere are increased by a factor of ten at the vicinity of Earth (see Scherer et al. 2008).

The paleo-heliosphere simulations discussed by Scherer et al. (2008) can be used as a guideline for astrospheric simulations. Nevertheless, a complete 3D MHD simulation including the interstellar neutrals is needed to give reliable results for the cosmic ray and hydrogen flux at the position of a habitable planet. The astrospheres around nearby stars are not directly observable; however, Wood et al. (2005) showed that astrospheric absorption in the stellar Lyman-$\upalpha$ emission line provides a sensitive tool to measure the stellar mass-loss rates of cool stars. In stars cooler than the Sun, the Lyman-$\upalpha$ emission line dominates the stellar UV spectrum (e.g., Linsky 2019); however, Lyman-$\upalpha$ photons traveling through the stellar astrosphere might be absorbed by the hydrogen wall, where inflowing neutral hydrogen piles up. Thus, for an observer looking at the star from the outside, the hydrogen wall absorption is blue-shifted relative to the interstellar medium absorption. The Lyman-$\upalpha$ photons further travel through the (local) interstellar medium where optically very thick hydrogen produces saturated absorption centered at the velocity of the interstellar gas. As the Lyman-$\upalpha$ photons pass through the heliosphere to an observer looking outward to the star, there is red-shifted absorption produced by the heliospheric hydrogen wall, which is slowed down relative to the interstellar medium (see upper panels of Fig. 4). Thus, by fitting the Lyman-$\upalpha$ absorption profile to the astrospheric models, one can gather information about the nearby astrospheres (lower panels of Fig. 4).

Fig. 4

The Lyman-$\upalpha$ emission line. Upper panels (second row): The left panel shows the Lyman-$\upalpha$ line as emitted by the star, the second panel shows that the astrospheric absorption is blue-shifted relative to the line of sight component of the interstellar flow velocity, the third panel represents the travel of Lyman-$\upalpha$ photons through the LISM finally being absorbed by interstellar neutral hydrogen at the projected flow velocity of the LISM. The heliospheric hydrogen wall produces red-shifted absorption, the right panel shows the effect of the three absorption components producing red shift in the emission line of the star. Figure adapted from Linsky and Wood (2014), licenced under the Creative Commons Attribution 4.0 License. The lower panels show a comparison of the observed Lyman-$\upalpha$ line profiles of selected stars (shaded in blue) and the reconstructed stellar line profiles (gray lines). Solid black lines represent the fitted absorption profile. Figure adapted from Wood et al. (2005) ©AAS. Reproduced with permission

3.2 $\text{H}\upalpha$ and the Projection of the Heliosphere at 1 pc

As a prototype of astrospheres around G-type stars, it is of the utmost interest to observe the heliosphere from afar. To generate 2D synthetic sky maps of the heliosphere that can be compared to actual observational data, the cronos-based 3D MHD results need to be line-of-sight integrated (Baalmann et al. 2020).

By rotating the cronos output about all three axes and shifting it to a desired position in the virtual sky, the orientation and location of the modeled star can be reproduced. In the case of the heliosphere, it is interesting to choose a close but non-zero distance, e.g., $d_{0}=1~\text{pc}$, and vary the rotational angles in order to get an outsider’s view from different positions around the heliosphere. At this distance, a spherical grid with a radius of $r_{\max}=900~\text{AU}$ would appear as a disk with a radius of $\phi _{\max}=0.25^{\circ}=15'$. Therefore, a projection of the heliosphere would require a pixel grid with edge lengths $30'\times 30'$.

Before the heliosphere can be projected onto this pixel grid, its optimal resolution must be calculated. Because cronos models are typically realized on spherical grids with constant radial and angular cell dimensions, $(\Delta r, \Delta \vartheta , \Delta \varphi )$, $\Delta \vartheta =\Delta \varphi$, the lateral distance between cells varies within the model. The largest distance between model cells, $(\Delta s)_{\max}$, lies at the outside of the model, $(\Delta s)_{\max} = \left (r_{\max}-\frac{1}{2}\Delta r\right ) \tan (\Delta \varphi )$. An example grid with $\Delta r=3~\text{AU}$ and $\Delta \varphi =3^{\circ}$ yields $(\Delta s)_{\max}\approx 47~\text{AU}$. To avoid empty pixels, the pixel grid resolution shall be coarser than the largest apparent distance between model cells,

$$ (\Delta \phi )_{\max} = 2\arctan \left ( \frac{\frac{1}{2}(\Delta s)_{\max}}{d_{0}-r_{\max}+\frac{1}{2}\Delta r} \right ), $$

(7)

which for the example grid leads to $(\Delta \phi )_{\max}\approx 0.79'$. To improve the grid resolution, it is expedient to linearly interpolate in the two angular directions by a small factor $k\in \mathbb{N}$, typically $k\leq 5$.

The model is projected by first calculating the emission of every model cell and then summing the emissions of all model cells that appear within the respective pixel. Let model cell $i$ appear at latitude $b_{i}$ and longitude $l_{i}$, $\Delta b$ and $\Delta l$ be the latitude and longitude increments of the pixel grid, i.e. the pixel resolution; typically $\Delta l=\Delta b\leq (\Delta \phi )_{\max}$. If $b_{\min}$ and $l_{\min}$ are the minimum latitude and longitude of the pixel grid’s edges, then the model cell appears in the $k_{b}$-th latitude and the $k_{l}$-th longitude pixel, where $k_{b}=\lceil (b_{i}-b_{\min})/(\Delta b)\rceil$ and $k_{l}=\lceil (l_{i}-l_{\min})/(\Delta l)\rceil$.

A typical cronos model of the heliosphere as described in Sect. 2 simulates the number density $n$, temperature $T$ and bulk velocity $\mathbf {v}$ of fully-ionized hydrogen as well as the magnetic field $\mathbf {B}$ in every model cell. With these parameters, some emission processes like bremsstrahlung, cyclotron radiation, or electronic transitions to high-energy states, $i\geq 2$, can be simulated, whereas other emission processes require a multi-fluid approach that includes neutral hydrogen, most notably the Lyman-series electronic transitions.

Some emission processes, for example, bremsstrahlung and cyclotron radiation, emit over a frequency spectrum. Because it is generally not known which emission processes compose an observational image, direct comparisons between synthetic and genuine data are virtually impossible. Thus, emission processes that are restricted to specific frequencies are favored when generating synthetic sky maps. The most commonly used spectral lines are $\text{H}\upalpha$, arising from the electronic transition $3\rightarrow 2$ of hydrogen. While $\text{H}\upalpha$ is not the expected primary emission line of the heliosphere, it can nonetheless serve as a general prototype for other emission lines. It is also an important channel of observations for astrospheres around massive stars (see, e.g., Brown and Bomans 2005, cf. Sect. 5.3), where the ISM is commonly fully ionized. Thus, it is interesting to generate $\text{H}\upalpha$ maps of the heliosphere for comparison.

The radiance of $\text{H}\upalpha$ of a model cell $i$ is calculated by

$$ L_{\Omega ,\text{H}\upalpha,i}=\alpha _{\mathrm{eff},3\rightarrow 2}(T_{i}) \cdot h\nu _{\text{H}\upalpha}n_{i}^{2} \cdot \frac{V_{i}}{4\pi d_{i}^{2}} \cdot \frac{1}{\Delta b \Delta l} \ , $$

(8)

where $n_{i}$ and $T_{i}$ are the number density and temperature of the model cell, $h$ is the Planck constant and $\nu _{\text{H}\upalpha}=456.81~\text{THz}$ the $\text{H}\upalpha$ frequency, $V_{i}$ and $d_{i}$ are the volume and distance of the model cell, $\Delta b \Delta l$ is the pixel’s solid angle, and $\alpha _{\mathrm{eff},3\rightarrow 2}(T_{i})$ is the temperature-dependent effective recombination rate coefficient, which is approximated by $\alpha _{\mathrm{eff},3\rightarrow 2}(T)\approx \alpha _{3}(T)+\sum _{n=4}^{16} \alpha _{n}(T) P_{n,3}$. Here, $\alpha _{i}(T)$ are the temperature-dependent recombination rate coefficients of state $i$ (taken from Mao and Kaastra 2016) and $P_{i,3}$ are the branching ratios from electron level $i$ through level 3, computed as in Baalmann et al. (2020) with data from Wiese and Fuhr (2009), cf. Baalmann et al. (2021).

The selection criterion that verifies if a model cell appears inside a pixel only considers its central coordinate without regard to its geometric shape or size. For example, this leads to errors when most of a model cell’s volume lies outside a line of sight, but its central coordinate appears within the respective pixel. Higher numbers of cells and finer pixel resolutions reduce the impact of this error, which at a count of more than two million cells for the above example grid is negligible.

Projections of the heliosphere at a distance of $d_{0}=1~\text{pc}$ are displayed in Fig. 5. The grid matches the example used in the preceding discussion, simulating the heliosphere with inner and outer radii of $r_{\min}=60~\text{AU}$ and $r_{\max}=900~\text{AU}$. The model grid consists of $280\times 60\times 120$ cells, resulting in radial and angular cell dimensions of $(\Delta r, \Delta \vartheta , \Delta \varphi )=(3~\text{AU}, 3^{\circ}, 3^{\circ})$. The ISM inflow comes from the positive $x$-axis. To increase the available pixel resolution, the model was linearly interpolated in both angular directions by a factor of $k=5$, leading to a division into $162~\text{px}$ per edge, which corresponds to a linear resolution of $\Delta \phi =0.19'$; roughly $3\times 10^{-9}~\text{sr}$ per pixel area.

Fig. 5

$\text{H}\upalpha$ projections of the heliosphere model at a distance of $1~\text{pc}$ face-on ($yx$-plane, left panel) and edge-on ($zx$-plane, centre panel) to the model’s ecliptic plane, and tail-on ($yz$-plane, right panel) to the Sun’s relative motion. The ISM inflow comes from the positive $x$-axis, that is from the top edge (left and middle panels) and into the page (right panel), respectively. $\text{H}\upalpha$ radiance in $\text{erg}\,\text{s}^{-1}\,\text{cm}^{-2}\,\text{sr}^{-1}$, angular extent in degrees. Figure after Baalmann et al. (2020)

The model has been projected at three different sets of angles, corresponding to the model grid’s coordinate planes. In the two left panels, the ISM inflow comes from the top edge, and in the right panel, it moves into the page. The maxima of the $\text{H}\upalpha$ radiance lie in the bow wave in front of the heliopause, where the number density is relatively high, and the temperatures are low. Radiance minima, especially evident in the middle panel, come from the heliotail, where the number density is very low and temperatures are high. The superposition of these two effects can be seen in the right panel, where the low-emission heliotail produces an elongated structure of roughly $1.5\times 10^{-12}~\text{erg}\,\text{s}^{-1}\,\text{cm}^{-2}\,\text{sr}^{-2}$ surrounded by the disk-like maximum of $2.5\times 10^{-12}~\text{erg}\,\text{s}^{-1}\,\text{cm}^{-2}\,\text{sr}^{-2}$ caused by the high-emission bow wave. The small-scale arcs, most apparent in the middle panel along with the $0^{\circ}$-axes, are Moiré patterns and, thus, non-physical structures.

4 Astrospheres of Cool Stars

Cool, low-mass stars are more common within the Galaxy than hotter, more massive ones. Their high number, their long main-sequence lifetime, and their low luminosity, caused by their low masses and small radii, make them perfect targets in the quest for habitable rocky (Earth-like) exoplanets (see, e.g., Gillon et al. 2017; Dittmann et al. 2017). Of particular interest are the smallest and coolest stars, the M-type stars.

Habitable zones (HZs), regions around stars where planetary temperatures are just right to sustain liquid water on the planetary surface, range between $0.03~\text{AU}$ (late-type M-dwarfs) and $0.5~\text{AU}$ (young M-dwarfs). According to West et al. (2004, 2015) and Mohanty et al. (2002) the more prominent the stellar convection envelope and the stellar rotation rates, the stronger the stellar activity. However, for mid-to late-type M-stars, the activity saturates at higher rotational velocities, while above M9, the activity levels decrease significantly (see, e.g., Kay et al. 2016). Furthermore, according to Vidotto et al. (2011), observations of surface magnetic field distributions suggest that young M-dwarfs host weak large-scale magnetic fields dominated by toroidal and non-axisymmetric poloidal configurations (see also Donati et al. 2006). Mid-M-dwarfs, on the other hand, are host to strong, mainly axisymmetric large-scale poloidal fields (Morin et al. 2008). According to Candelaresi et al. (2014), for the majority of these stars, a significantly higher stellar activity compared to the Sun is observed.

Consequently, the exoplanetary radiation environment around active cool stars most likely is much harsher compared to what we know from the Sun (see, e.g., Herbst et al. 2019c). Because of their long lifetimes on the main-sequence, their long stellar activity periods (West et al. 2004), and the small planet-star separations potentially habitable exoplanets in the vicinity of cool stars could be exposed to an enhanced stellar radiation field over millions of years. The latter could significantly affect exoplanetary habitability, for example, due to a hazardous flux of stellar energetic particles (SEPs) influencing its atmospheric evolution, climate and photochemistry (see, e.g., Grenfell 2012; Scheucher et al. 2018, 2020a,b) as well as the altitude-dependent atmospheric radiation dose (see, e.g. Yamashiki et al. 2019; Atri 2020).

Thus, detailed knowledge of the stellar radiation and particle environment and its impact on the (exo)planetary atmospheric chemistry, climate, and induced atmospheric particle radiation field is crucial to assess its habitability and, in particular, potential atmospheric biosignatures. However, up to now, the impact of the astrospheres and the modulation of GCRs within have generally been neglected in such studies.

4.1 Characteristics of Cool Stellar Environments

4.1.1 Properties of Active Cool Stars

The X-ray and UV observations by Chandra, the X-ray Multi-Mirror Mission (XMM-Newton), the Hubble Space Telescope, the Kepler instrument, and most recently the Transiting Exoplanet Survey Satellite (TESS) mission not only gave new insight on the diversity of stars but also on the evolution of Sun-like stars during different evolutionary stages. With this information, we are now on the verge of detecting the magnetic properties of a diversity of cool stars for the first time. Studying stellar atmospheric properties of stars during different stellar phases of evolution will shed light on the history of our star, the Sun.

According to Johnstone et al. (2015c,a,b) and Airapetian and Usmanov (2016) the stellar rotation velocity and the magnetic stars, in particular young Sun-like stars, strongly correlate with age. Thereby, the younger the star, the higher its rotation velocity, magnetic activity, and mass loss rate (see, e.g., Cleeves et al. 2016). Also commonly accepted is the fact that younger solar-type stars have much larger (magnetic) starspots concentrated at higher latitudes than our Sun (see, e.g., Aulanier et al. 2013; Herbst et al. 2021b).

The stellar magnetic activity in the form of X-ray flares is correlated with the stellar rotation velocity (see, e.g., Güdel 2007). Young Sun-like stars with high rotation velocities further show XUV fluxes that are two to three orders of magnitude higher than at the modern Sun (see, e.g., Ribas et al. 2005). However, because of the stellar wind and CMEs during their temporal evolution, the stellar rotation slows down (see, e.g., Weber and Davis 1967).

Stellar rotation-activity-age relations have been studied with the help of UV- and X-ray observations from space (see Zahnle and Walker 1982; Pallavicini et al. 1981, respectively). The following relations have been established so far:

• 1. Rotation period – X-ray luminosity: The higher the rotation period, $P$, the lower the stellar X-ray luminosity, $L_{\mathrm{X}}$. Thereby, $L_{\mathrm{X}} \propto P_{\mathrm{rot}}^{-2.7}$ for a given mass on the main sequence. Since the X-ray flux is mainly generated by the stellar magnetic field, and thus the stellar differential rotation, this relation hints to a close connection between the internal magnetic dynamo and the stellar surface rotation period (see, e.g., Wright et al. 2011, and references therein). However, it is worth mentioning that in the case of stellar rotation rates shorter than a couple of days $L_{\mathrm{X}}$ seems to saturate at a luminosity around $10^{-3} L_{\mathrm{bol}}$ (see Wright et al. 2011). The underlying cause, however, up to now, is poorly understood. A reasonable explanation may, however, be an internal magnetic dynamo threshold. Further, according to Pizzolato et al. (2003), a unified relation in the form of $P_{\mathrm{rot}}^{\mathrm{sat}} \approx 1.2 \left (L_{\mathrm{bol}}/L_{\odot}\right )^{-\frac{1}{2}}$ for the saturation exists for all cool main-sequence stars.

• 2. X-ray luminosity – stellar age: The older the star the lower its X-ray luminosity. Thereby, for Sun-like stars $L_{\mathrm{X}} \propto t^{-1.5}$ (see, e.g., Güdel et al. 1997), which is a direct consequence of the stellar spin-down.

• 3. Stellar age – stellar spin-down: The older the star, the slower its rotation period (Skumanich 1972). Based on cluster samples, for Sun-like stars on average $P_{\mathrm{rot}}\propto t^{0.6}$ (see, e.g., Ayres 1997).

Johnstone et al. (2015c,a), Tu et al. (2015), and most recently Johnstone et al. (2021) re-analyzed the data available by adopting a numerical solar-wind model to stars with different levels of stellar activity and magnetic fluxes. The left panels of Fig. 6 show predicted rotational evolution tracks of cool stars of different stellar masses (i.e. $1~\text{M}_{\odot}$, $0.75~\text{M}_{ \odot}$, $0.5~\text{M}_{\odot}$, and $0.25~\text{M}_{\odot}$; from top to bottom) and the corresponding evolutionary tracks of the $L_{\mathrm{X}}$ values in the right panels. As can be seen, in the case of a Sun-analog ($1~\text{M}_{\odot}$; upper panels) the X-ray emission of the slow rotator tracks (in blue) rapidly decay, while those of a fast rotating Sun-analog decrease much later in its evolution (purple curve). However, the X-ray luminosity values converge at later stellar evolutionary stages (after several hundred Myr) because the stellar rotation rates converge (see left panel). It also shows that the most extensive $L_{\mathrm{X}}$-value spread between slow and fast rotating Sun-analogs occurs between a few tens to a few hundred Myr. However, for the lower mass cases this difference is much smaller because the saturation threshold is located at a much lower rotation rates. In particular for M-stars the slow and fast rotator tracks for $L_{\mathrm{X}}$ are almost identical until 2 Gyr ($0.25~\text{M}_{\odot}$, lower right panel).

Fig. 6

Left panels: Predicted rotational evolution tracks of cool stars with stellar masses of $1~\text{M}_{\odot}$, $0.75~\text{M}_{\odot}$, $0.5~\text{M}_{\odot}$, and $0.25~\text{M}_{\odot}$ (from top to bottom). Fast rotator tracks are highlighted in purple, medium rotator tracks in magenta, and slow rotator tracks in blue. The lower boundaries correspond to the evolution of the envelope and upper boundaries to the core rotation. The rotators are defined as the observed 5th, 50th, and 95th percentiles. Observed rotation rates are shown as open circles. Right panels: Corresponding evolutionary tracks of $L_{\mathrm{X}}$. Shaded areas correspond to one standard deviation. For Sun-analogs (upper panels) measured $L_{\mathrm{X}}$ values are given as crosses (see Tu et al. 2015). Figures adapted from Johnstone et al. (2021), reproduced with permission © ESO

4.1.2 Stellar Winds of Cool Stars

To theoretically describe the stellar XUV emission 3D MHD modeling of stellar coronae and stellar winds is mandatory. Such models, however, up to now, are in their infancy: the first results based on a data-driven model of our Suns twin star $\kappa ^{1}$ Ceti most recently were published by Airapetian et al. (2019). One of the most important parameters for such models is stellar magnetism.

The first study of the average large-scale magnetic field of Sun-like stars utilizing the Zeeman Doppler imaging (ZDI, see, e.g., Donati and Brown 1997) method was presented by Vidotto et al. (2014). Based on a sample of 73 late F-, G-, K- and M-stars, in the pre-main and the main-sequence they demonstrated that the large-scale magnetic field $\left <|B|\right >$ of Sun-like stars decays with age and rotation period (see panels of Fig. 7). Statistically speaking, thereby the average large-scale magnetic field strength scales with

$$\begin{aligned} \left < |B|\right > &\propto t^{-0.655\pm 0.045}\\ &\propto P_{\mathrm{rot}}^{-1.32\pm 0.14}, \end{aligned}$$

(9)

where both provide constraints on the evolution of the average magnetic field strength.

Fig. 7

Left panel: Correlation between the average large-scale magnetic field strength and stellar age, for the non-accreting stars in the sample of 73 late F-, G-, K- and M-stars. Right panel: Correlation between the average large-scale magnetic field strength and the rotation period. Figures 2 and 3 from Vidotto et al. (2014)

A detailed review on the evolution of the solar wind and applications for the stellar wind of Sun-like stars most recently has been published by Vidotto (2021). Further stellar wind modeling of our nearest neighbor Proxima Centauri has been performed by Garraffo et al. (2016) showing that the stellar wind dynamic pressure might be three to four orders of magnitude higher than the solar wind pressure at Earth, which in the end would have a significant impact on the exoplanetary atmosphere of Proxima Centauri b.

4.2 Modeling the Astrospheres of Planet-Hosting Cool Stars

Utilizing cronos (see Sect. 2), Herbst et al. (2020b) modeled the astrospheres of the three M-stars V374 Peg, Proxima Centauri, and LHS 1140 for the first time. In particular, the latter two are hosts of potential Earth-like exoplanets within the habitable zone (HZ) of their host star. It showed that the astrospheres of cool stars could be rather diverse. While V374 Peg drives an astrosphere of several thousands of AU, the astrosphere of LHS 1140 would fit well within the orbit of Neptune. As shown in Fig. 8, Herbst et al. (2020b) further found the astrosphere of our nearest neighbor Proxima Centauri to be well-comparable to the heliosphere.

Fig. 8

Modeled density distribution within the equatorial plane of the heliosphere (upper) and Proxima Centauri (lower), and LHS1140 (on top). Adapted from Herbst et al. (2020b) ©AAS. Reproduced with permission

4.2.1 Role of the Azimuthal Component of the Stellar Wind Magnetic Field on the Global Structure of a Cool-Star Astrosphere

The azimuthal component of the stellar magnetic field can seriously impact the global shape of a stellar astrosphere. Based on the 3d kinematic mhd model Korolkov and Izmodenov (2021) most recently have investigated how the astrospheric shape depends on the stellar magnetic field, whereas Opher et al. (2015) and Drake et al. (2015), for the heliosphere, showed that this effect is essential for the plasma flows in the inner heliosheath.

Thereby the magnetic force $\mathbf {F}_{\mathrm{mag}} = ([\nabla \times \mathbf {B}] \times \mathbf {B})/(4 \pi )$ in the region beyond the TS, considering the problem of the hypersonic spherically symmetric stellar wind interaction with the surrounding interstellar medium at rest, can be estimated. In the subsonic region beyond the TS this leads to $V \sim 1/R^{2}$, $\rho \sim \rho _{\infty}$, and $p \sim p_{\infty}$, where $V$, $\rho$, $p$ are the velocity, density and pressure of the stellar wind, respectively, whereas $\rho _{\infty}$ and $p_{\infty}$ are the interstellar gas density and pressure, respectively. The latter can be used to calculate the frozen-in magnetic field in the kinematic approximation by solving $\nabla \times [\mathbf {V} \times \mathbf {B}]=0$ and assuming that the magnetic field is parallel to the stellar wind velocity vector. The solution above for $R< R_{\mathrm{TS}}$ has been obtained by Parker (1958). For the subsonic region beyond the TS ($R > R_{\mathrm{TS}}$), the following solution can be found: $B_{R} \sim 1/R^{2}$, $B_{\phi}\sim R\sin \theta$, and $B_{\theta}= 0$, where $\theta$ is the polar angle counted from the stellar rotational axis (x-axis) and $\phi$ is the azimuthal angle. As can be seen, the plasma $\beta$ decreases with distance and, therefore, a strong influence of the magnetic field on the plasma should be expected.

A detailed study on the effects of the azimuthal magnetic field component for the considered simplified case has been performed by Golikov et al. (2016). A more general investigation taking into account the motion of the LISM with respect to the star most recently has been published by Korolkov and Izmodenov (2021). In this study, the stellar wind at the inner boundary of the computational domain has been assumed to be spherically-symmetric and supersonic, the azimuthal component of the stellar magnetic field at the inner boundary has been assumed to follow the Parker spiral solution given as

$$ B_{\varphi}= B_{\varphi ,\mathrm{E}} \left (\frac{R_{\mathrm{E}}}{R} \right ) \sin \theta , $$

(10)

where $B_{\varphi ,\mathrm{E}}$ is the strength of the magnetic field at the distance of $R_{\mathrm{E}}$. The IMF and the radial component $B_{R}$ are assumed to be zero. The latter is not critical since $B_{R}\sim 1/R^{2}$ is small compared to $B_{\varphi}$ at large heliocentric distances. Based on Eq. (10), Korolkov and Izmodenov (2021) could show that the solution depends only on three dimensionless parameters: (1) the gas-dynamic Mach number in the ISM, $M_{\infty}$, (2) the Alfvénic Mach number in the stellar wind, $M_{\mathrm{A}}$, and (3) the parameter $\gamma$, the ratio of specific heat fluxes. The results of this study based on $M_{\mathrm{A}} = 12$ (that is close to the solar case) and various values of $M_{\infty}$ are shown in the panels of Fig. 9. The $x$-axis is directed against the incoming flow, while the $z$-axis is the stellar rotation axis. As can be seen in the left panel, for the small values $M_{\infty}$ the stellar wind flow is confined in the slightly distorted tube. The bigger the Mach number the larger the distortion of the tube (middle panel). At a critical value of the Mach number $M_{\text{crit}, \infty}$ the flow pattern changes its structure: for $M_{\infty} < M_{\text{crit}, \infty}$ the pattern has a tube-like structure, while for $M_{\infty}> M_{\text{crit}, \infty}$ the astropause has an open structure in the tail region (right panel).

Fig. 9

The isolines of density and pressure for $M_{\mathrm{A}} = 12$ and different values of $M_{\infty}$ are shown in the plane containing the freestream velocity vector and the star’s axis of rotation. The streamlines of the stellar wind are marked in white, and the incoming stream in black. The contact of the lines can determine the shape of the astropause. The interstellar medium flows from right to left. Figure 3 (1 – a, b, c) from Korolkov and Izmodenov (2021)

Further, the value of $M_{\text{crit}, \infty}$ increases with the decrease of the Alfvénic Mach number that corresponds to the effective increase of the stellar magnetic field. At $M_{\infty} =1$ another bifurcation of the flow pattern appears. The bow shock is formed in the upwind part of the interstellar medium, and the Mach disk is formed in the tail region.

4.2.2 Role of the ISM Field Direction on the Shape of a Cool-Star Astrospheres

One of the open scientific questions is the impact of the interstellar magnetic field direction on the shape of astrospheres and how observations can correctly be interpreted.

Utilizing cronos, in a first step we modeled the astrosphere of Proxima Centauri with the same stellar wind and ISM parameters discussed in Herbst et al. (2020a) (see Sect. 4.2) but with different orientation of the interstellar magnetic field (see also Ratkiewicz et al. 1998). Exemplarily, Fig. 10 shows the model efforts for an ISM field direction of $(\vartheta ,\varphi )=(90^{\circ},45^{\circ})$ (middle panel) and $(\vartheta ,\varphi )=(90^{\circ},135^{\circ})$ (right panel).

Fig. 10

Influence of the LISM field direction on the astrosphere of Proxima Centauri. Left panel: The direction of the magnetic field (shown as colored arrows) of two selected cases $(\vartheta ,\varphi )=(90^{\circ},45^{\circ})$ (in red) and $(90^{\circ},135^{\circ})$ (in blue). Middle and right panel: Model results of the corresponding astrospheres of Proxima Centauri

The direction of the interstellar magnetic field plays a crucial role in the structure of an astrosphere, thus interpreting astrospheric observations. However, to get the complete picture, more detailed studies are mandatory. More advanced investigations are currently in progress.

4.3 Modulation of Galactic Cosmic Rays in Cool Star Astrospheres

Most recently, the impact of cosmic rays and space weather on exoplanetary atmospheres and their effect on atmospheric chemistry and, thus, observable biosignatures and planetary habitability came into the focus of exoplanetary studies (e.g. Herbst et al. 2019b; Airapetian et al. 2020; Atri 2020; Yamashiki et al. 2019). With small astrospheres like those of LHS1140, the question arose whether GCRs might play a much more significant role than previously thought.

In the heliosphere, the turbulent HMF modulates the energy spectrum of GCRs below $\sim 40~\text{GeV}$ (see Engelbrecht et al. 2022). However, GCRs may play an essential role within, for example, Earth-like exoplanetary atmospheres because of their potential impact on atmospheric ionization, chemistry, and thus habitability.

Based on an analytic approach, Sadovski et al. (2018) found that GCRs with energies below $1~\text{TeV}$ are not present at the orbit of Proxima Centauri b. Similar results are derived when the so-called Force-Field approximation (e.g., Caballero-Lopez and Moraal 2004) is applied. Thus, in theory, GCRs should not play a significant role in studying the impact of CRs on exoplanetary atmospheres. However, previous investigations employed a much stronger AMF as it might be the case. Therefore, as a first step, Herbst et al. (2020b) utilized the stellar wind speed and magnetic field distributions along the inflow line provided by the cronos efforts shown in Fig. 8 in order to numerically investigate the modulation of GCRs within the ASPs of V374 Peg, Proxima Centauri, and LHS 1140.

Determining the modulation of GCRs within the heliosphere usually is based on solving the transport equation by Parker (1965). Since for astrophysical environments little is known about the diffusion coefficients and magnetic field, in Herbst et al. (2020b) a 1D (drift-less) version of the transport equation was solved numerically employing stochastic differential equations (SDEs, see Strauss and Effenberger 2017, and references therein).

Therefore, a distribution undergoing only diffusion and convection is solved. The transport equation can then be reduced to the Fokker-Planck equation (see, e.g., Caballero-Lopez and Moraal 2004):

$$ \frac{\partial f}{\partial t} = - v_{\mathrm{{sw}}} \frac{\partial f}{\partial r} + \frac{1}{r^{2}} \frac{\partial}{\partial r} \left (r^{2} \kappa _{rr} \frac{\partial f}{\partial r}\right ) + \frac{1}{3 r^{2}} \frac{\partial}{\partial r}\left ( r^{2} v_{\mathrm{{sw}}}\right ) \frac{\partial f}{\partial \ln p}, $$

(11)

where $f$ is the GCR distribution function, $V$ is the convection velocity, $\kappa _{rr}$ is the radial diffusion coefficient and $p$ is the momentum. The Ito stochastic differential equations can now be written as (e.g. Engelbrecht and Di Felice 2020)

$$ \Delta r = \left (\frac{\partial \kappa _{rr}}{\partial r} - V\right ) \Delta t + \sqrt{2\kappa _{rr}}\Delta W, $$

(12)

for the time-backward case, and for the change in momentum

$$ \Delta p = \frac{p}{3}\frac{\partial V}{\partial r} \Delta t. $$

(13)

The Wiener process is given by the standard normal distribution, $\mathrm{d}w_{r}$, with a mean of zero and a standard deviation of one (e.g. Zhang 1999)

$$ \Delta W = \sqrt{\Delta t}\lbrace \mathrm{d}w_{r},0\rbrace . $$

(14)

The time constraint is chosen so that $V\Delta t_{i} \ll \sqrt{\kappa _{rr} \Delta t_{i}}$, where $t_{i} = t_{i-1} + \Delta t_{i}$ (Strauss and Effenberger 2017). Further, the radial diffusion coefficient $\kappa _{rr}$ can be written in terms of a perpendicular diffusion mean free path (MFP) $\lambda _{\perp}$

$$ K_{rr}=\frac{v \lambda _{\perp}}{3}=\frac{v}{3}\frac{C}{B}, $$

(15)

where as a first approach a Bohm-like diffusion coefficient can be assumed (see, e.g., Zank et al. 2004; Hussein and Shalchi 2014), with $C$ a constant and $B$ the magnetic field of the star, similar to what has been assumed in previous, heliospheric studies (for further discussion see, e.g., Ferreira and Potgieter 2004; Manuel et al. 2014; Ferreira and Manuel 2014; Manuel et al. 2015; Caballero-Lopez et al. 2019).

The magnetic field used in this cosmic ray transport model is calculated by including the stellar magnetic field as a boundary condition in the fluid (MHD) model. An unmodulated boundary spectrum, the so-called local interstellar spectrum (LIS), is required to solve the above equations. For example, in Herbst et al. (2019a) it is assumed that this spectrum is described by the expression employed by Strauss et al. (2011). However, it should be noted that such a spectrum, constructed with near-heliospheric conditions in mind, may not be an accurate description of the LIS in the vicinity of all astrospheres, given the highly likely spatial and energy variations in the galactic distribution of GCRs (see, e.g., Amato and Blasi 2018, and references therein).

Utilizing the method discussed above, Herbst et al. (2020b) investigated the GCR modulation within the three M-star astrospheres. As an example, the GCR flux derived at LHS 1140 b (blue band) computed using the local interstellar spectrum by Strauss et al. (2011) compared to the analytical solution based on the Force-Field approximation (blue line) is shown in Fig. 11, alongside the GCR flux computed for heliospheric conditions at Earth between solar minimum and maximum conditions (purple band). As can be seen, substantial deviations between the utilized methods occur, which from heliospheric studies is not entirely unexpected (see, e.g., Engelbrecht and Di Felice 2020).

Fig. 11

GCR proton spectra for LHS 1140 b computed using the 1D SDE approach (blue band, according to Herbst et al. 2020b) and the force-field approximation (blue line). The computed solar cycle variation of the GCR flux at Earth is shown as a purple band

Other approaches to derive the GCR flux in Sun-like and cool star astrospheres based on theoretic modeling of the stellar wind and a 1D transport code are currently under investigation (see Rodgers-Lee et al. 2021; Mesquita et al. 2021, respectively).

4.4 Turbulence and Cosmic Ray Transport Coefficients in Astrospheres

In order to more fully understand the modulation of GCRs in astrospheres, it is necessary to model, in as realistic a manner as possible, the various processes governing their transport. From studies of the heliospheric modulation of these particles, it has become clear that diffusion both parallel and perpendicular to the heliospheric magnetic field, as well as drift, plays a key role in their transport (see, e.g., Engelbrecht and Burger 2015b; Qin and Shen 2017; Ghanbari et al. 2019), processes that require careful, 3D modeling. Many theories have been proposed to describe the parallel and perpendicular diffusion of charged particles in turbulent plasmas (see, e.g., Shalchi 2009; Ruffolo et al. 2012; Qin and Zhang 2014; Shalchi 2020), and much progress has been made in recent years towards incorporating these theories, along with advancements in our understanding of heliospheric turbulence, into GCR transport models (Oughton and Engelbrecht 2021, and references therein). In many cases, the results of these theories are complicated and mathematically intractable, often taking into account observed details in heliospheric turbulence that may be difficult to ascertain in an astrospherical context (see, e.g., Shalchi et al. 2010; Engelbrecht and Burger 2015a; Dempers and Engelbrecht 2020). Furthermore, in the heliosphere, drift effects due to gradients and curvatures in the heliospheric magnetic field, as well as along the heliospheric current sheet, have long been known to significantly affect the degree to which GCR spectra are modulated (see, e.g., Jokipii and Levy 1977; Jokipii and Thomas 1981; Burger et al. 1985). Although gradient and curvature drift velocities can in principle be calculated from astrospheric magnetic field profiles computed from MHD modeling, it is not yet observationally clear whether astrospherical equivalents of the heliospheric current sheet exist. As an additional complication, drift effects are also known from theory and simulation to be strongly influenced by magnetic turbulence (e.g., Bieber and Matthaeus 1997; Minnie et al. 2007b; Burger and Visser 2010; Tautz and Shalchi 2012; Engelbrecht et al. 2017).

In modeling GCR drift and diffusion coefficients, a core difficulty would be that we cannot measure the various required turbulence or related parameters in the different astrospheres we aim to consider. A potential solution to this would be to calculate them from first principles using our knowledge of the behavior of observable parameters, like the stellar magnetic field, and derived stellar wind parameters, relative to that observed in the heliosphere, for which extensive observations exist (see, e.g., Bruno and Carbone 2016). The first approach to this problem would be to scale turbulence quantities in the astrospheres relative to their behavior in the heliosphere. One way to do this would be to scale astrospheric magnetic variances up or down relative to their heliospheric equivalents in the same manner that the corresponding astrospheric magnetic fields scale up or down relative to the heliospheric magnetic field. This is motivated by the heliospheric observation that these quantities are related (see, e.g., Matthaeus et al. 1986; Zhao et al. 2018; Engelbrecht and Wolmarans 2020). Such an approach would allow us to estimate the basic turbulence quantities required as inputs for the parallel and perpendicular diffusion coefficients mentioned above, as well as the turbulent drift reduction factor. The spatial dependences for such turbulence quantities can be modeled after those yielded by advanced turbulence transport models in the heliosphere (e.g. Usmanov et al. 2016; Wiengarten et al. 2016; Zank et al. 2017; Adhikari et al. 2017), as well as heliospheric observations of these quantities (e.g. Zank et al. 1996; Smith et al. 2001; Pine et al. 2020a,b).

As a demonstration of a first approach to modeling GCR transport coefficients in the manner described above, we consider transport coefficients calculated for the astrosphere of LHS 1140, as treated by Herbst et al. (2020b). Here we consider only the astrospheric magnetic field as a function of radial distance in the ecliptic plane towards the nose of the astrosphere. The square of this quantity is shown in the top left panel of Fig. 12 (solid red line), compared with the square of the magnitude of the heliospheric Parker (1958) field (solid blue line). Note that the astrosphere of LHS 1140 is considerably smaller than the heliosphere, with a TS at $\sim 8.1~\text{AU}$, and an astrosheath extending out to $\sim 28~\text{AU}$. In order to estimate magnetic variances and correlation scales, we scale values observed in the heliosphere at $1~\text{AU}$ by the ratio of the square of the astrospheric magnetic field magnitude and heliospheric magnetic field magnitude at that radial distance, which is here equal to $3.44\times 10^{-7}$. Therefore, for slab and 2D magnetic correlation scales, we apply this scaling to the observations for these quantities reported by Weygand et al. (2011). In contrast, for the total magnetic variance, we apply this scaling to a total magnetic variance of $12~\text{nT}^{2}$, following observations reported by Smith et al. (2006). As a very first approach, slab and 2D variances are calculated under the assumption of a 20:80 variance anisotropy (Bieber et al. 1994). It is doubtful whether this assumption is very accurate, as this ratio is known to vary considerably in the heliosphere (see, e.g., Oughton et al. 2015). A smaller value is probably more realistic, given that the considerably smaller astrospheric magnetic field would be expected to have a smaller influence on spectral transfer than the considerably stronger Parker field (see Shebalin et al. 1983; Oughton et al. 1994). These effects, however, are challenging to ascertain. The total transverse magnetic variance and correlation scales are here spatially scaled using simple power laws, such that $\delta B^{2} \sim r^{-2.4}$, and it is assumed that both MFPs $\lambda _{\mathrm{2D}}$ and $\lambda _{\mathrm{sl}} \sim r^{0.5}$, following the trend of Voyager observations of these quantities in the heliosphere reported by Zank et al. (1996) and Smith et al. (2001). These scalings are demonstrated, for both the heliosphere and LHS 1140, in the top panels of Fig. 12. Note the considerably smaller magnetic variances and correlation scales for LHS 1140 relative to the heliospheric values. Such a scaling neglects the influence of turbulence generated by the formation of pickup ions (see, e.g., Williams and Zank 1994; Isenberg 2005; Cannon et al. 2014), which may influence the transport of low-energy cosmic rays (Engelbrecht 2017), and as such are similar to radial dependences yielded by various turbulence transport models when pickup ion effects are not considered (e.g. Breech et al. 2008; Engelbrecht and Burger 2015b; Adhikari et al. 2017). As a first approach, these radial scalings are not altered beyond the astrospheric termination shock. Although this is motivated by the strong decrease in transverse turbulence observed in the heliosheath (see, e.g., Burlaga et al. 2014, 2018), the present study does not take into account compressive turbulence observed in the heliosheath (see, e.g., Burlaga and Ness 2009), which can arise in the vicinity of the termination shock (see Gutynska et al. 2010; Fichtner et al. 2020) and beyond (e.g. Zhao et al. 2020). This omission is motivated by current limitations in our understanding of the transport and evolution of compressive turbulence in the inner heliosheath and our understanding of the influence of such turbulence on the transport coefficients of charged energetic particles.

Fig. 12

Top left panel: Square of the magnitude of the Parker (1958) heliospheric magnetic field, compared with that of the astrospheric magnetic field of LHS 1140 computed using the cronos MHD code, as employed by Herbst et al. (2020b). Also shown are scalings for the total magnetic variance employed here for the heliosphere and LHS 1140. Top right panel: scalings employed here to model slab and 2D correlation scales in the heliosphere and for LHS 1140. Bottom panels: Diffusion mean free paths and drift length scales calculated for heliospheric parameters (left panel) as well as those assumed for LHS 1140 (right panel). See text for details

The very limited information available to us as to the nature of astrospheric turbulence also places limitations on how GCR transport coefficients can be modeled, as the various scattering theories employed to do so require, as basic inputs, various details as to the underlying turbulence, and the turbulence power spectrum. Detailed observations of heliospheric turbulence can, in principle, be taken into account when deriving diffusion coefficients and often have a large influence on these quantities. Examples of this include turbulence anisotropy (e.g. Bieber et al. 1994), the low-wavenumber behavior of the turbulence power spectra (see, e.g., Shalchi et al. 2010; Qin and Shalchi 2012; Engelbrecht and Burger 2015a; Engelbrecht 2019), potential non-axisymmetry of turbulent fluctuations (e.g. Ruffolo et al. 2008; Strauss et al. 2016), the high-wavenumber behaviour of the turbulence power spectrum (e.g. Engelbrecht and Burger 2010, 2013), and potential effects of dynamical turbulence (e.g. Teufel and Schlickeiser 2003; Shalchi et al. 2004b; Gammon and Shalchi 2017; Dempers and Engelbrecht 2020). Lacking such detailed information, we employ a relatively simple, tractable set of expressions for GCR proton diffusion and drift length scales that have been shown to lead to GCR intensities in good agreement with spacecraft observations in the heliospheric context (see Moloto et al. 2018; Moloto and Engelbrecht 2020; Engelbrecht and Moloto 2021). The MFP parallel to the astrospheric magnetic field is modeled using the quasilinear theory expression derived by Teufel and Schlickeiser (2003) for an assumed slab turbulence power spectral form consisting of a wavenumber-independent energy-containing range, and an inertial range (see also Bieber et al. 1994), given by (Burger et al. 2008)

$$ \lambda _{\parallel} = \frac{3s}{(s-1)} \frac{R^{2}}{k_{\mathrm{m}}} \frac{B^{2}_{0}}{\delta B^{2}_{\mathrm{sl}}} \left [ \frac{1}{4\pi} + \frac{2R^{-s}}{\pi (2-s)(4-s)} \right ], $$

(16)

where $R=R_{\mathrm{L}} k_{\mathrm{m}}$, with $R_{\mathrm{L}}$ the maximal particle gyroradius and $k_{\mathrm{m}} = 1/\lambda _{\mathrm{sl}}$ the wavenumber at which the inertial range, with Kolmogorov spectral index $-s=-5/3$, commences. The quantity $B_{0}$ denotes the magnitude of the background astrospheric magnetic field, while $\delta B^{2}_{\mathrm{sl}}$ denotes the slab variance. The perpendicular MFP is modeled using the nonlinear guiding center (NLGC) theory (Matthaeus et al. 2003) expression derived by Shalchi et al. (2004a)

$$ \lambda _{\perp}= \left [ \alpha ^{2} \sqrt{3\pi} \frac{2\nu -1}{\nu} \frac{\Gamma (\nu )}{\Gamma (\nu -1/2)} \lambda _{\mathrm{2D}} \frac{\delta B^{2}_{\mathrm{2D}}}{B^{2}_{0}} \right ]^{2/3} \lambda ^{1/3}_{\parallel}, $$

(17)

where $\alpha ^{2}=1/3$ (Matthaeus et al. 2003), $\nu =5/6$, and the subscript ‘2D’ denotes a turbulence quantity pertaining to the 2D turbulence power spectrum. Note that Eq. (17) is derived assuming a 2D turbulence spectral form similar to that assumed in the derivation of Eq. (16). Although the assumption of such a form does not entirely reflect what is observed in the heliosphere (see, e.g., Bieber et al. 1993; Goldstein and Roberts 1999), it does lead to simpler analytical forms for $\lambda _{\perp}$. The above expressions also have the advantage of yielding MFPs in reasonably good agreement with numerical test particle simulations for low levels of turbulence (Minnie et al. 2007a). To model the influence of turbulence on particle drift lengthscales, we employ an expression derived by Engelbrecht et al. (2017) from the theory of Bieber and Matthaeus (1997), given by

$$ \lambda _{\mathrm{A}} = \left [1 + \frac{\lambda ^{2}_{\perp}}{R_{\mathrm{L}}^{2}} \frac{\delta B^{2}_{\mathrm{T}}}{B^{2}_{0}}\right ]^{-1}R_{ \mathrm{L}}, $$

(18)

with $\delta B^{2}_{\mathrm{T}}=\delta B^{2}_{\mathrm{sl}}+\delta B^{2}_{\mathrm{2D}}$ denoting the total transverse magnetic variance. This expression also yields results in good agreement with the numerical test particle simulations of Minnie et al. (2007b) and Tautz and Shalchi (2012). For very low levels of turbulence, Eq. (18) reduces to the weak-scattering result, where $\lambda _{\mathrm{A}}=R_{\mathrm{L}}$ (see Forman et al. 1974). For the purposes of comparison, we also consider the expression for the perpendicular MFP employed by Herbst et al. (2020b). In that study, it is assumed that $\lambda _{\perp}$ scales as the inverse of the astrospheric magnetic field magnitude, normalised to a value $\lambda _{0} = a {B_{\mathrm{E}} (1~\text{AU})}/{B_{1}}$, with $B_{\mathrm{E}}$ and $B_{1}$ the magnitudes of the heliospheric and astrospheric magnetic fields at $1~\text{AU}$, and $a$ a free parameter, chosen here to be equal to $0.01~\text{AU}$, similar to previous approaches in heliospheric as well as prior astrospheric modelling (e.g. Ferreira and Potgieter 2003; Scherer et al. 2015). This expression is given by

$$ \lambda _{\perp} = \lambda _{0} \left ( \frac{1\,\mathrm{nT}}{B_{0}} \right ) \left (\frac{P}{P_{0}}\right )^{1/3} $$

(19)

where $P$ the particle rigidity, in units of GV, $P_{0} =1~\text{GV}$, and $B_{0}$ in units of nT.

The resulting GCR proton diffusion MFPs and drift and drift lengthscales are shown in the bottom panels of Fig. 12 at a rigidity of $1~\text{GV}$, as a function of radial distance, with the left panel displaying results calculated for heliospheric parameters, and the right panel for parameters estimated and calculated for LHS 1140. For comparison, quantities calculated for heliospheric parameters are only shown out to $28~\text{AU}$, the extent of the astrosphere of LHS 1140. The parallel MFP (red line) increases steadily with radial distance in the heliosphere, reflecting the steady decrease in the slab variance. The NLGC perpendicular MFP (Eq. (17), solid blue line) displays only a minor increase with radial distance, as the decrease in 2D variance with radial distance is matched by this quantity’s $\lambda _{\parallel}$ dependence. This behavior is markedly different from that displayed by the perpendicular MFP as modeled using Eq. (19) (dashed blue line), which, as expected, displays the same radial dependence as the maximal Larmor radius (solid green line). This, however, leads to this quantity already assuming values more than an order of magnitude larger than those yielded by the NLGC expression. The turbulence-reduced drift coefficient (dashed green lines) does not significantly differ from the weak scattering value ($R_{\mathrm{L}}$) due to the relatively small perpendicular MFP yielded by Eq. (17) at this rigidity within $\sim 10~\text{AU}$, as well as the low turbulence levels beyond this radial distance. Turning to LHS 1140, extremely low turbulence levels lead to a parallel MFP considerably larger than in the heliosphere, showing no evidence of the increase in $B$ associated with the astrospheric termination shock. The radial dependencies of both perpendicular MFP expressions considered here follow that of the inverse of the magnetic field magnitude, or its square, respectively, with that of the NLGC expression tempered by the low values of the turbulence inputs. Equation (19) yields values roughly an order of magnitude larger than Eq. (17) within the termination shock and several orders of magnitude larger than the NLGC perpendicular MFP in the astrosheath. This latter quantity remains small for all radial distances shown, despite its dependence on the large parallel MFP, due to the low turbulence levels. The latter also influence the turbulence-reduced drift scale so that it yields values essentially identical to the weak scattering value, which, intriguingly, is consistently larger than both perpendicular MFPs considered here by several orders of magnitude.

Given that the behavior of turbulence quantities in the astrosphere of LHS 1140 is modeled relying on estimates based on the behavior of these quantities in the heliosphere, the results discussed above should be viewed as highly tentative. Nevertheless, they raise some interesting questions. The differences in the NLGC perpendicular MFP, when compared to that employed by Herbst et al. (2020b), would lead to apparent differences in modulated GCR differential intensity spectra, even when these are modeled using a 1D approach. A comparison of the differences in transport coefficients for LHS 1140, relative to those modeled for the heliosphere, would imply that modulation in astrospheres, particularly LHS 1140, would be very different from what is expected in the heliosphere. For example, we would expect a consistent drift dominated modulation in LHS 1140 because it has such a large drift scale relative to the perpendicular MFP, as opposed to the alternating drift versus diffusion dominated behavior seen going from solar minimum to solar maximum in the heliosphere (see, e.g., Moloto and Engelbrecht 2020). Furthermore, the large parallel MFP for LHS 1140 might mean that parallel diffusion plays a more significant role in the modulation of GCRs than in the heliosphere. As such, this would further imply that, in order to take these effects into account, astrospheric GCR modulation models would need to be solved in 3D, as is done in the heliospheric context.

5 Astrospheres of Massive Stars – The Bow-Shocks of Runaway

Massive stars ($\ge 8~\text{M}_{\odot}$) represent a few percent of all stellar objects. However, their rare occurrence in star-forming processes is inversely proportional to their importance in the cycle of matter in the ISM of Galaxies (Langer 2012). The many-body gravitational interactions in stellar clusters in which they form can eject a significative fraction of all massive stars through the ISM (Blaauw 1961; Gies 1987; Hoogerwerf et al. 2001), and, by wind-ISM interaction, runaway massive stars produce bow shocks which can be observed. Several features distinguish bow shocks around massive stars from their lower-mass counterparts, such as their vast (tens of) parsec-scale size, induced by the strong stellar winds, making any stellar surface magnetic field dynamically unimportant in the nebulae’s shaping and their time-dependence, reflecting the evolution of the stellar surface properties.

5.1 Distorted Circumstellar Wind Bubbles Around Runaway OB Stars

The surroundings of static massive stars produce bubble nebulae (Weaver et al. 1977), which are the finger print of the star’s past evolution onto their local medium (Dwarkadas 2007; Freyer et al. 2006; van Marle et al. 2015). The massive star community, therefore, defines bow shocks as circumstellar bubbles distorted by stellar motion. Because of the peculiar wind properties of massive stars, those bow shocks are qualitatively similar but quantitatively different from the heliosphere. Their careful characterization evolved independently from that of the heliophysics, driven mainly by observations (see, e.g., the optical O[III] arc around the massive runaway OB star $\zeta$ Ophiuchi; Gull and Sofia 1979). Progress in their comprehension are both made of serendipitous observational discoveries and numerical modeling. They constitute a multiwavelength laboratory for observers and numerics, which study permits to constrain both stellar evolution and the local properties of the ISM.

Vocables for stellar wind bubbles persisted in the description of bow shocks around massive stars, made of an inner reverse shock (the termination shock), an outer forward shock (the bow shock), encompassing the wind-ISM interface, which is often referred to as the contact discontinuity but is in the MHD notation a tangential discontinuity. Because both types of discontinuities exist in the MHD framework, the correct phrase for the LISM-wind interface is tangential discontinuity.

Based on Eq. (1), Parker (1963) deduced for a supersonic inflowing interstellar medium with no magnetic field for the termination shock distance $R$:

$$\begin{aligned} R = r_{0} \sqrt{ \frac{\rho _{\mathrm{SW}} v_{\mathrm{SW}}^{2}}{\rho _{\mathrm{ISM}} v_{\mathrm{ISM}}^{2}}}, \end{aligned}$$

(20)

where $\rho$ and $v$ are the density and speed, respectively, the index $\mathrm{SW}$ denotes the solar (stellar) wind, and $\mathrm{ISM}$ the respective values in the ISM, based on the Bernoulli law, where the total pressure along a streamline is constant. In the HD case, there is a streamline passing through the star and parallel to the inflow, on which the closest distance of the termination shock, the stagnation point (or closest distance of the astropause), and the closest distance of the bow shock in a stellar rest frame are located. Equation (20) can easily be reformulated by replacing the stellar wind density by the stellar mass loss rate $\dot{M}$ leading to the formula deduced by (Baranov and Malama 1993) or Wilkin (1996):

$$ R = \sqrt{\frac{ \dot{M} v_{\mathrm{w}} }{ 4 \pi n_{\mathrm{ISM}} v_{\star}^{2} } }, $$

(21)

The above equation determines the TS distance, while the stagnation distance can only be determined by further assuming that potential theory can describe the flow. The bow shock distance can – so far – not be determined analytically.

In the contemporary MHD theories and multifluid models, the TS, AP, and BS distances cannot be determined in general. Note that in ideal 3D MHD, there must not exist a streamline through the star, the closest TS distance, stagnation point, and if it exists the closest BS distance (e.g., Nickeler et al. 2014; Scherer et al. 2020). The reason is that the ISM magnetic field destroys the axis-symmetry of the pure HD flow, as shown, for example, in Fig. 10.

Bow shocks of massive stars permit an independent estimate of the mass-loss rate of massive stars (Kobulnicky et al. 2010 and Gvaramadze et al. 2012a). A growing number of BS around massive stars have been observed since then, either serendipitously (Kaper et al. 1997) or by inspecting the surroundings of stellar clusters, from which many runaway stars are ejected (Gvaramadze and Bomans 2008; Gvaramadze et al. 2011; Gvaramadze and Gualandris 2011; Gvaramadze et al. 2012b; Peri et al. 2015). A multi-wavelength, high-resolution picture of bow shock nebulae (BSN) around OB stars then arose, from X-ray (López-Santiago et al. 2012), ultraviolet (Le Bertre et al. 2012) and radio measures (Benaglia et al. 2010). Nevertheless, most observations were taken in the near infrared waveband, which led to several catalogues (van Buren and McCray 1988; van Buren et al. 1995; Noriega-Crespo et al. 1997b; Peri et al. 2012, 2015; Kobulnicky et al. 2016, 2017, 2018).

Bow shocks of massive stars motivated analytic (Dgani et al. 1996b,a) and numerical gasdynamical simulations which explored their internal structure and non-linear instabilities (Blondin and Koerwer 1998), mostly with multi-purpose Eulerian codes such as pluto (Mignone et al. 2007, 2012; Vaidya et al. 2018), although Lagrangian models exist (Mohamed et al. 2012). Furthermore, full 3D MHD simulations of runaway stars have been performed recently. The results on $\lambda$-Cephei have been discussed in Scherer et al. (2020) and Baalmann et al. (2020) and those on -Cas by Katushkina et al. (2018). Most recently IRC-10414 and Betelgeuse have been modeled by Meyer et al. (2021a).

Notably, the role of electronic heat conduction at the TS involving large temperature gradients ($>10^{6}~\text{K}$) is important (Comerón and Kaper 1998; Meyer et al. 2016, 2017). The penetration of ISM dust in the infrared emission motivated detailed semi-analytic investigations of dust grain physics in the astrospheres of runaway stars (Henney and Arthur 2019a,b,c; Henney et al. 2019). To directly compare observations with models, gasdynamical outputs are often post-processed with Monte-Carlo radiative transfer tools such as the radmc-3d code (Dullemond 2012), see Meyer et al. (2017). Lastly, stellar wind bow shocks have been long suspected to be the site of non-thermal processes such as particle acceleration (Benaglia et al. 2010; del Valle et al. 2015; del Valle and Pohl 2018; Benaglia et al. 2021), although contradictory results exist (De Becker et al. 2017; Binder et al. 2019).

5.2 Astrospheres of Fast-Moving Evolved Massive Stars

The stellar wind properties of high-mass stars are the consequence of complex nuclear reactions at work in their core, engendering much more violent and by far quicker ($\sim \text{Myr}$) changes than that of low-mass stars ($\sim\text{Gyr}$) (Maeder and Meynet 2012). Mainly according to their mass (Ekström et al. 2012), massive stars evolve to a short, cold M-type red supergiant phase, which leads to changes in wind properties modifying their nebula morphology. The new-born red supergiant stellar wind is blown into the freely-expanding main-sequence wind, forming a circular shell by wind-wind interaction, eventually colliding with the main-sequence bow shock (Brighenti and D’Ercole 1995b,a). If the star moves sufficiently fast, it interacts directly with the ambient medium, after a short-lived Napoleon’s hat (Wang et al. 1993) and a dusty (van Marle et al. 2011) red supergiant bow shock is formed (Gvaramadze et al. 2014). The latter can be modeled by time-dependently updating the stellar wind boundary conditions (Meyer et al. 2014b).

While all (runaway) massive stars evolve, only three runaway red supergiant stars with a BS have been reported in the literature so far: IRC-10414, $\upalpha$ Ori (Betelgeuse), and μ Cep, raising the question of the apparent missing bow shock problem. The BS of IRC-10414 (Gvaramadze et al. 2014) is smooth, as a consequence of an external ionising photon field (Meyer et al. 2014a), see Fig. 13. Similarly, the stable appearance of bow shocks of Betelgeuse (Noriega-Crespo et al. 1997a) results in both stabilisation by an external ionization stellar field (Mackey et al. 2014) and of the damping of HD instabilities by the native ISM magnetic field of the Orion arm (van Marle et al. 2014). Its roundness indicates its young age, i.e., Betelgeuse just underwent a stellar phase transition (Mohamed et al. 2012). A mysterious bar, almost parallel to the direction of stellar motion, probably of interstellar origin (Decin et al. 2012)³ completes this picture.

Fig. 13

Optical $\text{H}\upalpha$ emission maps of 3D MHD simulations of the bow shock of the red supergiant star IRC-10414 (top) and Betelgeuse (bottom), using the pluto and radmc-3d codes. Figures 14 and 15 from Meyer et al. (2021a)

Some massive stars $\ge 35~\text{M}_{\odot}$ can further evolve to the violent, hot-winded Wolf-Rayet phase. These unsteady stellar wind nebulae exhibit spectra with line emission (van Marle et al. 2005, 2007; Meyer et al. 2020a; Reyes-Iturbide et al. 2019). Eventually, massive stars can die as a supernova explosion, leaving behind a supernova remnant carrying the imprint of the progenitor’s BS (Meyer et al. 2015, 2020b, 2021b).

5.3 Projections of Model Astrospheres Around Massive Stars

In analogy to Sect. 3.2 the model astrospheres around massive stars must be projected to generate synthetic images that can be compared to actual data. Like with heliospheric simulations, the $\text{H}\upalpha$ radiance of a model cell is calculated by Eq. (8). A similar approach yields the radiances of other emission lines of ionized hydrogen. The total radiance of free-free emission (Rybicki and Lightman 1979) for model cell $i$ follows

$$\begin{aligned} L_{\Omega ,\mathrm{ff,}i}&=\left ( \frac{2\pi k_{\mathrm{B}} T_{i}}{3m_{\mathrm{e}}}\right )^{1/2} \frac{2^{5}\pi e^{6}}{3h m_{\mathrm{e}}c^{3}} n_{i,\mathrm{e}} n_{i, \mathrm{H}^{+}} g_{\mathrm{ff}}(T_{i}) \frac{V_{i}}{4\pi d_{i}^{2}} \frac{1}{(\Delta \phi )^{2}} \end{aligned}$$

(22)

$$\begin{aligned} &\approx 1.435\times 10^{-27}~\frac{\text{erg}\,\text{cm}^{3}}{\text{s}\,\sqrt{\text{K}}} \cdot \sqrt{T_{i}}n_{i}^{2} g_{\mathrm{ff}}(T_{i}) \frac{V_{i}}{4\pi d_{i}^{2}} \frac{1}{(\Delta \phi )^{2}} \ , \end{aligned}$$

(23)

where $n_{i}=n_{i,\mathrm{e}}=n_{i,\mathrm{H}^{+}}$, $T_{i}$, $V_{i}$, $d_{i}$ are the number density, temperature, volume, and distance to the observer of model cell $i$, respectively; $\Delta \phi$ is the linear pixel resolution; $k_{\mathrm{B}}$, $m_{\mathrm{e}}$, $e$, $h$, $c$ are the Boltzmann constant, electron mass, elementary charge, Planck constant and vacuum speed of light, respectively; and $g_{\mathrm{ff}}(T)$ is the temperature-dependent Gaunt factor for free-free transitions, which can be numerically calculated by, e.g., van Hoof et al. (2014). An often-made assumption for the Gaunt factor is $\bar{g}_{\mathrm{ff}}=1.2$; this is in line with $g_{\mathrm{ff}}\in [1.0, 1.5]$ for typical astrospheric temperatures (Baalmann et al. 2021).

Unlike in the heliosphere or astrospheres around other cool stars, gas and dust are coupled in the astrospheres around OB stars (Henney and Arthur 2019b). Therefore, even astrosphere models that do not include dust can be projected in dust observables. To this effect, Henney and Arthur (2019c) have evaluated the emissivities of a variety of grain sizes and species around OB stars of different spectral types, as generated by the plasma physics code cloudy (Ferland et al. 2017). The emission curves for these different scenarios have been found to align very well, allowing a direct relationship between the emissivity of dust at $70~\upmu \text{m}$ and the radiation field, dependent only on the host star’s luminosity, its distance to the emitting dust, and the plasma density (Henney and Arthur 2019c; Baalmann et al. 2021).

Figure 14 shows projections of an astrosphere at different inclinations in $\text{H}\upalpha$, free-free radiation, and 70 μm dust emission, calculated as stated above. The inclination⁴ is varied from $90^{\circ}$ (viewing the astrosphere head-on) to $0^{\circ}$ (viewing the ecliptic plane face-on) in steps of $15^{\circ}$. The astrosphere around λ Cephei, a runaway blue supergiant of spectral type O6.5 at a distance of $617~\text{pc}$ to Earth, has been chosen as a prototype for an astrosphere around massive stars due to its large and visible bow shock (see Baalmann et al. 2021, for more detail on the simulation). Because of a comparably weak interstellar magnetic field, the astrosphere is symmetric about its central axis, which is why only the inclination was varied for the projections. As can be seen in Fig. 14, no astrospheric structure is observable in the head-on view of the astrosphere, where the relative motion between the star and the surrounding ISM aligns with the observer’s line of sight. The more the ecliptic plane aligns with the observer’s image plane (from left to right), the more apparent the astrospheric structure becomes, forming the bright arc that typically is identified as the observed BS. All three calculated types of emission come predominantly from the outer astrosheath, the volume between the bow shock and the astropause, where the density is high and the temperature is comparably low. The latter is explained by the fact that the $\text{H}\upalpha$ radiance is proportional to the square of the density and to the effective recombination rate coefficient (cf. Eq. (8)), which decreases with increasing temperatures. Thus, the outer astrosheath is exceptionally bright in $\text{H}\upalpha$. Free-free emission, which is also proportional to the square of the density but also to the square root of the temperature (cf. Eq. (22)), is brightest close to the astropause close to the star. Similarly, the infrared dust emission is brightest in the region of the outer astrosheath closest to the star because it depends not only on the density but also on the distance from the star.

Fig. 14

Projections in $\text{H}\upalpha$ (top row), free-free radiation (centre row), $70~\upmu \text{m}$ dust emission (bottom row) of the astrosphere of λ Cephei at a distance of $617~\text{pc}$ at ecliptical inclinations decreasing from $90^{\circ}$ to $0^{\circ}$ in steps of $15^{\circ}$ (from left to right). The ISM inflow comes from the left at $0^{\circ}$. $\text{H}\upalpha$ radiance and free-free radiance in $\text{erg}\,\text{s}^{-1}\,\text{cm}^{-2}\,\text{sr}^{-1}$, spectral dust radiance in $\text{Jy}\,\text{sr}^{-1}$, angular extent in degrees. After Baalmann (2021)

While the models mentioned above and their projections are a faithful representation of what an astrosphere in a perfectly homogeneous ISM with a perfectly symmetric outflow of its stellar wind would be, they fail to accurately reproduce the distorted shapes of most observed bow shocks (see, e.g., Cox et al. 2012). With reasonable insight into the specific stellar surroundings, it is possible to accurately model the astrosphere of a particular star (see, e.g., Gvaramadze et al. 2018); however, the necessary insight is often not available. It is, therefore, useful to apply simple perturbations to an astrospheric model and to examine the resulting distortions of the synthetic observations in order to understand the origins of similar distortions in authentic images (see, e.g., Baalmann et al. 2021).

Examples of the projections of perturbed astrospheres are given in Fig. 15. The unperturbed model is identical to that of Fig. 14; the geometric orientation is identical to the right column of that figure. The left two columns of Fig. 15 show the same perturbed model at different times, a simple perturbation of the inflow speed that causes a cavity in the ISM. Because the $\text{H}\upalpha$ radiance and the dust emissivity are sensitive to different plasma properties, the locations of their emission maxima do not align. Therefore, a comparable genuine observation may be misidentified as a detached bow shock even though gas and dust are coupled in this model. The more severe distortion of the right column stems from a much more massive perturbation of the density (i.e., a clump of dense material), generating a filament in $\text{H}\upalpha$ and a spiral-like structure in the infrared. Similar genuine observations may give the impression that relative speed between the star and its environment may point to the bottom left even though there is no vertical component to its orientation (see Baalmann et al. 2021, for more information).

Fig. 15

Projections in $\text{H}\upalpha$ (top row) and 70 μm dust emission (bottom row) of different perturbations of the astrosphere of λ Cephei at a distance of $617~\text{pc}$, viewing the ecliptic plane face-on. The perturbed ISM inflow comes from the left. $\text{H}\upalpha$ radiance in $\text{erg}\,\text{s}^{-1}\,\text{cm}^{-2}\,\text{sr}^{-1}$, spectral dust radiance in $\text{Jy}\,\text{sr}^{-1}$, angular extent in degrees. After Baalmann et al. (2021)

5.4 Filamentary Structures of Astrospheres of Runaway Stars

One of the most effective ways to observe and study astrospheres of massive wind-blowing stars are observations in the infrared where dust particles absorb and radiate. Many hundreds of new astrospheres were revealed recently with new high-resolution telescopes like the Spitzer Space Telescope, the Wide-field Infrared Survey Explorer, and the Herschel Space Observatory. Some of them have a quite specific filamentary structure (see Fig. 16), the nature of which, so far, is not well understood. Recently, Katushkina et al. (2017) proposed a physical mechanism possibly responsible for the origin of the filamentary structure of astrospheres around runaway stars by additionally modeling the spatial distribution of the interstellar dust in the astrospheres on top of the 3d kinematic mhd model results. The modeling has been performed for different magnitudes and directions of the interstellar magnetic field perpendicular and parallel to the velocity vector of the interstellar flow. It has been shown that the alternating minima and maxima of the dust density occur between the astrospheric BS and the AP due to periodical gyromotion of the dust grains. These filamentary structures appear in the model for particles with a gyration period comparable to the characteristic time of the dust motion between the BS and the AP.

Fig. 16

Spitzer 24 μm (left panel) and WISE 22 μm (middle and right panels) images of astrospheres associated with the three early B stars  Cas,  Car and  CMa, respectively. The orientation of the images is the same. Figure 2 from Katushkina et al. (2017)

In the subsequent paper by Katushkina et al. (2018), the 3d kinematic mhd model has been applied to study the astrosphere of  Cas, to produce a synthetic map of the thermal dust emission at $24~\upmu \text{m}$ and explain its observed filaments. Results of the modeling efforts are shown in Fig. 17, demonstrating how different maps could depend on the type of grains and their size distribution. It has been found that distinct filaments would appear in the emission map only if quite large (μm-sized) dust grains are prevalent in the local ISM, while it is not observable in the warm phase of the ISM continuous power-law size distribution of dust because individual filaments merge due to the influence of small grains.

Fig. 17

Examples of the model results obtained for two narrow ranges of grain radii (SO1: $[1.3, 1.7]~\upmu \text{m}$ and SO2: $[1.8, 2.2]~\upmu \text{m}$). Plots (a) and (d) present dust number density (in dimensionless units). Plots (b), (c), (e), and (f) present intensity maps for graphite and carbon. The termination shocks, astropause and the bow shock are shown by white lines. Figure 11 from Katushkina et al. (2018)

The model with large (1.3–2.2 μm) graphite and pure carbon dust grains reproduces the observational data quite well if the temperature of these grains in the region where the filaments are formed is about 40 and $75~\text{K}$, respectively. Comparison of the observed distance from the star to the brightest arc in the astrosphere with the model results allows an estimation of the ISM number density in the order of $3\text{--}11~\text{cm}^{-3}$. The local interstellar magnetic field strength is also constrained, and it should reach $18\text{--}35~\upmu \text{G}$, exceeding the typical field strength in the warm phase of the LISM.

This mechanism for the origin of the filamentary structure of astrospheres is not unique, and various processes might produce the observed filaments. For example, they could arise because of the rippling effect in radiative shocks, time-dependent wind velocity variations (Decin et al. 2006) or might be caused by AP instabilities.

5.5 GCR Modulation in the Astrospheres of Massive Stars

Based on the procedure discussed in Sect. 4.3, Fig. 18 shows the computed and normalized differential intensities for 1 TeV (black line), 1 GeV (purple line), and 10 MeV (blue line) of a massive O-type star. The results are normalized to 1.0 at the outer boundary to assure a direct comparison. Radiative cooling is included, and an ISM magnetic field of 0.3 nT is assumed in the astrospheric model. As can be seen, particles with an energy of 1 TeV and above show no modulation, while the computations for 1 GeV show modulation by a factor of four. In the case of particles with energies of 10 MeV, the interstellar differential intensity is reduced by a factor of five, reflecting the energy dependence of the cosmic ray modulation in massive star astrospheres.

Fig. 18

The computed and normalized differential intensity as a function of radial distance for energies of 1 TeV (black line), 1 GeV (purple line), and 10 MeV (blue line), normalized for a direct comparison. Radiative cooling is included while a ISM magnetic field of 3 μG is assumed in the astrospheric model

Similar results have been published in Scherer et al. (2015) using the energy dependence for the diffusion coefficients from Büsching and Potgieter (2008). However, it was found that even particles with energy as high as 100 TeV could be modulated.

6 Astrospheres of Relativistic Objects: Pulsar Wind Nebulae

Astrospheres of rapidly spinning, strongly magnetized neutron stars (pulsars) are fascinating objects, vastly different in physical, morphological, and spectral appearance. These astrospheres are powered at the expense of the rotational energy $E_{rot}$ of the pulsars, known to be the fast rotators with very stable rotation periods. If the period $P$ and its time derivative $\dot{P}$ are known from observations, the rotation energy loss (spin-down power) is given by $\dot{E}_{rot} = 4\pi ^{2} I \dot{P}/P^{3} \simeq 3.6\times 10^{38} \dot{P}_{-12}\cdot P_{10}^{-3}~\text{erg}\cdot \text{s}^{-1}$ and ranges between $5\times 10^{38}~\text{erg}\cdot \text{s}^{-1}$ for the young and energetic Crab pulsar with $P=33~\text{ms}$ and $\dot{P}=4.2\times 10^{-13}$, down to a few $10^{28}~\text{erg}\cdot \text{s}^{-1}$ for the weakest pulsars known today (here $I\simeq 1.1\times 10^{45}~\text{g}\cdot \text{cm}^{2}$ is a moment of inertia of a fiducial pulsar of mass $1.4~\text{M}_{\odot}$ and radius 10 km, $P_{10}$ is the period scaled in 10 milliseconds, and $\dot{P}_{-12}$ is a period derivative scaled in $10^{-12}~\text{s/s}$). To feed a prominent astrosphere, appearing in X-rays as a bright synchrotron nebula around the star, pulsars have to be powerful enough: $\dot{E}_{rot} \gtrsim 4\times 10^{36}~\text{erg/s}$ (e.g Gaensler and Slane 2006). However, a high spin-down is not a guarantee: sometimes quite powerful pulsars may not have an observable X-ray nebula at all or have one, but with very low X-ray efficiency (e.g. Chevalier 2000). The vast majority of the rotational energy losses goes into the pulsar wind. Significant progress in studying pulsar wind nebulae was achieved in the last two decades with the Chandra X-ray observatory which provided sensitive X-ray observations with the sub-arcsecond spatial resolution of a few dozen PWNe.

At variance with the solar wind, the pulsar wind consists predominantly of electrons and positron $(e^{\pm})$, relativistic, and strongly magnetized at its base. Being supersonic and magnetized, the wind decelerates at some termination surface where demagnetized plasma within a striped sector of the wind likely forms a TS. The TS can be strongly aspherical since the wind energy flux is latitudinally anisotropic (though the wind is asymptotically radial). The flux scales as $\sin ^{2} \theta$, being maximum in the region of the rotational equator of the pulsar and decreasing progressively toward the polar axis, from which the colatitude $\theta$ is counted (Lyubarsky 2002; Bogovalov 1999; Cerutti and Giacinti 2020). The shock shape of such a wind resembles a butterfly in its meridional cross-section, as can be seen in Fig. 19. In 3D, the shock represents a flattened spheroid, symmetrical about the axis of the pulsar’s rotation; its working surface bulges further out from the pulsar in the plane of the rotational equator and dives strongly towards the pulsar at the axis, creating a kind of cusps.

Fig. 19

Topology of the toroidal magnetic field $B_{\varphi}$ in a single-torus (left), and double-torus (right) MHD models of pulsar wind nebulae, as they appear in our setup. Shown are the poloidal slices of the nebulae, with the z-axis aligned with the rotational axis of the pulsar. The $B_{\varphi}$ strength [μG] is color-coded in the bars above the panels. Left panel: Single-torus PWN model with $\alpha =45^{\circ}$, $\sigma _{0} = 3$, evolved at low numerical viscosity from the beginning of the simulation. Right panel: The douple-torus PWN model with $\alpha = 80^{\circ}$, $\sigma _{0} = 0.03$, evolved with a high numerical viscosity upto entering a self-similar expansion, and with a low viscosity afterwards. In both panels, the black contour delineates the approximate position of the wind, terminations shock (here $\gamma = 9$)

Downstream the TS, the wind slows down, thermalizes, produces high energy non-thermal particle population, and inflates a plasma bubble, where charged wind leptons $(e^{\pm})$ begin to gyrate in the local magnetic field, emitting continuous non-thermal radiation (synchrotron and inverse-Compton) (Kennel and Coroniti 1984; Arons 2012; Sironi and Spitkovsky 2011). This radiation around the pulsar is bright enough to be seen as a bright nebulosity, the pulsar wind nebula (PWN). In combination with the wind, the PWN forms the pulsar’s astrosphere. Its properties cannot be studied without addressing the physics of the magnetospheres of pulsars, where the most $e^{\pm}$-pairs are self-produced via the cascade processes in the superstrong magnetic field of the neutron star. In pulsars with prominent PWNe, the inferred surface magnetic fields range from $B_{p}\gtrsim 10^{4}~\text{T}$ in millisecond pulsars and up to $10^{11}~\text{T}$ in magnetars (e.g. Gaensler and Slane 2006). These estimates presume a dipolar field geometry of the pulsar, in which case $B_{s} = 3.2\times 10^{19} (P\dot{P})^{1/2} \simeq 3.2\times 10^{8} (P_{10} \cdot \dot{P}_{-12})^{1/2}~\text{T}$, where in the right-hand equality $P$ is scaled in 10 ms, and $\dot{P}$ in $10^{-12}~\text{s/s}$. Recently, the particle-in-cell (PIC) simulations of pulsar magnetospheres have been extended out to a distance of $50R_{LC}$ from the pulsar (Cerutti et al. 2020), with $R_{LC}=cP/2\pi$ being the radius of the ‘light cylinder’ – an imaginary cylindrical surface centered on the spin axis of the pulsar. At this surface, the corotating speed approaches the speed of light and the magnetosphere gets open, and the pulsar wind sets off.

The $e^{\pm}$ plasma carries along an frozen-in magnetic field. As the plasma trails the pulsar’s rotation, at a large distance from the pulsar, the field becomes asymptotically toroidal $(B\equiv B_{\varphi})$. Around the rotational equator, $B_{\varphi}$ changes its polarity twice per rotational period, since the magnetic moment $\boldsymbol{\mu}$ of the star usually makes some angle $\alpha$ to its rotational axis $\boldsymbol{\Omega}$. So the equatorial sector of the wind, subtended by the angle $\theta = \pi /2 \pm \alpha$, is filled by narrow stripes of alternating magnetic polarity separated by an oscillating current sheet. As the striped wind approaches the TS and compresses at it, its alternating field is forced to reconnect and dissipate almost entirely (Lyubarsky 2003; Cerutti and Giacinti 2020). So the equatorial $e^{\pm}$ outflow enters the nebula, being barely magnetized (and subsonic). Away from the striped zone, the wind is well described by a rotating splitted monopole-like solution (Michel 1973; Bogovalov 1999). The unstriped wind outside the sector $\pm \alpha$ enters the nebula at higher latitudes across the oblique parts of the TS and remains strongly magnetized (and relativistic withal). As a result, in the inner nebula just downstream the shock, the strength of frozen in toroidal $B_{\varphi}$ is expected to peak at mid-latitudes and drop towards the pole and equator.

The energy of the wind is shared between the kinetic and field components. Their partitions can be expressed in terms of a magnetization parameter $\sigma$. The parameter represents the magnetic to particle energy densities ratio in the wind, $\sigma =B^{2}/(4\pi m_{e} n_{e} c^{2} \gamma ^{2})$, with the plasma’s co-moving density $n_{e}$ and field $B$ and the Lorenz factor $\gamma$ of the bulk outflow. The wind changes its magnetization from $\sigma \gg 1$ at the light cylinder, near which the wind originates and accelerates, to $\sigma \lesssim 1$ at the TS. Taking the partition into account, the power supplied by the pulsar to the wind is $\dot{E}_{rot} = k \dot{N}_{\mathrm{GJ}}\, m_{e} c^{2} \,\gamma (1 + \sigma )$, with $\dot{N}_{\mathrm{GJ}} \simeq 3\times 10^{30} \mu _{30}/P^{2}$ being the Goldreich-Joulian rate – the number of electrons extracted per second from the surface of the star by the surface electric field ($P$ is in seconds, and $\boldsymbol{\mu}$ is in $10^{30}~\text{G}\cdot \text{cm}^{3}$). Further, the multiplicity $k$ is the number of $e^{\pm}$ pairs self-generated in the cascade processes in the magnetosphere, per one extracted electron.

The basic properties of the pulsar wind – its initial magnetization $\sigma _{0}$, density $n_{e}$, and an opening angle $2\alpha$ of the striped zone – is difficult to deduce from the direct observations, even in nearby pulsars, since the wind is cold and radiatively inefficient until it terminates at the shock (e.g. Kirk et al. 2002; Amato 2020). In the famous Crab Nebula, for instance, the wind zone appears in X-rays as a markedly underluminous region inside the bright ring, which is plausibly identified with the TS (Weisskopf et al. 2000). At present, the primary wind properties, both at origin and in transit to the shock, can be guessed only indirectly – either via PIC modeling from first principles (e.g. Cerutti and Giacinti 2020), or via spectral and morphological modeling of the observed properties of the radiation, which the shocked wind gives off in PWNe (e.g. Komissarov and Lyubarsky 2004; Del Zanna et al. 2004; Porth et al. 2017; Amato 2020).

6.1 On the Diversity of PWNe Environments

Like the astrospheres of cool stars, pulsar astrospheres can develop in a wide variety of environments. For the $\text{few}\times 10^{4}$ years, a PWN usually expands inside the remnant of the supernova that gave birth to its parent pulsar. Later, the remnants can dissolve into the ISM, or a pulsar can leave the remnant. From this moment on, the PWN develops surrounded by the ISM.

Typical velocities that pulsars acquire at birth are hundreds of kilometers per second, but in some cases, the pulsars can make up to $\sim 1000~\text{km/s}$. These speeds are much higher than the typical sonic speed in the ISM, $c_{s} \sim 10\text{--}100~\text{km/s}$. So the motion of pulsars in the ISM can be strongly supersonic, in which case their PWNe become bullet-shaped and drive a bow shock ahead. For a recent review on the bow-shock PWNe with observational examples, see Kargaltsev et al. (2017), Bykov et al. (2017), Bucciantini et al. (2020). 3D models of the BS structure can be found in Vigelius et al. (2007) in a non-relativistic HD approach and by Olmi and Bucciantini (2019) in a relativistic MHD approach.

The internal structure of the bow-shock PWNe in MHD models has much in common with that of the heliosphere displayed in the sketch in Fig. 1. The TS of the pulsar wind is separated from the BS of the nebula by the contact discontinuity. A stand-off distance of the apex of the bow $r_{\mathrm{BS}} \sim 2\times 10^{16} \dot{E}^{1/2}_{35} n_{\mathrm{ism}}^{-1/2} V^{-1}_{200}~\text{cm}$ is determined by the pressure balance between the pulsar wind and the ISM. Here $\dot{E}_{35}$ is the pulsar spin-down power measured in units of $10^{35}~\text{erg}\,\text{s}^{-1}$, $V_{200}$ is the pulsar velocity relative to the ambient medium measured in units of $200~\text{km}\,\text{s}^{-1}$. In the tailward direction, the TS acquires a form of the Mach disk – an extended segment of the shock working surface which is nearly normal to the radial streamlines of the wind; there, the wind can slow down to subsonic velocities. In strongly supersonic PWNe, the Mach disk forms much farther from the pulsar and is followed by a sequence of weak alternating oblique and normal shocks – Prandtl-Meyer expansion-compression waves and the secondary Mach disks. This wave pattern develops as the post Mach disk outflow tries to balance its pressure with the surroundings since it enters the nebula, being overcompressed with respect to the neighboring streams that have crossed the oblique side surface of the TS.

As the pulsar wind is strongly magnetized and non-isotropic, the bow-shock PWNe observables may differ considerably from those predicted by simplified analytical and numerical models based on HD and/or 2D approaches. Varying combinations of magnetic stress and wind geometry, as well as the ISM density gradients, result in a fairly rich zoo of PWNe – with different wideness of the heads, with filled-in or edge-brightened tails, sometimes with a few jet and tail-like structures (Pavan et al. 2016; Posselt et al. 2017). The recent 3D rmhd models of bow-shock PWNe (Olmi and Bucciantini 2019; Barkov et al. 2019; Bucciantini et al. 2020) show fairly good agreements with observations. The observed X-ray emission is non-thermal, and therefore kinetic treatment of the spectral evolution of the synchrotron emitting electrons and positrons (see, e.g., Bykov et al. 2017) is necessary to disentangle the wind parameters, which determine both the inner PWN structure and the fascinating extended synchrotron structures like those observed by Pavan et al. (2016) in the Lighthouse nebula. Moreover, in some cases, the high efficiency of radiation and magnetic field amplification by kinetic instabilities due to the non-thermal population may affect the PWN dynamics beyond the MHD models.

Astrospheres of relatively young and slowly moving pulsars that have not yet escaped the remnants develop within the matter ejected by supernovae. The state of ejecta depends on the evolutionary stage of the remnant and the history of the mass loss of the massive progenitor on a pre-supernova stage. The properties of ejecta change over time as the remnant evolves and the pulsar traverses the remnant. As a result, spectro-morphological properties of astrospheres of pulsars also change over time. Although to a lesser extent than in bow-shock nebulae, in slowly moving nebulae the structure is modulated by the external medium as well. The modulation depends on the inhomogeneity of the ejecta, on its temperature and velocity profiles, not to mention its motion relative to the pulsar. The best resolved PWNe in supernova remnants – the Crab and Vela nebulae – can serve as brief illustrations of the variety of external conditions for astrospheres of pulsars within the remnants.

The supernova that gave birth to the Crab pulsar likely exploded in low-density surroundings since neither forward nor reverse shocks were found in the remnant. The Crab’s astrosphere evolves in a freely expanding ejecta with a typical velocity profile $v_{ej} \propto r$. At the pulsar’s location, the dense ejecta expands slower than the low-density PWN inflates, and the boundary between these two becomes the subject to the Rayleigh-Taylor (RT) instability. The expanding RT bubbles sweep up the local ejecta, and some of it becomes compressed into the thin shell-like structures. These structures are then photoionized by hard radiation of the nebula and appear around it as a network of thin, optically bright filaments (Hester 2008). About $\sim 26\%$ of its current energy losses the Crab pulsar converts into the synchrotron radiation of its PWN, a comparable fraction of energy goes into the work done on the acceleration of the filaments, and some energy is spent on the acceleration of nonthermal particles (Hester 2008). At present, the Crab nebula is known as the most efficient accelerator in Galaxy, capable of accelerating cosmic rays up to PeV energies.

The parent supernova of the Vela pulsar exploded in a relatively dense surrounding with a noticeable density gradient. Subsequent asymmetric expansion of the ejecta led to an asymmetric reverse shock being launched into the remnant. The shock swept over the Vela pulsar and displaced its original astrosphere, which is now observed in radio, X-rays, and gamma-rays as an extended relic Vela X nebula. The reverse shock heated the ejecta and initiated behind a weakly supersonic flow of Mach number $\simeq 1.3$. Since then, the Vela pulsar has inflated a new relativistic X-ray PWN which we now observe as the Vela pulsar wind nebula (Helfand et al. 2001; Pavlov et al. 2001; Chevalier and Reynolds 2011). It is much younger than the pulsar itself and originated in conditions different from the Crab nebula, since it developed from the very beginning in an oncoming flow. Relative to this flow, the new Vela nebula turns out to be transonic, even though the proper (physical) velocity of the Vela pulsar ($\simeq 80~\text{km/s}$) is much slower than the sonic speed in the hot ejecta at the pulsar’s location ($c_{s} \simeq 560~\text{km/s}$; Chevalier and Reynolds 2011).

6.2 On Morphology and Magnetic Field Topology in PWNe

The vast difference of PWNe in spectra and morphologies indicates the sensitivity of their observables to the differences in intrinsic parameters of pulsars and their environments – the pulsar’s spin-down $\dot{E}_{rot}$ and magnetic inclination $\alpha$, the magnetosphere’s efficiency $k$ in pair production, the wind’s flux anisotropy $F(\theta )$ and initial magnetization $\sigma _{0}$, etc. The ambient matter, in turn, can be characterized by density, temperature, magnetic field, and gas velocity relative to the pulsar. As the number of these parameters is fairly large, and some of them have a complex dependence on (co)latitude $\theta$, spectro-morphological models of PWNe often suffers from degeneracy, especially when nebulae are heavily distorted by the ram pressure of an oncoming flow (for examples of degeneracy, see, e.g., Vigelius et al. 2007; Bühler and Giomi 2016). For the nebulae moving relatively slowly with respect to their surroundings, the models are less susceptible to this problem provided the motion is properly accounted (Ponomaryov et al. 2020). The more peculiar the nebulae morphology and the finer it is resolved in space, the better the chance to break the degeneracy and thus get insight into the fundamental physics of pulsars astrosphere.

Many of PWNe in slow relative motion are visible in X-rays as jet-torus structures (see, e.g., Reynolds et al. 2017). A bright X-ray torus is usually located in the nebula’s equatorial plane (the pulsar’s rotational equator), while a jet and a counter-jet originate in the TS funnels and stretch along the symmetry axis of the torus. The toroidally shaped X-ray emission is a consequence of the interplay of the latitudinal anisotropy of the energy flux of the wind and wind’s toroidal field (e.g. Lyubarsky 2002). The jets form due to the hoop stress of the toroidal field, which is especially strong within the shock funnels, where the strongly magnetized plasma of the unstriped wind enters. The Crab and Vela nebulae mentioned above are also jet-torus objects, yet their X-ray manifestations are markedly different. The Crab nebula is shaped as a single torus (Weisskopf et al. 2000), while the Vela PWN exhibits a peculiar double-torus structure, as shown in Fig. 20 (Helfand et al. 2001; Pavlov et al. 2001). Both objects have been thoroughly studied for decades. Their multiband spectra are well established, their structures are finely resolved (down to a few light days), and their dynamics are known down to a week time scale. Both nebulae are further used as archetypes for visualizing in Fig. 19 the characteristic topology of the toroidal magnetic field in the pulsars astrospheres.

Fig. 20

Left panel: Chandra X-ray image (0.5–8 keV) of the Vela pulsar wind nebula; the first observations were reported by Helfand et al. (2001) and Pavlov et al. (2001). A bright double-torus inner nebula is visible, with two jets and with faint diffuse emission around (outlined in red), extending to south-southwest. The Vela pulsar position is marked with a cross; the pulsars’s proper velocity direction makes $\simeq 8^{\circ}\pm 5^{\circ}$ west of the PNW’s northwest jet direction. The contact discontinuity to the north from the pulsar (implied from the radio morphology, see for details Chevalier and Reynolds 2011; Bykov et al. 2017) is delineated in black. The apparent morphology of the diffuse emission can be understood in the picture where the nebula meets the large-scale transonic flow of Mach number $\sim 1.3$ from the north. Right panel: Synthetic synchrotron X-ray map of the PWN model with $\alpha =80^{\circ}$ and $\sigma _{0}=0.03$ (see text for details). The model is shown under the aspect of how the Vela nebula is seen from Earth: the northwest end of the PWN axis makes $\simeq 120^{\circ}$ to the line of sight and $\simeq 310^{\circ}$ east of north. The model is inflated within a transonic ambient matter flow of Mach number 1.3; the arrow shows the direction of the oncoming flow adopted in the model, as seen in the pulsar’s rest frame

The left panel of Fig. 19 depicts a Crab-like nebula. Shown is the field topology in the axisymmetric ideal rmhd model of the pwn with $\alpha =45^{\circ}$ and $\sigma _{0} = 3$ (close to $\alpha =45^{\circ}$ and $\sigma _{0}=1$ declared as the best for reproducing the X-ray morphology of the Crab nebula, see ref. in Amato 2000). The model is run at low numerical viscosity from the beginning. In this case, fast and highly magnetized outflows of opposite magnetic polarity, running along the arched part of the TS (black line) in the upper and lower hemispheres of the nebula, are focused towards the equatorial plane. There they twist and tangle so that the equatorial belt turns out to be the most turbulent and most magnetized (besides the jets) region within the nebula. When synthetic synchrotron X-ray emission is computed from this $\mathbf {B}$-field topology according to commonly used recipes, a thick toroidal emission structure comes about on the X-ray map of the nebula (e.g. Porth et al. 2017).

The right panel of Fig. 19 depicts a Vela-like nebula.⁵ There we again show the field topology based on a axisymmetric ideal rmhd model. The model stands for the nebula of an almost orthogonal rotator (pulsar) with a weakly magnetized wind: $\alpha = 80^{\circ}$, $\sigma _{0}=0.03$. The wide opening angle of the striped ector of the wind ($2\alpha$) allows the stripes to fill most of the wind volume and effectively consume the Poynting flux due to magnetic reconnection. As a consequence, the post-shock outflow turns out to be barely magnetized at the equatorial latitudes of the nebula. The low value of $\sigma _{0}$ adopted in the model helps to dampen the shock perturbations, which are the driving factor of the strong turbulence in the shock downstream (Camus et al. 2009; Porth et al. 2017). In the model the shock perturbations are also suppressed by a high numerical viscosity at the stage before the entry into a self-similar expansion regime. The right panel of Fig. 19 illustrates that the model creates a nebula with a barely magnetized equatorial belt and two highly-magnetized, persistent toroidal vortices flanking this belt. If synthetic synchrotron radiation were superimposed on this field topology, the X-ray map of the nebula would show two bright coaxial tori with an underluminous area in between. Alternatively one can get the double-torus structure typical for Vela-like PWNe by imposing the external plasma flow around the pulsar instead of somewhat artificial numerical approach which is using a high numerical viscosity at the initial stage of the inflation. The synchrotron map of a simulated nebula with $\alpha = 80^{\circ}$ and $\sigma _{0}=0.03$, inflated with a low numerical viscosity in the external flow of Mach number 1.3, is shown on the right panel of Fig. 20. The external flow smoothens the plasma inhomogeneities at the periphery of the compact PWN and pulls some of the plasma with magnetic vortices out from the nebula. This prevents the turbulence from building up and backreacting on the shock wave geometry. In this way, the flow can calm the shock and contribute to the formation of a low-magnetized equatorial region and thus a nebula with a double-torus structure. For a comparison to the simulated double-torus nebula the observed X-ray image of the Vela nebula obtained by the Chandra telescope (Pavlov et al. 2001) is shown in the left panel of Fig. 20.

Some of the cosmic ray spectral features revealed recently in the GeV-TeV regime with space-borne and ground-based detectors, can be produced by bow shock astrospheres. Bow shock wind nebulae can be accelerators of high energy cosmic rays by the Fermi type mechanism in the converging and colliding flows (see, e.g., Bykov et al. 2017; Pittard et al. 2020). It is exceptionally efficient in the case of BS pulsar wind nebula in which the wind is relativistic.

The Hubble Space Telescope and Chandra detected an extended astrosphere of a scale size about $2\times 10^{16}~\text{cm}$ (see Rangelov et al. 2016, and the references therein) created by the nearest old recycled millisecond pulsar PSR J0437-4715. This source can be responsible for both the enhancement of the positron fraction above a few GeV detected by PAMELA and AMS-02 spectrometers and for the TeV range spectral breaks in the cosmic ray fluxes observed with Fermi, NUCLEON, CALET, DAMPE and by the Cherenkov atmospheric telescopes H.E.S.S. (the High Energy Stereoscopic System) and VERITAS (the Very Energetic Radiation Imaging Telescope Array System) (Bykov et al. 2019; Rankin et al. 2022). Moreover, the interaction of relativistic flows produced by a pulsar wind with the equatorial disk of a Be star wind in gamma-ray binaries like PSR J2032+4127 can provide bright flares of PeV regime photons and neutrinos (Bykov et al. 2021).

7 Final Remarks

Over the past decades, we have collected a large amount of information on the Sun and the dominating physical processes like solar flares, coronal mass ejections, and the solar wind within the heliosphere. Being formed by the interaction of the solar wind (and radiation) with the interstellar gas, great detail is known about both quantities. However, our present knowledge only gives insight into how the Sun is behaving now. To understand the heliosphere’s past and future, we rely on detailed observations of other stars at different evolutionary stages. The latter, unfortunately, is challenging: while stars are moving at supersonic speed drive stellar bow shocks that are observable with the help of, e.g., Herschel, H.E.S.S. or the Aristarchos telescope, astrospheres of smaller cool stars, in particular, have not been observable so far due to instrumental restrictions.

However, with the help of state-of-the-art numerical codes like cronos, the kinematic mhd model, pluto, and 3d rmhd pwn it is possible to model the interaction of the stellar wind and the interstellar medium. Their output can be compared to observations and – in the long run – might help to predict the mandatory observations to study stellar astrospheres in more detail. Thus, computational and observational astronomy can support each other to understand stars and their environments better.

This review focused on the astrospheres produced by stellar-mass objects of different genesis: from the astrospheres of cool diverse M-stars, the bow-shocks of hotter runaway stars to relativistic pulsar wind nebulae, which are as diverse as their stars. In addition, the winds of young massive stars and supernovae are producing structures hundreds of parsecs large around the starforming regions in galaxies (Krause et al. 2020), which are efficient sources of cosmic rays, high energy neutrinos, and non-thermal radiation (see, e.g. Bykov et al. 2020, for a recent review). However, there are also very impressive astrospheres produced by massive black holes or starburst regions produced by a huge energy release. Predehl et al. (2020) reported recently an eROSITA SRG finding of soft-X-ray-emitting bubbles extending both above and below the Galactic center to 14 kpcs. The estimated energy needed to create the large-scale structures is about $10^{56}~\text{ergs}$.

We hope to have shown in this review that significant progress in understanding the formation and evolution of stellar astrospheres and bow shocks has been made. However, modeling and observing astrospheres of, in particular, planet-hosting stars will be of utmost importance in the upcoming years. With the newly-emerging field of studying the impact of cosmic rays on exoplanetary habitability (e.g., Herbst et al. 2019b,c, 2021a; Mesquita et al. 2021; Rodgers-Lee et al. 2021) many new exciting questions have opened up, and answering them will only be possible with the help of close collaborations between diverse research fields.