© L. R. Baalmann et al. 2022

1 Introduction

Stellar bow shocks occur when the relative motion between a star and its environment is super(fast magneto)sonic. This speed criterion can be met either by a star at rest experiencing an ambient flow, which is referred to as the weather vane scenario (Povich et al. 2008), or, more commonly, by a runaway star moving through the interstellar medium (ISM; e.g. van Buren et al. 1995). While a wide variety of stellar objects can produce bow shocks (e.g. van Buren & McCray 1988), for example red super-giants and AGB stars (Cox et al. 2012), a significant number of bow shock observations are associated with OB stars (e.g. Kobulnicky et al. 2017, and references therein).

The interactions between a respective star’s stellar wind (SW) and its surrounding ISM are typically modelled numerically, using two-dimensional (2D) hydrodynamics (HD; e.g. Matsuda et al. 1989; Comeron & Kaper 1998; Scherer et al. 2016, and references therein), 2D magnetohydrodynamics (MHD; van Marle et al. 2014; Meyer et al. 2017), or, more recently, three-dimensional (3D) MHD (e.g. Scherer et al. 2020; Meyer et al. 2021; Fraternale et al. 2021, and references therein). From a computational perspective, the runaway and the weather vane scenarios are generally indistinguishable.

The interaction region of the SW and the ISM consists of three important features overall, all of which are MHD discontinuities (cf., e.g. Goedbloed et al. 2010, Sects. 20.2 and 20.3): the bow shock (BS), at which the ISM inflow is decelerated from faster than a characteristic speed to slower than a characteristic speed; the astropause (AP), which separates the stellar and interstellar fluids; and the termination shock (TS) and Mach disk (MD), which are the counterparts of the BS for the SW. In MHD, the characteristic speed is the fast magnetosonic speed, υ_f, whereas in HD, it is the sonic speed, υ_c. If the ISM inflow is not faster than the characteristic speed, namely, if it is sub-fast or subsonic, no BS occurs. The BS, AP, TS, and MD are all surfaces of infinitesimal thickness; the BS, TS, and MD are genuine MHD shocks, whereas the AP is a tangential discontinuity. The BS is shaped like a paraboloid; its apex is pointed towards the direction from which the inflow comes, which is called the upwind direction, and lies in front of the star. The AP is bullet-shaped; its apex lies in front of the star but downwind from the BS, that is, in between the star and the BS. The TS is shaped like a stylised tulip (see, e.g. Figs. 1 and 5) with a similar orientation and position of its apex, downwind from the AP. The MD is in good approximation shaped like a spherical cap and it serves as the counterpart to the TS in the downwind direction. Its apex is pointed downwind and it meets the TS at a roughly circular line that is, for historical reasons, known as the triple point (TP; Bleakney & Taub 1949). From the TP, another tangential discontinuity (TD), which is roughly cylindrical in shape (e.g. Scherer et al. 2016, Sect. 4.1), extends into the downwind direction.

The region between the TS & MD and the AP is referred to as the inner astrosheath, while the one between the AP and the BS is known as the outer astrosheath. Some authors, such as Zank (2015), refer to these regions as the astrosheath and the very local interstellar medium (VLISM), respectively; the outer astrosheath is occasionally also referred to as the bow shell (e.g. Henney & Arthur 2019a). What is often called a bow shock in observational studies is typically emission from the outer astrosheath (Baalmann et al. 2020). The region inside the TD is called the astrotail, though this term often more generally refers to the entire region downwind of the MD.

If the magnetic field is weak compared to HD effects, the astrospheric structure is rotationally symmetric about the axis that is parallel to the upwind and downwind directions and passes through the star, namely, the central axis. In HD models, this line is typically referred to as the stagnation line because an ISM flow along this axis would have to stagnate some distance after crossing the BS in order to preserve the rotational symmetry. In MHD a stagnation point, where the fluid comes to rest, may occur as well but generally does not lie on the central axis because flow lines are bent by the magnetic field (e.g. Scherer et al. 2020, Sect. 2.1).

Within this study the astrosphere around the runaway O-star λ Cephei has been used as a basis for the modelling procedure (cf. Scherer et al. 2020; Baalmann et al. 2020, 2021). Following the classification scheme proposed by Henney & Arthur (2019a,b,c), this astrosphere features a wind-supported BS, that is, the dynamics are determined by the MHD flows instead of the stellar radiation, and the interstellar dust is highly coupled to the plasma. This justifies the use of single-fluid non-radiative MHD to model the astrosphere, taking the stellar radiation and the interstellar dust into account only as energy gains and losses via heating and cooling effects, respectively (cf. Sect. 2.1).

The structure of this paper is as follows. Section 2 gives a brief overview of the applied methodology, Sect. 3 discusses the results, and Sect. 4 provides a summary and draws conclusions.

2 Methodology

This section gives an overview of the methodology used in this investigation. The computational model is described in Sect. 2.1 and details of this approach can be found in Scherer et al. (2020) and Baalmann (2021, Sect. 2.4). The projection method is summarised in Sect. 2.2. It was first introduced in Baalmann et al. (2020) and expanded upon in Baalmann et al. (2021). A detailed formulation is given in Baalmann (2021, Sect. 3). The setup of this investigation is described in Sect. 2.3.

2.1 Computational model

The model astrospheres were computed as described by Baalmann et al. (2021, Sect. 2.1), using the semi-discrete finite-volume code CRONOS (Kissmann et al. 2018) to solve the single-fluid ideal MHD equations: $$\frac{\partial n}{\partial t}+\nabla \cdot (nu)=0,$$(1) $$\frac{\partial }{\partial t}(mnu)+\nabla \cdot (mnu\otimes u)+\nabla p+\frac{1}{{\text{μ}}_0}B\times (\nabla \times B)=0,$$(2) $$\frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+p+\frac{1}{2{\text{μ}}_0}{\left|B\right|}^2\right)u-\frac{1}{{\text{μ}}_0}(u\cdot B)B\right]=\Gamma -\Lambda ,$$(3) $$\frac{\partial B}{\partial t}-\nabla \times (u\times B)=0,$$(4)

