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

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 Hα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\text{H}\upalpha$\end{document} projections to 1pc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$1~\text{pc}$\end{document}, 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.


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  (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 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).
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: where ρ, 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. 1 The importance of neutrals and other species like pick-up ions is discussed in Sokół et al. (2022). continuity, momentum, and energy equations with the Euler-Maxwell equations is given by: Here, j gives the particle species, v j represents its velocity, ρ j its mass density, p j its pressure, and the total energy is with the inner energy e j . Furthermore, I gives the unit tensor,σ the stress tensor, F an external force per unit mass and volume, Q the heat flow, while S x j gives sources and sinks caused by charge exchange in the continuity (x = c), energy (x = e), and momentum (x = m) equation. ⊗ describes the dyadic/tensorial product. To account for cooling effects R L represents a cooling function. The subscript rad accounts for radiation transport of momentum and energy coupling, while CR accounts for cosmic rays. In addition, P 1 , P 2 , P 3 , and P 4 are constants, while 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 In the case of ideal HD modeling the following conditions must apply: ρ = 0, v = 0, E = 0, P 1 = P 2 = P 3 = P 4 = 0, R L = S x j = F rad = E rad = 0, as well as B = 0. In case of ideal MHD modeling the same criteria apply, however B = 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).

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. 2015Scherer et al. , 2016aKleimann 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: 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 ). The computational grid was arranged in spherical coordinates (r, ϑ, ϕ), producing a spherical polyhedron with equidistant intervals in radius r and equiangular intervals in the polar (ϑ ) and azimuthal (ϕ) 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 .
Both codes have been utilized to produce results discussed in Sects. 3, 4, and 5.

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 with F r and F g as the forces of solar radiation pressure and gravitation, respectively, f p (r, w H ) as the local Maxwellian distribution function of protons, σ HP ex (u) as the effective charge-exchange cross section with u as the relative atom-proton velocity, and ν ph = 1.67 × 10 −7 (R E /R) 2 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.

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.

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.

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. 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 AU around the Sun, including values of the number density n, the thermal pressure p therm , as well as the vectors of the flow velocity v and the magnetic flux density B. The values are normalized to r 0 = 1 AU, Furthermore, the temperature is given by T = T 0 · p therm /n, with T 0 = v 2 0 · m p /(2k B ) ≈ 28811 K, k B the Boltzmann constant, and the normalized quantities p therm and n. The total pressure p tot is composed of the thermal pressure p therm , the dynamic pressure p dyn , and the magnetic pressure p mag . The latter two are given by p mag = (B 0 B) 2 /(2μ 0 p 0 ) and p dyn = 1/(2p 0 ) · m p · n 0 · n · (v 0 v) 2 , respectively, where μ 0 = 4π · 10 −7 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 f . This result, 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 however, strongly depends on the utilized model. According to Scherer and Fichtner (2014), a BS will occur when a multi-fluid simulation is applied.
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 AU and 125 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 AU and the HP at 122 AU, while Voyager 2 passed the TS at 84 AU and the HP at 119 AU. Thus, the model results are in good agreement with the observations (see Richardson et al. 2022).

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. 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-α emission line provides a sensitive tool to measure the stellar mass-loss rates of cool stars. In stars cooler than the Sun, the Lyman-α emission line dominates the stellar UV spectrum (e.g., Linsky 2019); however, Lyman-α 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-α 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-α 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-α absorption profile to the astrospheric models, one can gather information about the nearby astrospheres (lower panels of Fig. 4).