using an HLL Riemann solver and a second-order Runge–Kutta scheme. In the above equations n, u, p, B, and e are the number density, fluid velocity, thermal pressure, magnetic induction, and total energy density, respectively; m and μ₀ are the proton mass and vacuum permeability, respectively. The dyadic product is represented by ⊗. Two additional equations close the above system: $$e=\frac{p}{\gamma -1}+\frac{1}{2}mu{\left|u\right|}^2+\frac{1}{2{\text{μ}}_0}{\left|B\right|}^2,$$(5) $$\nabla \cdot B=0,$$(6)

assuming a mono-atomic ideal gas, leading to the polytropic index γ = 5/3. A source term was added to the right-hand side of the energy Eq. (3), incorporating gains, Г, and losses, Λ, of the total energy density by heating and cooling effects, respectively. Heating was incorporated as per Kosiński & Hanasz (2006), including photoionisation (α n²), and photoelectric heating by dust, Coulomb collisions with cosmic rays, and the dissipation of interstellar turbulence (α n). Radiative cooling was incorporated following Schure et al. (2009) within the temperature range log(T [K]) ∈ [3.8, 8.16]; for lower temperatures, cooling effects were neglected; higher temperatures did not occur. As discussed in detail in Baalmann et al. (2021, Sect. 3.1), the terms referred to as heating and cooling, which change the total energy density as per Eq. (3), have only an insignificant influence on the temperature of the plasma due to the additional energy input being continuously and self-consistently redistributed between the constituent terms of Eq. (6).

The computational grid takes the shape of a full sphere, described by spherical coordinates, (r,ϑ,φ), where r ∈ [r_SW, r_ISM] is the radial coordinate, that is, the distance from the coordinate centre (r = 0), which corresponds to the position of the star; r_SW > 0 and r_ISM ≫ r_SW are the inner and outer boundaries of the model (Table 1), at which the respective boundary conditions (top segment of Table 2) for the SW and the ISM are set. More precisely, these are kept fixed for the inner (r_SW) and the inflow-upwind half of the outer (r_ISM) spherical boundary, whereas at the outflow-downwind half of the outer boundary all values are extrapolated linearly into the boundary cells. Additionally, no material is allowed to flow back into the computational domain through the downwind half-sphere by prohibiting any inward motion of material there, u_r ≥ 0. The stellar parameters of λ Cephei are given in the bottom segment of Table 2. The angles ϑ ∈ [−90°, 90°] and φ ∈ [−180°, 180°] are the polar and the azimuthal angle, respectively. Model cells are arranged with equidistant intervals in r and equiangular intervals in ϑ, φ; thus, the spatial extents of the cells perpendicular to r as well as their volumes increase with r. The coordinates (r_i, ϑ_i, φ_i) of model cell i describe the central point of that cell. The numbers of cells in the {r, ϑ, φ}-direction are referred to as N_{{r,ϑ, φ}}; the radial and angular extents of the cell sizes are referred to as ∆{r,ϑ, φ} (Table 1).

For the 3D models the cell sizes in both angular directions are set to be identical, ∆ϑ = ∆φ, 2N_ϑ = N_φ. The angles of the homogeneous ISM fluid inflow are set to ϑ_u = 90° and φ_u = 180°, aligning the inflow with the x-axis. Thus, the x-axis points in the upwind direction, which is the direction from where the ISM inflow comes, or, conversely, in which the star moves. The direction in which the negative x-axis points is referred to as the downwind or tail direction. The SW is set up as an isotropic outflow. The homogeneous interstellar magnetic field is angled with ϑ_B and φ_B; the SW magnetic field is compliant with Parker’s spiral (Parker 1958) for the case of a purely radial, constant, monopolar field strength at the boundary surface, as implemented by Scherer et al. (2020, Sect. 3.1.1), $$B_{\text{SW}}(r)=B_0\frac{r_{\text{SW}}^2}{r^2}\left(e_r+\frac{r{\Omega }_{\text{SW}}}{u_{\text{SW}}}\sin \vartheta e_{\phi }\right),$$(7)

where B₀ is a prefactor chosen to comply with $r_{\left\{{\text{BS,AP,MD}}_{\left\{\max ,\min \right\}}\right\}}$; e_{r,φ} are the respective unit vectors, and Ω_SW = 2π/P_* is the angular frequency of stellar rotation.

For the 2D models the φ-direction is covered by a single cell, N_φ = 1, ∆φ = 360°, forcing cylindrical symmetry about the z-axis. The ISM inflow comes from the positive z-direction, ϑ_u = 0°. To maintain symmetry the interstellar magnetic field must be parallel to the z-axis, ϑ_B ∈ {0°, 180°}, and the SW magnetic field must be zero, B_SW = 0.

The simulation times of all models are normalised to the time it takes to travel 1 pc at the Alfvén speed in a medium with a number density of n = 1 cm⁻³ and a magnetic field of 1 nT, which is υ_A,0 = 1 nT · (1 cm⁻³μ₀m_p)^−1/2 ≈ 21.812 km s⁻¹, resulting in a time normalisation of $$t_{\text{pc}}=\frac{1\text{pc}}{{\upsilon }_{\text{A,0}}}\approx 44829\text{yr.}$$(8)

2.2 Projection method

In order to generate synthetic observational images, the models were projected onto a virtual sky as described by Baalmann et al. (2020). To this effect the model was shifted to a desired position, (d₀, b₀, l₀), in Galactic coordinates, where d₀ is the distance modulus of the star, b₀ its Galactic latitude, and l₀ its Galactic longitude. For all projections in this work b₀ = l₀ = 0 was chosen because these coordinates merely function as an offset of the image on the projection grid. The 2D models were extruded to cylindrically symmetric 3D models with 2N_ϑ = N_φ before projection. Because the models have the shape of a sphere, they appear as a disk with a radius of $${\varphi }_{\max }=\arcsin \left(\frac{r_{\max }}{d_0}\right),$$(9)

where r_max = r_ISM − ∆r/2 is the radius of the outermost model cell’s centre. The largest spatial distance between two adjacent model cells is $${(\Delta s)}_{\max }=\max \left\{\Delta r,2r_{\max }\sin \left(\frac{\Delta \vartheta }{2}\right)\right\},$$(10)

corresponding to the largest apparent distance between two model cells within the projection, $${(\Delta \varphi )}_{\max }=2\arctan \left(\frac{\frac{1}{2}{(\Delta s)}_{\max }}{d_0-r_{\max }}\right).$$(11)

The square projection grid of size ϕ_max × ϕ_max is divided into N × N pixels, N ∈ {2, 3}. The pixels are recursively refined by the same division until the pixel size would be smaller than (∆ϕ)_max. Because the pixel resolution, ∆ϕ, is limited by the spatial distance of the model’s outermost cells, (∆s)_max, due to the constant angular cell extent, ∆ϑ, it cannot adequately resolve the model’s inner structure. Therefore, the model is linearly interpolated in both angular coordinates by a factor of k = 5 before projection.

The model cells are assigned to the pixel in which their Galactic coordinates appear, the projection value of each pixel is computed by summing the Hα radiance of all its assigned model cells, $$L_{\Omega ,i}={\alpha }_{\text{eff,3}\to \text{2}}(T_i)h{v}_{\text{H}\alpha }{\eta }_i^2\frac{V_i}{4\pi d_i^2}\frac{1}{{(\Delta \varphi )}^2},$$(12)

where V_i and d_i are the volume and distance of model cell i, h is the Planck constant, v_Hα = 456.81 THz is the frequency of Hα photons, and α_eff,3→2(T_i) is the temperature-dependent effective recombination rate coefficient, approximated by $${\alpha }_{\text{eff,3}\to \text{2}}(T_i)\approx {\alpha }_3(T_i)+{\displaystyle \sum _{q=4}^{16}{\alpha }_q(T_i)P_{q,3},}$$(13)

where α_q(T_i) are the temperature-dependent recombination rate coefficients of state q, taken from Mao & Kaastra (2016), and P_q,3 are the branching ratios from state q through state 3, computed from data by Wiese & Fuhr (2009) as described by Baalmann et al. (2020). The highest radiance is expected to come from the outer astrosheath, where the density is highest and the temperature is comparably low (Baalmann et al. 2020).

2.3 Investigative setup

The objective of this examination was to compare otherwise identical astrospheric models with differing ISM inflow speeds, u_ISM, including a model that would not generate a BS, in either 2D or 3D. As a prototype the astrosphere of the runaway O-star λ Cephei was chosen, cf. Baalmann et al. (2020, 2021). The stellar parameters of λ Cephei are given in the bottom segment of Table 2; the initial boundary conditions of the model are listed in the top segment of the same table. The relative speed between λ Cephei and its surrounding ISM is roughly u_ISM = 80 km s⁻¹. The model astrosphere with this u_ISM is identical to the hires model presented by Baalmann et al. (2021).

The sonic and Alfvénic speeds, υ_{c,A}, for the surrounding ISM are $${\upsilon }_{\text{c,ISM}}=\sqrt{\frac{\gamma p_{\text{ISM}}}{{\rho }_{\text{ISM}}}}\approx 15.74\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1},$$(14) $${\upsilon }_{\text{A,ISM}}=\frac{B_{\text{ISM}}}{\sqrt{{\text{μ}}_0{\rho }_{\text{ISM}}}}\approx 6.58\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1},$$(15)

where p_ISM = 2n_ISMk_BT_ISM is the thermal pressure of the ISM, k_B is the Boltzmann constant, factor 2 comes from the quasi-neutrality and pressure equilibrium of protons and electrons, while ρ_ISM = (m_p + m_e)n_ISM ≈ m_pn_ISM is the mass density of the ISM. The characteristic speed for MHD shocks is the fast magnetosonic speed, $$\begin{array}{c}{\upsilon }_{\text{f,ISM}}={\sqrt{\frac{1}{2}\left({\upsilon }_{\text{A}}^1+{\upsilon }_{\text{c}}^2+\sqrt{\left({\upsilon }_{\text{A}}^2+{\upsilon }_{\text{c}}^2\right)-4{\upsilon }_{\text{A}}^2{\upsilon }_{\text{c}}^2{\cos }^2{\vartheta }_{\upsilon ,B}}\right)}|}_{\text{ISM}}\\ \in \left[\max \{{\upsilon }_{\text{c,ISM}},{\upsilon }_{\text{A,ISM}}\},\sqrt{{\upsilon }_{\text{c,ISM}}^2+{\upsilon }_{\text{A,ISM}}^2}\right]\\ \approx [15.74,17.06]km{s}^{-1},\end{array}$$(16)

which depends on the angle between the bulk velocity and the magnetic field, ϑ_υ,B. The speed of the ISM inflow speed, u_ISM, was incrementally decreased by a factor of two until it reached a subfast value, yielding u_ISM ∈ {80, 40, 20, 10} km s⁻¹. Because the size of the shock structure is dependent on, for instance, the ISM inflow speed, models with smaller u_ISM require larger outer boundaries of the simulation grid. To keep computational times low, the inner boundaries were adapted similarly. The radii of the inner and outer boundaries, r_{SW,ISM}, for the different ISM inflow speeds, as well as the radial resolutions, ∆r, are given in Table 1.