Hα 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 .
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 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 AU would appear as a disk with a radius of φ max = 0.25 • = 15 . Therefore, a projection of the heliosphere would require a pixel grid with edge lengths 30 × 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, ( r, ϑ, ϕ), ϑ = ϕ, the lateral distance between cells varies within the model. The largest distance between model cells, ( s) max , lies at the outside of the model, ( s) max = r max − 1 2 r tan( ϕ). An example grid with r = 3 AU and ϕ = 3 • yields ( s) max ≈ 47 AU. To avoid empty pixels, the pixel grid resolution shall be coarser than the largest apparent distance between model cells, which for the example grid leads to ( φ) max ≈ 0.79 . To improve the grid resolution, it is expedient to linearly interpolate in the two angular directions by a small factor k ∈ N, typically k ≤ 5. A typical CRONOS model of the heliosphere as described in Sect. 2 simulates the number density n, temperature T and bulk velocity v of fully-ionized hydrogen as well as the magnetic field B in every model cell. With these parameters, some emission processes like bremsstrahlung, cyclotron radiation, or electronic transitions to high-energy states, i ≥ 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 Hα, arising from the electronic transition 3 → 2 of hydrogen. While Hα 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 Hα maps of the heliosphere for comparison.
The radiance of Hα of a model cell i is calculated by Here, α 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  with data from Wiese andFuhr (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 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 AU and r max = 900 AU. The model grid consists of 280 × 60 × 120 cells, resulting in radial and angular cell dimensions of ( r, ϑ, ϕ) = (3 AU, 3 • , 3 • ). 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 px per edge, which corresponds to a linear resolution of φ = 0.19 ; roughly 3 × 10 −9 sr per pixel area.
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 Hα 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 × 10 −12 erg s −1 cm −2 sr −2 surrounded by the disk-like maximum of 2.5 × 10 −12 erg s −1 cm −2 sr −2 caused by the high-emission bow wave. The small-scale arcs, most apparent in the middle panel along with the 0 • -axes, are Moiré patterns and, thus, non-physical structures.

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 AU (late-type M-dwarfs) and 0.5 AU (young M-dwarfs). According to West et al. (2004West et al. ( , 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. 2018Scheucher et al. , 2020a 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.

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 X . Thereby, L X ∝ P −2.7 rot 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 X seems to saturate at a luminosity around 10 −3 L 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 sat rot ≈ 1.2 (L bol /L ) − 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 X ∝ 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 rot ∝ t 0.6 (see, e.g., Ayres 1997). and the corresponding evolutionary tracks of the L X values in the right panels. As can be seen, in the case of a Sun-analog (1 M ; 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 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 X are almost identical until 2 Gyr (0.25 M , lower right panel).

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 κ 1 Ceti most recently were published by . 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 |B| of Sunlike stars decays with age and rotation period (see panels of Fig. 7). Statistically speaking, thereby the average large-scale magnetic field strength scales with where both provide constraints on the evolution of the average magnetic field strength.
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.

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.

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  (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 F mag = ([∇ × B] × B)/(4π) 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 ∼ 1/R 2 , ρ ∼ ρ ∞ , and p ∼ p ∞ , where V , ρ, p are the velocity, density and pressure of the stellar wind, respectively, whereas ρ ∞ and p ∞ 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 ∇ × [V × B] = 0 and assuming that the magnetic field is parallel to the stellar wind velocity vector. The solution above for R < R TS has been obtained by Parker (1958). For the subsonic region beyond the TS (R > R TS ), the following solution can be found: where θ is the polar angle counted from the stellar rotational axis (x-axis) and φ is the azimuthal angle. As can be seen, the plasma β 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 where B ϕ,E is the strength of the magnetic field at the distance of R E . The IMF and the radial component B R are assumed to be zero. The latter is not critical since B R ∼ 1/R 2 is small compared to B ϕ at large heliocentric distances. Based on Eq. (10), Korolkov and Izmodenov (2021)  Further, the value of M crit,∞ increases with the decrease of the Alfvénic Mach number that corresponds to the effective increase of the stellar magnetic field. At M ∞ = 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.

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  (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 (ϑ, ϕ) = 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.

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 ∼ 40 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 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): where f is the GCR distribution function, V is the convection velocity, κ 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) for the time-backward case, and for the change in momentum The Wiener process is given by the standard normal distribution, dw r , with a mean of zero and a standard deviation of one (e.g. Zhang 1999) 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 The time constraint is chosen so (Strauss and Effenberger 2017). Further, the radial diffusion coefficient κ rr can be written in terms of a perpendicular diffusion mean free path (MFP) λ ⊥ 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;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).
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 Mesquita et al. 2021, respectively).

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;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;. 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. Wiengarten et al. 2016;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  (1958) field (solid blue line). Note that the astrosphere of LHS 1140 is considerably smaller than the heliosphere, with a TS at ∼ 8.1 AU, and an astrosheath extending out to ∼ 28 AU. In order to estimate magnetic variances and correlation scales, we scale values observed in the heliosphere at 1 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 × 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 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., . 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 δB 2 ∼ r −2.4 , and it is assumed that both MFPs λ 2D and λ sl ∼ 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. 2014Burlaga et al. , 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;) 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.
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. 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 where R = R L k m , with R L the maximal particle gyroradius and k m = 1/λ 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 δB 2 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) where α 2 = 1/3 (Matthaeus et al. 2003), ν = 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 λ ⊥ . 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  from the theory of Bieber and Matthaeus (1997), given by with δB 2 T = δB 2 sl + δB 2 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 λ A = R 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 λ ⊥ scales as the inverse of the astrospheric magnetic field magnitude, normalised to a value λ 0 = aB E (1 AU)/B 1 , with B E and B 1 the magnitudes of the heliospheric and astrospheric magnetic fields at 1 AU, and a a free parameter, chosen here to be equal to 0.01 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 where P the particle rigidity, in units of GV, P 0 = 1 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 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 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 λ 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 L ) due to the relatively small perpendicular MFP yielded by Eq. (17) at this rigidity within ∼ 10 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.