3 Results

The results of the investigation are discussed in this section. The 2D models and the 3D models are analysed in Sects. 3.1 and 3.3, respectively. In order to better understand the formation of the structures displayed in these sections, the evolution of the models is expounded at the example of a 2D model in Sect. 3.2. Synthetic observations are presented in Sect. 3.4. Further aspects are discussed by Baalmann (2021, Sect. 4.4).

3.1 2D models

Slice planes of the 2D models are presented in Fig. 1. The vertical axes are the z-axis of the respective model; the ISM inflow comes from the top of the respective panel. Because 2D models are symmetric about the central (z-)axis, the horizontal axes can be either the x- or y-axis or any other axis within the xy-plane. Figure 1 shows contour plots in the logarithmic number density, log(n), measured in [cm⁻³]. The bottom-right panel is a linearly rescaled closeup of the top-right panel. Due to their differences in size and ISM inflow speed, the models have taken different times to reach stationarity. The model with u_ISM = 80 km s⁻¹ has reached stationarity of its upwind structure after 4t_pc ≈ 179 kyr, however, its downwind structure vacillates (discussed below). The two top-left panels depict the astrosphere in its most distended (left panel) and most compressed (centre panel) state at 19.96t_pc ≈ 895 kyr and 18.96t_pc ≈ 850 kyr, respectively. The model with u_ISM = 40 km s⁻¹ (top right panel) is depicted at 30t_pc ≈ 1.3 Myr, although its structure within the innermost 8 pc already reaches stationarity at 23t_pc ≈ 1.0 Myr. The model with u_ISM = 20 km s⁻¹ (bottom left panel) has reached stationarity at 300t_pc ≈ 13.4 Myr; its inner structure has already reached stationary after 70t_pc ≈ 3.1 Myr. The model with u_ISM = 10 km s⁻¹ (bottom centre panel) has not reached stationarity after 400t_pc ≈ 17.9 Myr, at which point the simulation was stopped, although the changes during the previous 150t_pc ≈ 7 Myr are minor (Baalmann 2021, Sect. 4.4.1). Compared to the age of λ Cephei, which Bouret et al. (2012) estimate to be 5 Myr, and the time since the star has been accelerated to its high relative speed 2.5 Myr ago according to Gvaramadze & Gualandris (2011), only the simulation times of the models with u_ISM ∈ {80, 40} km s⁻¹ are reasonable.

The homogeneous region of n = 11 cm⁻³ at the top of the panels is the ISM. The neighbouring region of higher densities is the outer astrosheath; the border between the two delineates the BS. The distances of the BS, the AP, and the TS along the central axis, z, are tabulated in Table 3. The model with u_ISM = 10 km s⁻¹ does not generate a BS because the ISM inflow is subfast. The strengths of the BSs also vary, which can be measured well with the compression ratio, s, the ratio of the number density behind the shock (i.e. in the outer astrosheath) to that in front of the shock (i.e. in the ISM). For the model with u_ISM = 80 km s⁻¹ the compression ratio is highest, s₈₀ = 3.64. The model with u_ISM = 40 km s⁻¹ has a notably weaker shock, s₄₀ = 2.73, and the shock of the model with u_ISM = 20 km s⁻¹ is weaker still, s₂₀ = 1.29. Another clear difference of the BS of the model with u_ISM = 20 km s⁻¹ as compared to the other two BSs is the opening angle. Whereas the BSs of the models with u_ISM ∈ {80, 40} km s⁻¹ are similar in shape, featuring an opening angle of roughly θ_{80,40} ≈ 80°, the BS opening angle of the model with u_ISM = 20 km s⁻¹ is almost θ₂₀ ≈ 150° (Baalmann 2021, Sect. 4.4.1).

All three superfast models feature a vacillating MD. This has been closely examined for the model with u_ISM = 80 km s⁻¹, where the MD vacillates between distances along the downwind z-axis of z_MD,80,min = −1.4 pc and z_MD,80,max = −3.4 pc (cf. Table 3). This periodic behaviour repeats every 2.70t_pc ≈ 121 kyr; the phase of distension takes 1.00t_pc ≈ 44.8 kyr and the phase of compression accordingly 1.70t_pc ≈ 79.2 kyr (Baalmann 2021, Sect. 4.4.1).

At the time of compression (top left panel of Fig. 1), there is a visible corner between the TS and the MD, the TP, which is located at (x, z) = (1.16 pc, −3.18 pc). From the TP, the TD extends downwind; it separates the astrotail with a number density of n ≈ 0.003 cm⁻³ from the inner astrosheath with n ≈ 0.005 cm⁻³. A reflected shock (RS) also begins at the TP at an angle of 30° to the vertical; across this TP, the number density jumps from n ≈ 0.004 cm⁻³ at the side closer to the AP to n ≈ 0.007 cm⁻³ at the side closer to the TD (Baalmann et al. 2021, Sect. 4.1). At the time of distension (top centre panel), no such corner between the TS and the MD is visible in the number density, and neither are the TD or the RS. In other parameters, such as the fast magnetosonic Mach number, M_f = u/υ_f, or the sonic Mach number, M_c = u/υ_c, both the TP and the TD are visible (Baalmann 2021, Sect. 4.4.1).

The MDs of the models with u_ISM ∈ {40, 20} km s⁻¹ vacillate as well, although their periodicity has not been examined more closely due to time constraints. However, for the model with u_ISM = 40 km s⁻¹ the TP and TD are always visible in the number density, even at the astrosphere’s most compressed state. The MD is located along the z-axis at z_MD,40,max = −4.37 pc at its most distended state after 33t_pc ≈ 1.48 Myr and at z_MD,40,min = −3.87 pc when most compressed after 37.8t_pc ≈ 1.69 Myr (cf. Table 3). A RS extends from the TP and is reflected multiple times by the AP and the TD, as can be seen in the bottom right panel of Fig. 1. For the model with u_ISM = 20 km s⁻¹, the astrosphere is most distended at 265t_pc ≈ 11.9 Myr with z_MD,20max = −5.6pc and most compressed at 288t_pc ≈ 12.9 Myr with z_MD,20,min = −4.5 pc (cf. Table 3). The model with u_ISM = 10 km s⁻¹ is in its most compressed state with z_MD,10,min = −4.5 pc for the first time at 282t_pc ≈ 12.6 Myr, and has not stopped its subsequent distension when the simulation was terminated at 400t_pc; it is therefore unknown whether its structure vacillates or not.

With the aim of probing the cause of the vacillating astrotails, a 2D model with u_ISM = 80 km s⁻¹ and the resolution of the respective 3D model was generated. The astrotail of this model did not vacillate, indicating that the stronger influence of numeric diffusion on more coarsely resolved grids inhibits the vacillating behaviour.

In Fig. 2, the stellar-centric distance of the discontinuities along the inflow axis of the 2D models as tabulated in Table 3 are plotted against the respective inflow speed. Because the model with u_ISM = 10 km s⁻¹ has not yet reached its most distended state, the distance of the distended MD is not plotted; there also is no BS for this model. In the dual logarithmic graph the distances of the TS, AP, and distended MD (MD_max) lie on straight lines for superfast inflows, namely, for u_ISM ∈ {20, 40, 80} km s⁻¹ with slopes of −0.96, −1.08, and −0.36, respectively. The respective distances of the subfast model, u_ISM = 10 km s⁻¹, do not lie on these lines. Similarly, no linear behaviour on the dual logarithmic graph can be observed for the distances of the BS and the compressed MD (MD_min). Unfortunately, the low number of available models does not allow to assign a high level of confidence to these results.

In Fig. 3, stellar-centric distances of the BS, AP, and MD at both its most compressed and most distended states are plotted against the distance of the TS (cf. Müller et al. 2006, Fig. 6). The ordinate axis for the distances of the BS and AP is scaled logarithmically, the ordinate axis for the two MD distances as well as the abscissa axis for the TS are scaled linearly. The distances of the BS appear to lie on a line on the semi-logarithmic axis system, as do the distances of the AP excluding the point for u_ISM = 80 km s⁻¹ (leftmost in the graph). The distances of the AP excluding the point for u_ISM = 10 km s⁻¹ (rightmost in the graph) also appear to lie on a line on the linear axis system, agreeing with Fig. 2. The distances of the MD appear concave on the linear axis system. This only partially agrees with the findings of Müller et al. (2006) based on a parameter study of 27 heliosphere-like models, who found linear correlations between the distances of the BS and AP to that of the TS.

Fig. 1 Slice planes of the 2D λ-Cephei-like models at different ISM inflow speeds. From left to right and top to bottom: uISM = 80 km s−1 in its most distended state, uISM = 80 km s−1 in its most compressed state, uISM = 40 km s−1, uISM = 20 km s−1, uISM = 10 km s−1; number density, log(n), in [cm−3]. Bottom right panel: close-up of the inner astrosheath for uISM = 40 km s−1; number density, n, in [10−3 cm−3]. Axes in [pc]. Taken with permission from Baalmann (2021, Fig. 4.16).

Fig. 2 Distances, ${\dot{M}}_{*}$, of the BS (black squares, solid line), AP (blue diamonds, dashed line), TS (purple circles, dotted line), and the MD in its most distended (MDmax, orange triangles, dash-dotted line) and compressed (MDmin, green pyramids, dash-dot-dotted line) positions along the central axis of the 2D models (Table 3) against the respective inflow speed, uISM. Both axes are scaled logarithmically.

Fig. 3 Distances, ${\left|B_{\text{SW}}\right|}_{r=r_{\text{SW}}}=B_{\text{SW}}$, of the BS, AP, and the MD in its most distended (max) and compressed (min) positions along the central axis of the 2D models (Table 3) against the respective TS distance, rTS, making use of both logarithmic (left ordinate axis, filled symbols) and linear (right ordinate axis, empty symbols) scaling; the abscissa has been scaled linearly. The AP has been plotted on both scales. Colours as in Fig. 2 with the exception of the linearly scaled AP (purple dotted line).

3.2 2D evolution

To facilitate the comprehension of the structures seen in the slice planes, the time evolution of the 2D model at u_ISM = 20 km s⁻¹ from its inception until reaching stationarity has been analysed. This is depicted for selected times in Fig. 4; the number density, n, is plotted linearly in order to better visualise the developing outer structure; the inner astrosphere cannot be resolved by the colour scale. At t = 0 (top left panel), the model shows the initial conditions: the ISM as a homogeneous domain of n = 11 cm⁻³ surrounding the radially decreasing gradient of the SW with n < 0.1 cm⁻³, which cannot be resolved by the linear colour scale. The levelling between these two domains is continuous and smooth. The BS and AP move outwards until t = 10t_pc ≈ 448 kyr (top centre-left panel); the MD and TS are unresolved. While these two inner shocks are still a single, perfectly spherical structure, the BS and AP are already distorted in the direction of the ISM inflow. After t = 10t_pc the BS continues moving outwards, whereas the AP is being pushed downwind by the ISM. This can be seen at the depicted slice plane for t = 40t_pc ≈ 1.79 Myr (top centre-right panel), where the BS has already moved out of the simulation grid in the downwind direction, the outer astrosheath has drastically grown in size, and plasma has piled up in front of the AP in the upwind direction. The ISM continues pushing the AP downwind but is balanced by the isotropic SW in the upwind direction. This causes a bulge on the AP close to the central axis, as can be seen at t = 50t_pc ≈ 2.24 Myr (top right panel); more material piles up in front of it. At t = 70t_pc ≈ 3.14 Myr (bottom left panel) the balance between the ISM and the SW has given the AP its bullet-like shape down to z ≈ −20 pc. A secondary BS moves from the AP outwards, whereas the original BS now moves inwards, cf. t = 150t_pc ≈ 6.72 Myr (bottom centre-left panel). At t = 200t_pc ≈ 8.97 Myr (bottom centre-right panel) the two BSs meet close to the central axis. The inner BS stays in position but the outer BS, farther away from the central axis, where it has not yet met the inner BS, moves further inwards. The AP has long since obtained its bullet-like shape within the entire simulation grid. At t = 300t_pc ≈ 13.4 Myr (bottom right panel), both BSs have aligned; cf. the bottom left panel of Fig. 1 for a logarithmic scaling at this time (Baalmann 2021, Sect. 4.4.2).