Astrospheres of Massive Stars -The Bow-Shocks of Runaway
Massive stars (≥ 8 M ) 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.

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 ζ 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: where ρ and v are the density and speed, respectively, the index SW denotes the solar (stellar) wind, and 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Ṁ leading to the formula deduced by (Baranov and Malama 1993) or Wilkin (1996): 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;). 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 andGvaramadze 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 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 Xray (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. , 2015Kobulnicky et al. 2016Kobulnicky et al. , 2017Kobulnicky et al. , 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(Mignone et al. , 2012Vaidya et al. 2018), although Lagrangian models exist . Furthermore, full 3D MHD simulations of runaway stars have been performed recently. The results on λ-Cephei have been discussed in  and  and those on κ-Cas by . 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 K) is important (Comerón and Kaper 1998;Meyer et al. 2016Meyer et al. , 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).

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 (∼ Myr) changes than that of low-mass stars (∼ 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 mainsequence 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 . 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, α Ori (Betelgeuse), and µ Cep, raising the question of the apparent missing bow shock problem. The BS of IRC-10414 ) 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 ). Its roundness indicates its young age, i.e., Betelgeuse just underwent a stellar phase transition . A mysterious bar, almost parallel to the direction of stellar motion, probably of interstellar origin (Decin et al. 2012) 3 completes this picture.
Some massive stars ≥ 35 M 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(van Marle et al. , 2007Meyer 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(Meyer et al. , 2020b(Meyer et al. , 2021b.

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 Hα 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 where n i = n i,e = n i,H + , T i , V i , d i are the number density, temperature, volume, and distance to the observer of model cell i, respectively; φ is the linear pixel resolution; k B , m e , e, h, c are the Boltzmann constant, electron mass, elementary charge, Planck constant and vacuum speed of light, respectively; and g ff (T ) is the temperature-dependent Gaunt factor for freefree transitions, which can be numerically calculated by, e.g., van Hoof et al. (2014). An often-made assumption for the Gaunt factor isḡ ff = 1.2; this is in line with g ff ∈ [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 µ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 Hα, free-free radiation, and 70 µm dust emission, calculated as stated above. The inclination 4 is varied from 90 • (viewing the astrosphere head-on) to 0 • (viewing the ecliptic plane face-on) in steps of 15 • . The astrosphere around λ Cephei, a runaway blue supergiant of spectral type O6.5 at a distance of 617 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 Hα 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 The ISM inflow comes from the left at 0 • . Hα radiance and free-free radiance in erg s −1 cm −2 sr −1 , spectral dust radiance in Jy sr −1 , angular extent in degrees. After Baalmann (2021) bright in Hα. 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.
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., . With reasonable insight into the specific stellar surroundings, it is possible to accurately model the astrosphere of a particular star (see, e.g., ; 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 Hα 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 Hα 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).

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 hun-

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) dreds 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. , 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  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.
In the subsequent paper by , 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 µ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.
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 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-11 cm −3 . The local interstellar magnetic field strength is also constrained, and it should reach 18-35 µG, exceeding the typical field strength in the warm phase of the LISM. 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 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.

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

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 deriva-tiveṖ are known from observations, the rotation energy loss (spin-down power) is given bẏ E rot = 4π 2 IṖ /P 3 3.6 × 10 38Ṗ −12 · P −3 10 erg · s −1 and ranges between 5 × 10 38 erg · s −1 for the young and energetic Crab pulsar with P = 33 ms andṖ = 4.2 × 10 −13 , down to a few 10 28 erg · s −1 for the weakest pulsars known today (here I 1.1 × 10 45 g · cm 2 is a moment of inertia of a fiducial pulsar of mass 1.4 M and radius 10 km, P 10 is the period scaled in 10 milliseconds, andṖ −12 is a period derivative scaled in 10 −12 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:Ė rot 4 × 10 36 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 ± ), 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 θ , being maximum in the region of the rotational equator of the pulsar and decreasing progressively toward the polar axis, from which the colatitude θ 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.
Downstream the TS, the wind slows down, thermalizes, produces high energy nonthermal particle population, and inflates a plasma bubble, where charged wind leptons (e ± ) 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 ± -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 10 4 T in millisecond pulsars and up to 10 11 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 × 10 19 (PṖ ) 1/2 3.2 × 10 8 (P 10 ·Ṗ −12 ) 1/2 T, where in the right-hand equality P is scaled in 10 ms, andṖ in 10 −12 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 , with R LC = cP /2π 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 ± 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 ≡ B ϕ ). Around the rotational equator, B ϕ changes its polarity twice per rotational period, since the magnetic moment μ of the star usually makes some angle α to its rotational axis . So the equatorial sector of the wind, subtended by the angle θ = π/2 ± α, 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 ± 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 ±α 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 ϕ 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 σ . The parameter represents the magnetic to particle energy densities ratio in the wind, σ = B 2 /(4πm e n e c 2 γ 2 ), with the plasma's co-moving density n e and field B and the Lorenz factor γ of the bulk outflow. The wind changes its magnetization from σ 1 at the light cylinder, near which the wind originates and accelerates, to σ 1 at the TS. Taking the partition into account, the power supplied by the pulsar to the wind isĖ rot = kṄ GJ m e c 2 γ (1+σ ), withṄ GJ 3×10 30 μ 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 μ is in 10 30 G · cm 3 ). Further, the multiplicity k is the number of e ± pairs self-generated in the cascade processes in the magnetosphere, per one extracted electron.
The basic properties of the pulsar wind -its initial magnetization σ 0 , density n e , and an opening angle 2α 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;. 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 indirectlyeither 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;.

On the Diversity of PWNe Environments
Like the astrospheres of cool stars, pulsar astrospheres can develop in a wide variety of environments. For the few × 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 ∼ 1000 km/s. These speeds are much higher than the typical sonic speed in the ISM, c s ∼ 10-100 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 ), 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 BS ∼ 2 × 10 16Ė 1/2 35 n −1/2 ism V −1 200 cm is determined by the pressure balance between the pulsar wind and the ISM. HereĖ 35 is the pulsar spin-down power measured in units of 10 35 erg s −1 , V 200 is the pulsar velocity relative to the ambient medium measured in units of 200 km 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 PWNewith 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 bowshock 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 ∝ 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 ∼ 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 gammarays as an extended relic Vela X nebula. The reverse shock heated the ejecta and initiated behind a weakly supersonic flow of Mach number 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 ( 80 km/s) is much slower than the sonic speed in the hot ejecta at the pulsar's location (c s 560 km/s; Chevalier and Reynolds 2011).

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 environmentsthe pulsar's spin-downĖ rot and magnetic inclination α, the magnetosphere's efficiency k in pair production, the wind's flux anisotropy F (θ) and initial magnetization σ 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 θ , 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.
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 α = 45 • and σ 0 = 3 (close to α = 45 • and σ 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 The douple-torus PWN model with α = 80 • , σ 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 γ = 9) 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 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. 5 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: α = 80 • , σ 0 = 0.03. The wide opening angle of the striped ector of the wind (2α) 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 σ 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 α = 80 • and σ 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  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 8 • ± 5 • 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 ∼ 1.3 from the north. Right panel: Synthetic synchrotron X-ray map of the PWN model with α = 80 • and σ 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 120 • to the line of sight and 310 • 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 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 lowmagnetized 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 × 10 16 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;). 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).

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 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. 2019bHerbst et al. ,c, 2021aMesquita 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.
Funding Note Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.