This time evolution is most emphatically not compliant with the physical evolution of an astrosphere but is presented to qualitatively illustrate the dynamical behaviour of the system and the timescales on which it occurs. It also serves as an aid to understand the origin of the visible structures, which will be used in Sect. 3.3. Similar behaviour can be observed for the models with u_ISM ∈ {80, 40} km s⁻¹, although their narrower outer astrosheaths drastically reduce the lifetime of the secondary BS.

Fig. 4 Slice planes of the 2D λ-Cephei-like models with an ISM inflow speed of uISM = 20 km s−1 at times t ∈ {0, 10, 40, 50, 70, 150, 200, 300} tpc (from left to right and top to bottom) in the number density, n, in [cm−3]; axes in [pc]. Taken with permission from Baalmann (2021, Fig. 4.17).

3.3 3D models

Cut slices of the corresponding 3D models are presented in Fig. 5, displaying the yx-plane (y-axis pointing rightwards, x-axis upwards) of the respective grid in the logarithmic number density, log(n), in the top row, and the linear number density, n, in the bottom row. As before, the linear colour scale cannot resolve the inner astrosphere (i.e. the TS and MD).

The 3D model with u_ISM = 80 km s⁻¹ (leftmost column), depicted at its time of stationarity, t = 5.84t_pc ≈ 262 kyr, agrees well with its 2D counterpart. Its MD does not vacillate due to the coarser angular resolution, which increases the effects of numeric diffusion and therefore inhibits instabilities. The positions of the discontinuities along the central axis, cf. Table 3, match those of the 2D model closely. Due to the comparatively coarse angular resolution, no RS is visible (Baalmann 2021, Sect. 4.4.3).

The 3D model with u_ISM = 40 km s⁻¹ (centre-left column) shows only minor differences with regard to its 2D counterpart as well. It is depicted at its time of stationarity, t = 35t_pc, which is similar to that of the 2D model. The comparatively diffuse AP and the lack of a RS are artifacts of the coarse angular resolution. The positions of the discontinuities along the central axis are given in Table 3 and agree fairly well with the 2D simulation, though it is notable that the MD lies further inward than its 2D counterpart even in its compressed state, and that the BS has a slightly larger opening angle (Baalmann 2021, Sect. 4.4.3).

The 3D model with u_ISM = 20 km s⁻¹ (centre-right column) is depicted at t = 160t_pc ≈ 7.17 Myr and not yet stationary. Its shows a similar double-BS structure than its 2D counterpart (Baalmann 2021, Sect. 4.4.3).

The 3D model with u_ISM = 10 km s⁻¹ (rightmost column), however, shows significant differences with regard to its 2D counterpart. It is depicted at 180t_pc ≈ 8.07 Myr, where the inner astrosphere is already stationarity but two bow waves still move outwards, akin to the two BSs of the previously examined models. The upwind structure of the 3D model lies considerably farther inwards compared to the 2D model, x_AP,10 = 6pc and x_TS,10 = 3.8 pc, compared to the 2D model’s z_AP,10 = 9.6 pc and z_TS,10 = 4.5 pc, whereas the MD is at a similar distance along the central axis, x_MD,10 = −5pc, compared to z_MD,10 = −4.8 pc (cf. Table 3). The entire astrotail is contained by the simulation grid, which is in stark contrast to the 2D analogue. An additional structure downwind of the astrotail is visible for x < −63.5 pc; the gradient of the number density significantly steepens its increase from n < 1cm⁻³ at x > −62 pc to n > 15 cm⁻³ at x < −68 pc. A structure of high density, n > n_ISM = 11 cm⁻³, lies farther downwind and continuously moves outwards. The material of this structure has been gathered by the second bow wave from the outer surface of the AP during the process of the AP’s tightening, cf. the time evolution during t ∈ [50, 150]t_pc of the 2D model with u_ISM = 20 km s⁻¹ (Sect. 3.2). This structure is unstable, both moving downwind and dispersing at the depicted time. Because this structure’s origin lies in the processes of the initial formation of the model astrosphere, which does not reflect the physical formation of a genuine astrosphere, the structure likely has no physical counterpart. Conceivably, a similar structure could be formed through different evolutionary phases of the SW and ISM. This structure does not exist in the 2D model; it is caused by the AP’s severe asymmetry before its compression due to the second bow wave, which in turn is an effect of the oblique magnetic field (Baalmann 2021, Sect. 4.4.3).

Due to the unphysical origin of this structure and the unrealistically high age of the astrosphere at its inception, a systematic study of this structure is unlikely to yield insightful results. However, stars with much longer lifetimes like M-type or even Sunlike stars may generate similar structures. It may be possible that such structures are common in astrospheric simulations but were not previously observed because most modelling approaches do not include the full astrotail (e.g. Katushkina et al. 2018; Meyer et al. 2020; Scherer et al. 2020; Herbst et al. 2020; Baalmann et al. 2020, but conversely Meyer et al. 2021).

Fig. 5 Slice of the yx-planes of the 3D λ-Cephei-like models at different ISM inflow speeds, uISM ∈ {80, 40, 20, 10} km s−1 (from left to right) in the number density, log(n) (top row) and n (bottom row); number density in [cm−3], axes in [pc]. The two left images in the bottom row feature number densities in the outer astrosheath above the values covered by the colour scale, up to n = 41.5 cm−3 and n = 30.6 cm−3, respectively. Taken with permission from Baalmann (2021, Fig. 4.18).

3.4 Synthetic observations

The different astrospheric models are projected in Fig. 6 in the Hα radiance as per Eq. (12), at a stellar distance of d₀ = 617 pc (Gaia Collaboration 2018). The four lines in the panels correspond to the four ISM inflow speeds, u_ISM ∈ {80, 40, 20, 10} km s⁻¹ (from top to bottom). The first column displays synthetic observations of the 2D models, which have been cyclindrically extruded to 3D models before line-of-sight-integrating; the geometry of the projection has been set up to view the slice planes of Fig. 1, namely, the yz-plane of the extruded model, face-on. The three right columns show synthetic observations of the respective 3D models, in the centre-left column face-on to the slice plane of Fig. 5, which is the yx-plane, namely, the ecliptic, of the model; in the centre-right column the models are projected edge-on to the slice plane, namely, face-on to the zx-plane, and in the right column tail-on to the astrosphere, namely, face-on to the zy-plane. All projections are subject to Moiré patterns, which are visible, for example, as lines along the coordinate axes.

For u_ISM = 80 km s⁻¹ (top row), the 2D projection and the face-on and edge-on projections of the 3D model agree well. This is as expected; because hydrodynamic effects are dominant over the magnetic field, as measured, for example, by the thermal and ram pressures versus the magnetic pressure, the 3D model must be similar to the 2D model and furthermore be reasonably symmetric about the central axis. Nevertheless, the opening angle of the 3D model’s BS is slightly larger due to the added magnetic pressure. The tail-on projection is rotationally symmetric about its centre; its radiance is highest within an outer ring with a radius of about 0.4°, where the lengths of the local lines of sight through the high-radiance outer astrosheath are longest (Baalmann 2021, Sect. 4.4.4).

A similar behaviour is apparent for u_ISM = 40 km s⁻¹ (second row). Again, the opening angle of the 3D model’s BS is slightly larger than its 2D counterpart. The ring of high radiance in the tail-on projection has a radius of about 1°, corresponding to the lines of sight with the longest components within the outer astrosheath (Baalmann 2021, Sect. 4.4.4).

The rotational symmetry is broken for u_ISM = 20 km s⁻¹ (third row), as the tail-on projection reveals. The ring of high radiance is no longer rotationally symmetric but elliptical. Its long axis is tilted by about 30° to the vertical axis. The face-on and edge-on projections are accordingly dissimilar; the opening angle of the edge-on view’s BS is larger than its face-on counterpart. The difference with regard to the 2D projection comes mostly from the difference in astrospheric evolution: the 2D model is already stationary, whereas the 3D model is still evolving, as can be seen, for example, by the double-BS structure that is apparent in both the face-on and the edge-on projections (Baalmann 2021, Sect. 4.4.4).

Some of the dissimilarities between the 2D and the 3D model with u_ISM = 10 km s⁻¹ (bottom row) are caused by the different times of projection as well. In the face-on and edge-on projections of the 3D model, the outbound bow waves are vaguely perceptible in the upwind direction, x > 0, as disks of slightly higher radiance with radii of about 4° and 8°, respectively; in the 2D projection, where the bow waves have already dispersed, the radiance within these regions is notably lower. The different shapes of the inner astrosphere, that is, those of the AP and the astrotail, are not caused by the difference in time, however. While the 2D structure is drop-shaped and extends out of the simulation grid in the downwind direction, the 3D structure varies in shape in the face-on and edge-on view and is considerably smaller. The high-density structure noted in Sect. 3.3 is visible as a high-radiance structure, appearing to be ring-like from the face-on perspective and irregular in the edge-on view. This asymmetry can also be seen in the tail-on projection, where the high-radiance structure is again tilted by 30° to the vertical axis. The irregular structure of the high-density object is, however, not readily apparent (Baalmann 2021, Sect. 4.4.4).

Fig. 6 Hα projections of λ-Cephei-like astrospheres with ISM inflow speeds (from top to bottom) uISM ∈ {80, 40, 20, 10} km s−1 onto the lb-plane; modelled in 2D (left), and in 3D projected face-on to the ecliptic (centre-left), edge on to the ecliptic (centre-right), and tail-on to the astrosphere (right). Hα radiance in [erg s−1 cm−2 sr−1], angular extent in degrees. Taken with permission from Baalmann (2021, Fig. 4.19).

4 Summary

Four models of λ-Cephei-like astrospheres with different inflow speeds of the ISM, u_ISM ∈ {80, 40, 20, 10}, were computed in both 2D and 3D. The 2D models were simulated without a SW magnetic field and an ISM magnetic field parallel to the ISM inflow, whereas the 3D models were simulated with a simplified Parker spiral for the SW and an oblique field for the ISM. The astrospheric evolution was analysed based on the example of the 2D model with u_ISM = 20 km s⁻¹. Synthetic observations in Hα were generated for all models.

It has been found that, due to the dominance of hydrodynamic effects over the magnetic field, the 2D and 3D models are in good agreement for u_ISM ∈ {80, 40} km s⁻¹, and that the 3D model is reasonably symmetric about its central axis; this is also apparent in the projections. However, the 2D models generated a vacillating MD and astrotail, which is inhibited by the greater influence of numeric diffusion in the more coarsely refined 3D models.

For the models with u_ISM = 20 km s⁻¹, where the ISM inflow is barely superfast and supersonic, this symmetry is broken; the 3D model astrosphere is distorted by the oblique ISM magnetic field. However, with the exception of this distortion, the 2D and 3D models still agree remarkably well. At this ISM inflow speed the evolution of the model contains a sizable double-BS structure for a significant time, which is apparent in the synthetic observations. The physical time until the simulation has reached stationarity typically exceeds the lifetime of the star.

The 2D and 3D models with u_ISM = 10 km s⁻¹, where the ISM inflow is subfast and subsonic and therefore does not generate a BS, are not in good agreement. The extent of the astrotail is much smaller in the 3D model, and a high-density structure develops at its downwind end. This structure is highly asymmetric and irregular; it dominates the synthetic Hα observations. The origin of this structure lies in the oblique magnetic field; it cannot be ruled out that such a structure is common at the astrotails of model astrospheres either with or without a BS.

The following conclusions can be summarised:

1. 2D MHD simulations can be sufficient to produce accurate models of the large-scale astrospheric structure if hydrodynamic effects are dominant over the ISM magnetic field (Sect. 3.3). However, no SW magnetic field can be simulated with this approach;

2. Fine-resolution models, which are less computationally expensive in 2D compared to 3D, feature a vacillating astrotail (Sect. 3.1). This is inhibited in coarse-resolution models due to the increased impact of numeric diffusion (Sect. 3.3);

3. Astrospheres around high-mass stars with comparably slow ISM inflows require extensive simulation times to reach stationarity, potentially much longer than the actual lifetimes of the modelled stars (Sect. 3.2). This implies that observed astrospheres of such short-lived stars must be considered to be in a transient state.

4. The model with an ISM inflow too slow to generate a BS shows significant differences between its 2D and 3D simulations (Sect. 3.3);

5. Within this model, the oblique magnetic field has generated a hitherto unobserved high-density structure at the downwind end of the astrotail, which dominates the synthetic Hα radiance map (Sect. 3.4).

To conclude, the use of 2D MHD can be expedient with regard to its more computationally expensive 3D counterpart as long as the importance of the stellar magnetic field and of the orientation of the ISM magnetic field is negligible. While the 2D approach predetermines cylindrical symmetry, additional features that stem from the finer resolution, such as a vacillating astrotail, can emerge at no additional computational cost with regard to the 3D approach. When modelling in 3D, it can prove vital to extend the simulation boundary to include the full astrotail in order to reproduce all important features of the model.

Acknowledgements

This work makes use of the ColorBrewer colour scales, designed by Cynthia A. Brewer, Geography, Pennsylvania State University (www.ColorBrewer.org). KS is grateful to the Deutsche Forschungsgemeinschaft (DFG), funding the project SCHE334/9-2. JK acknowledges financial support through the Ruhr Astroparticle and Plasma Physics (RAPP) Center, funded as MERCUR project St-2014-040.

References

All Tables

All Figures

Fig. 1 Slice planes of the 2D λ-Cephei-like models at different ISM inflow speeds. From left to right and top to bottom: uISM = 80 km s−1 in its most distended state, uISM = 80 km s−1 in its most compressed state, uISM = 40 km s−1, uISM = 20 km s−1, uISM = 10 km s−1; number density, log(n), in [cm−3]. Bottom right panel: close-up of the inner astrosheath for uISM = 40 km s−1; number density, n, in [10−3 cm−3]. Axes in [pc]. Taken with permission from Baalmann (2021, Fig. 4.16).

In the text

Fig. 2 Distances, ${\dot{M}}_{*}$, of the BS (black squares, solid line), AP (blue diamonds, dashed line), TS (purple circles, dotted line), and the MD in its most distended (MDmax, orange triangles, dash-dotted line) and compressed (MDmin, green pyramids, dash-dot-dotted line) positions along the central axis of the 2D models (Table 3) against the respective inflow speed, uISM. Both axes are scaled logarithmically.

In the text

Fig. 3 Distances, ${\left|B_{\text{SW}}\right|}_{r=r_{\text{SW}}}=B_{\text{SW}}$, of the BS, AP, and the MD in its most distended (max) and compressed (min) positions along the central axis of the 2D models (Table 3) against the respective TS distance, rTS, making use of both logarithmic (left ordinate axis, filled symbols) and linear (right ordinate axis, empty symbols) scaling; the abscissa has been scaled linearly. The AP has been plotted on both scales. Colours as in Fig. 2 with the exception of the linearly scaled AP (purple dotted line).

In the text

	Fig. 4 Slice planes of the 2D λ-Cephei-like models with an ISM inflow speed of uISM = 20 km s−1 at times t ∈ {0, 10, 40, 50, 70, 150, 200, 300} tpc (from left to right and top to bottom) in the number density, n, in [cm−3]; axes in [pc]. Taken with permission from Baalmann (2021, Fig. 4.17).
In the text

Fig. 5 Slice of the yx-planes of the 3D λ-Cephei-like models at different ISM inflow speeds, uISM ∈ {80, 40, 20, 10} km s−1 (from left to right) in the number density, log(n) (top row) and n (bottom row); number density in [cm−3], axes in [pc]. The two left images in the bottom row feature number densities in the outer astrosheath above the values covered by the colour scale, up to n = 41.5 cm−3 and n = 30.6 cm−3, respectively. Taken with permission from Baalmann (2021, Fig. 4.18).

In the text

Fig. 6 Hα projections of λ-Cephei-like astrospheres with ISM inflow speeds (from top to bottom) uISM ∈ {80, 40, 20, 10} km s−1 onto the lb-plane; modelled in 2D (left), and in 3D projected face-on to the ecliptic (centre-left), edge on to the ecliptic (centre-right), and tail-on to the astrosphere (right). Hα radiance in [erg s−1 cm−2 sr−1], angular extent in degrees. Taken with permission from Baalmann (2021, Fig. 4.19).

In the text