© ESO 2020

1 Introduction

The interactions of interstellar and stellar plasmas result in distinct structures around the host stars, as first analytically modelled by Weaver et al. (1977). They have been simulated in one-dimensional (1D) hydrodynamics (HD; e.g. Arthur 2012), 2D HD (e.g. Scherer et al. 2016a; del Valle & Pohl 2018, and references therein), 2D magnetohydrodynamics (MHD; e.g. Meyer et al. 2017), and more recently also in 2.5D MHD (e.g. van Marle et al. 2014), 3D HD (e.g. Mohamed et al. 2012), and 3D MHD (e.g. Gvaramadze et al. 2018; Scherer et al. 2019). The resulting models have been successfully compared to observational images and data in various observables (e.g. Kobulnicky et al. 2017; Sánchez-Ayaso et al. 2018, and references therein). The shapes of these structures strongly depend on the physical parameters of the star and on the surrounding interstellar medium (ISM; see e.g. Gvaramadze et al. 2012), as well as on the orientation of the stellar motion relative to the observer’s line of sight (LOS; see e.g. Tarango-Yong & Henney 2018). Because Earth is located inside the heliosphere, numerical models of the latter can be compared to observational data with much higher spatial resolution than for other astrospheres. Recently, both numerical 3D MHD models and observational studies have been used to explore various aspects of the heliosphere in detail (see e.g. Opher et al. 2016; Dialynas et al. 2017; Kornbleuth et al. 2018; Guo et al. 2019). Because observational images are 2D projections of a 3D structure, 3D models must be treated similarly to allow comparisons between simulation and observation. Such a 3D model can be rotated at will before projecting, simulating different orientations of the modelled structure to the observer’s LOS.

We here use 3D single-fluid MHD simulations of astrospheres as in Scherer et al. (2016b) as a basis for projections in various observables. The available physical parameters of the simulations are used to calculate further physical values in each model cell that can then be projected through a modified LOS integration to yield observables such as the radiation fluxes of bremsstrahlung, cyclotron radiation and Hα emission, the column density, or the rotation measure. The simulation volume can be rotated about all three axes and shifted to a chosen position in Galactic coordinates before integration.

The astrospheres of two notable stars, λ Cephei and the Sun, have been examined for example by Scherer et al. (2016b). The Sun is interesting because it is our host star; Earth is embedded in the solar astrosphere, the heliosphere (Zank 1999; Opher 2016; Pogorelov et al. 2017). λ Cephei is the brightest runaway O star in the sky (Scherer et al. 2016a), and therefore is the object most likely to have the most visible astrosphere. The astrosphere of λ Cephei has been proposed as a sink for cosmic rays (Scherer et al. 2015), which has been examined further by Scherer et al. (2016c).

The shock structures of the models are assumed to be comparable to the shock structure described by Scherer et al. (2016b): a termination shock (TS) close to the star, where the stellar wind jumps from supersonic to subsonic speeds, an astropause (AP) separating the stellar and interstellar fluids, and a bow shock (BS) in the direction of the incoming ISM in front of the star, where the interstellar fluid jumps from supersonic to subsonic speeds. In the downwind direction, that is, the direction opposite to the interstellar inflow, the stellar wind jumps from supersonic to subsonic speeds at the Mach disc (MD). A BS may not exist in all models, depending on the relative speed between the star and the surrounding ISM (see e.g. McComas et al. 2012; Zank et al. 2013). The region between the TS and the AP is called the inner astrosheath, and the region between the AP and the BS is the outer astrosheath. The outer astrosheath is divided intoa hot outer astrosheath (HOA) close to the BS and a cold outer astrosheath (COA) close to the AP.

Details of the simulations, including the parameters of the stars, are presented in Sect. 2. An overview of the available physical parameters is given in Sect. 3, and the projection methods are outlined in Sect. 4. The results are presented in Sect. 5 and compared with observational data in Sect. 6. Conclusions are drawn in Sect. 7.

2 Astrosphere models

The astrospheres were modelled in 3D single-fluid MHD with the CRONOS code (Kissmann et al. 2018), a semi-discrete finite-volume code that advances Eqs. (1)–(4) in time, opting for an HLL Riemann solver and a second-order Runge-Kutta scheme, $$\begin{array}{l}\frac{\partial n}{\partial t}+\nabla \cdot \left(nu\right)=0,\\ \frac{\partial }{\partial t}\left(mnu\right)+\nabla \cdot \left(mnu\otimes u\right)+\nabla p+\frac{1}{{\mu }_0}B\times \left(\nabla \times B\right)=0,\\ \frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+p+\frac{1}{2{\mu }_0}{\left|B\right|}^2\right)u-\frac{1}{{\mu }_0}\left(u\cdot B\right)B\right]=0,\\ \frac{\partial B}{\partial t}=-\nabla \times E=\nabla \times \left(u\times B\right)\text{},\end{array}$$

where n, u, p, B, e, E are the number density, fluid velocity, thermal pressure, magnetic induction, total energy density, and electric field. m and μ₀ are the proton mass and the vacuum permeability, respectively; ⊗ represents the dyadic product. The system is closed by the equations $$\begin{array}{l}e=\frac{p}{\gamma -1}+\frac{1}{2}mn{\left|u\right|}^2+\frac{1}{2{\mu }_0}{\left|B\right|}^2,\\ \nabla \cdot B=0,\end{array}$$

where γ is the polytropic index. Instead of the thermal pressure, we usually considered the temperature T = 3p∕(4nk_B), where k_B is the Boltzmann constant. Factor 1∕2 stems from the assumption of thermodynamic equilibrium between electrons and protons (p = p_e + p_p = 2p_p).

Additional terms for heating follow Reynolds et al. (1999) and were treated as in Kosiński & Hanasz (2006), and cooling following Schure et al. (2009) was added at the right-hand side of Eq. (3), allowing photoelectric heating by dust, dissipation of interstellar turbulence, and Coulomb collisions with cosmic rays, and in particular, radiative cooling in the temperature range T ∈ [10⁴, 10⁷] K. The polytropic index was set to γ = 1.66667 ≈ 5∕3.

The model was set up as the interior of a near-spherical polyhedron with cells in spherical coordinates (r, ϑ, φ) that were equidistant in the radial direction e_r and equiangular in directions e_ϑ,φ perpendicular to e_r. The entire angular ranges of ϑ ∈ [−90°, 90°] and φ ∈ [−180°, 180°] are covered, while the radial range extends over the interval r ∈ [r_SW, r_ISM], r_SW ≠ 0. The increments and therefore the radial and angular resolutions depend on the model.

Boundary values for the stellar wind (SW) and the ISM were set for the inner- and outermost cells (at r_SW,ISM): both v and B of the ISM are homogeneous with angles ϑ_v,B, φ_v,B, while v of the SW is spherically symmetric with only a radial component. B of the SW was set up as a spiral field following Parker (1958). Likewise, n and T were set and kept fixed at r_SW,ISM for the SW and ISM.

By adjusting the boundary values to observed or derived parameters of stars and their surroundings, the astrospheres of these stars were simulated. While any star may be examined this way, the two objects of particular interest, the Sun and λ Cephei (see Sect. 1), have been closely investigated. We briefly examine the astropheres of two other stars, those of Proxima Centauri and V374 Pegasi. The stellar winds of these stars have been modelled by Garraffo et al. (2016) and Vidotto et al. (2011), respectively. The boundary values of the corresponding models are given in Table 2; for λ Cephei, see Scherer et al. (2016a, 2019), and for the Sun, see Pogorelov et al. (2017). The values for the ISM of Proxima Centauri were adapted from the Sun; the stellar magnetic field was taken from Reiners & Basri (2008), and the density and wind speed were selected to produce a large shock structure. For V374 Pegasi, the ISM was assumed to be a warm ionized medium (WIM, see e.g. McKee & Ostriker 1977) with corresponding parameters, while values for the SW were taken from Vidotto et al. (2011). The numbers of simulated cells also differ between the models (see Table 1). A more detailed description of the simulation can be found in Scherer et al. (2016a).

3 Physical parameters

Of the dynamic variables, only the number density is directly related to an observable: to the column density. To investigate the projections in other observables, their unprojected physical values must be calculated first. In the cases of the total fluxes of bremsstrahlung (Sect. 3.1) and synchrotron or cyclotronradiation (Sect. 3.2), these values are the total energy loss rates, while the Hα flux (Sect. 3.3) corresponds to the number of recombinations per time and volume, that is, the recombination rate. The rotation measure (Sect. 3.4) requires a different approach (see Sect. 4).

3.1 Bremsstrahlung

When we assume that the electron speed follows a Maxwellian distribution, the total energy loss rate for thermal bremsstrahlung is $$-{\left(\frac{\text{d}E}{\text{d}t}\right)}_{\text{brs}}\left[\frac{\text{erg}}{{\text{cm}}^3\text{s}}\right]=1.435\times {10}^{-27}{Z}^2\sqrt{T}{g}_{\text{ff}}(T)n{n}_{\text{e}},$$(7)

with Z the mass number (here Z = 1), T the temperature in K, n, n_e the number densities of the ions and electrons in cm⁻³ (here n = n_e), and g_ff(T) the temperature-dependent Gaunt factor for free–free transitions (Longair 1992). A common approximation for g_ff (T) is its temperature-averaged value of ḡ = 1.2. Gaunt factors for any desired non-relativistic temperature were calculated here following van Hoof et al. (2014) and result in values g_ff (T) ∈ [1.0, 1.5) for the temperature range of the examined astrosphere models.

The emission spectrum of bremsstrahlung was assumed to be flat up to a critical frequency $${\nu }_{\max }=\frac{4\pi m_{\text{e}}{v}^2}{h},$$(8)

with m_e the electron mass, v the (thermal) electron speed, and h the Planck constant (Longair 1992). Because the total spectrum of the complete model is the weighted sum of the spectra of all cells, the total spectrum will be flat up to the lowest critical frequency. This is useful for calculating the spectral flux F_ν by dividing the total flux F by the frequency range, that is, the critical frequency. This is only a preliminary result, however, because self-absorption and similar processes will heavily modify a flat emission spectrum. These processes are currently not considered.

3.2 Synchrotron and cyclotron radiation

The total energy loss rate of a single electron by synchrotron radiation is $$-{\left(\frac{\text{d}E}{\text{d}t}\right)}_{\text{syn}}={\sigma }_{\text{T}}c\frac{B^2}{4\pi }{\gamma }^2{\left(\frac{v}{c}\right)}^2{\sin }^2\alpha ,$$(9)

with σ_T the Thomson cross-section, γ the Lorentz factor, c the speed of light, B, v the absolute values of the magnetic induction and the electron speed, and α the pitch angle between v and B (Blumenthal & Gould 1970). With γ ≈ 1, Eq. (9) becomes the total energy loss rate of cyclotron radiation for non-relativistic electrons. The electron speed contains both the absolute value of the fluid velocity and the thermal speed calculated through $$v_{\text{therm}}=\sqrt{\frac{8k_{\text{B}}T}{\pi m_{\text{e}}}}$$(10)

(Condon & Odishaw 1958), so that the angle between the electron motion and the magnetic field is no longer equal to α. The factor sin²α must be substituted by its average value of 2∕3 in Eq. (9).

The emitted power of the ultra-relativistic synchrotron spectrum is of the form $$P_{\text{em}}(\nu )=\frac{\sqrt{3}{q}_{\text{e}}^{3}B}{m_{\text{e}}{c}^2}\frac{\nu }{{\nu }_{\text{c}}}{\displaystyle {\int }_{\nu /{\nu }_{\text{c}}}^{\infty }{K}_{5/3}}(s)\text{d}s,$$(11)

with ν the frequency, q_e the elementary charge, ν_c the critical frequency, and K_5∕3 the modified Bessel function of the second kind (Blumenthal & Gould 1970). The critical frequency is $${\nu }_{\text{c}}=\frac{3}{2}{\nu }_{\text{g}}{\gamma }^2\sin \alpha ,$$(12)

with ν_g = eB∕(2πm_e) the gyration frequency. Again, the factor sinα must be substituted by π∕4 because α is no longer the angle between v and B. Integrating over the entire spectral range ν ∈ [0, ∞) yields Eq. (9). If the electron speed is non-relativistic, the total power is emitted at ν_g. At higher speeds, radiation is emitted at the gyrofrequency harmonics ${\nu }_l=l{\nu }_{\text{g}},l\in \mathbb{N}$; the power decreases by a factor of ${(v/c)}^2$ for every harmonic. The harmonics were ignored in the calculations because the electron speed remains lower than 0.1c in all examined astrospheres. By multiplying the total power with the number density, we calculated the power density.

3.3 Hα emission

Following Reynolds (1984), the photon flux of Hα is connected to the recombination rate by $$I_{\text{H}\alpha }\left[\frac{\text{photons}}{{\text{cm}}^2\text{s}\text{sr}}\right]=\frac{\epsilon {\tilde{r}}_{\text{rr}}}{4\pi },$$(13)

with ε the average number of Hα photons produced per recombination, and ${\tilde{r}}_{\text{rr}}$ the recombination rate along the LOS in cm⁻². The radiation flux was calculated by multiplying with the Hα photon energy hν_Hα, where ν_{H α} = 456.81 THz is the Hα frequency. Because ${\tilde{r}}_{\text{rr}}$ is the LOS-integrated case B recombination rate, $${\tilde{r}}_{\text{rr}}={\displaystyle {\int }_s{r}_{\text{B}}}\text{d}s,$$(14)

the emitted power density of a selected volume is $$P_{\text{H}\alpha }=h{\nu }_{\text{H}\alpha }\epsilon r_{\text{B}}.$$(15)

The case B recombination rate r_B describes recombinations to all atomic states n except for the ground state n = 1; recombination to the latter would result in subsequent ionization and therefore in no net recombination. The recombination rate can be calculated by the electron and ion number densities $n_{\text{e},{\text{H}}^{\text{+}}}$ and the case B recombination rate coefficient α_B, which is equivalent to the sum over all recombination rate coefficients α_n except for that of the ground state: $$r_{\text{B}}=n_{\text{e}}{n}_{{\text{H}}^{\text{+}}}{\alpha }_{\text{B}}=n_{\text{e}}{n}_{{\text{H}}^{\text{+}}}{\displaystyle \sum _{n=2}^{\infty }{\alpha }_n}.$$(16)

The temperature-dependent coefficients α_n(T) were taken from Mao & Kaastra (2016) and are valid for the entire observed temperature range T ∈[10¹, 10⁸] K for states n ∈ [1, 16]. The coefficients were not only resolved by the principal quantum number n, but also by the azimuthal quantum number l; the α_n were calculated by summation over all applicable l.

Because no estimation for ε is given by Reynolds (1984), a different approach was chosen to calculate εα_B. The emission coefficient for the atomic transition i → k is given by Draine (2011) as $$j_{i\to k}=\frac{n_e{n}_{{\text{H}}^{+}}}{4\pi }\frac{A_{i\to k}h{\nu }_{i\to k}}{{\sum }_l{A}_{i\to l}}\cdot \left[{\alpha }_i+{\displaystyle \sum _{l,l>i}{\alpha }_l}{P}_{\text{B},l,i}\right]\text{},$$(17)

where A_i→k is the transition probability for spontaneous emission of transition i → k, taken from Wiese & Fuhr (2009), and P_B,i,k is the branching ratio, that is, the probability that an atom in state i decays through state k, $$\begin{array}{l}{P}_{\text{B},i,k}=\left(A_{i\to k}+{\displaystyle \sum _{j<n<i}{A}_{i\to n}}{P}_{\text{B},n,k}\right)/{\displaystyle \sum _{n<i}{A}_{i\to n}}\\ P_{\text{B},i,i+1}=A_{i\to i+1}/{\displaystyle \sum _{n<i}{A}_{i\to n}}\text{}.\end{array}$$

An effective recombination rate coefficient α_eff,i→k was introduced to shorten Eq. (17) to $$j_{i\to k}=\frac{n_e{n}_{{\text{H}}^{+}}}{4\pi }{\alpha }_{\text{eff},i\to k}h{\nu }_{i\to k}.$$(20)

It directly ensues from the above equations that α_eff,Hα = εα_B. Because Hα corresponds to the transition 3 → 2 and the α_n quickly decreases for higher n, it follows that $${\alpha }_{\text{eff},3\to 2}\approx {\alpha }_3+{\sum }_{n=4}^{16}{\alpha }_n{P}_{\text{B},n,3}.$$(21)

Comparison of calculated values of α_eff,Hα and α_B yields $\epsilon \approx {0.2643}_{-0.019}^{+0.022}$ in the temperature range T ∈ [10¹, 10⁸] K.

3.4 Faraday rotation

As a measure for the Faraday rotation, the Faraday depth was used (Van Eck et al. 2017): $$\varphi (d)\left[\frac{\text{rad}}{{\text{cm}}^2}\right]=0.812{\displaystyle {\int }_d^0\left(\frac{n_{\text{e}}}{{\text{cm}}^{-3}}\right)}\left(\frac{B}{\mu \text{G}}\right)\cdot \left(\frac{\text{d}l}{\text{pc}}\right),$$(22)

with d the distance from the source. If d is equal to the distance between observer and source, the Faraday depth is equal to the rotation measure. Instead of using the magnetic induction and the line element as vectors, B may be substituted by its component parallel to the LOS B_|| and dl by its scalar counterpart dl.

4 Projection methods

The model projections were performed in Galactic coordinates (d, b, l), with the modulus d, the Galactic latitude b, and the Galactic longitude l; b and l cover the ranges b ∈ [−90°, 90°], l ∈ [−180°, 180°], with b = 0 being the Galactic plane (i.e. the xy-plane, z = 0). The model cells (i, j, k) are denominated by the Cartesian coordinates to their centres O(i, j, k); henceforth (κ) = (i, j, k) is used as a shorthand. The model can be rotated by angles α, β, γ around the x-, y-, and z-axes via rotation matrices A(α), B (β), C (γ) and shifted to a new position of the centre of the computational domain D(d₀, b₀, l₀), see Fig. 1 (upper row). The new position vectors to the model cells are $$N(\kappa )=A(\alpha )\cdot B(\beta )\cdot C(\gamma )\cdot O(\kappa )+D(d_0,b_0,l_0),$$(23)

with the rotation matrices $$\begin{array}{lll}A(\alpha )\hfill & =\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & \sin \alpha \\ 0& -\sin \alpha & \cos \alpha \end{array}\right),\hfill & \hfill \\ B(\beta )\hfill & =\left(\begin{array}{ccc}\cos \beta & 0& -\sin \beta \\ 0& 1& 0\\ \sin \beta & 0& \cos \beta \end{array}\right),\hfill & \hfill \\ C(\gamma )\hfill & =\left(\begin{array}{ccc}\cos \gamma & \sin \gamma & 0\\ -\sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right).\hfill & \hfill \end{array}$$

The shifting vector D(d₀, b₀, l₀) arises as theresult of transforming from Galactic to Cartesian coordinates: $$D(d_0,b_0,l_0)=\left(\begin{array}{l}{d}_0\cos b_0\cos l_0\hfill \\ d_0\cos b_0\sin l_0\hfill \\ d_0\sin b_0\hfill \end{array}\right).$$(27)

The scalar values n(κ), T(κ) of the simulation are given for each cell (κ) and need no further modification, whereas the vectorial values must be rotated. Both the fluid velocity and the magnetic induction are available in spherical coordinates v_s(κ), B_s (κ), as is a transformation matrix T (κ) from spherical to Cartesian coordinates. The new vectors ξ ∈{v, B} follow $$\xi (\kappa )=A(\alpha )\cdot B(\beta )\cdot C(\gamma )\cdot T(\kappa )\cdot {\xi }_{\text{s}}(\kappa ).$$(28)

The Galactic coordinates (d(κ), b(κ), l(κ)) of the new position of the cell centres N(κ) were calculated through the usual transformation $$\begin{array}{l}d(i,j,k)=\sqrt{N_x^2(\kappa )+N_y^2(\kappa )+N_z^2(\kappa )},\\ b(i,j,k)=\arcsin \left(\frac{N_z(\kappa )}{d(\kappa )}\right),\\ l(i,j,k)=\arctan \left(\frac{N_y(\kappa )}{N_x(\kappa )}\right).\end{array}$$

Because each model is spherical in nature, its radius R is half the difference of the largest and smallest Cartesian x-coordinate, $$R=\frac{1}{2}\left(\max (N_x(\kappa ))-\min (N_x(\kappa ))\right).$$(32)

The model was then projected onto the surface of a sphere around the origin, where the observer is located. The exact coordinates (b, l) of the projected point of a cell centre on the 2D surface is identical to the angular items of the 3D coordinates of the cell centre. However, the surface of the sphere was divided into a grid that set the resolution of the projection (see below).

If the distance of the computational domain centre to the origin d₀ is shorter than the model radius R, that is, if the origin is inside the model, the model projection will occupy the entire sky (the sphere of projection). If d₀ > R, the spherical model is projected as a disc with an angular diameter ⊘ < 180°. The angular radius Δb follows fromsimple geometry as $$\Delta b=\frac{1}{2}\text{⊘}=\arcsin \left(\frac{R}{d_0}\right);$$(33)

the symbol Δb was chosen because the angular radius is the angular interval of latitude between the disc centre and its edge in every geometric configuration.

If the latitude of the computational domain centre b₀ (henceforth called the elevation) is high or low enough for the projection disc to extend to or over the poles of the projection sphere (|b₀ |+ Δb ≥ 90°), the disc is spread out over all longitudes l ∈ [−180°, 180°]. For an elevation b₀ ≠ 0° and a smaller angular radius |b₀ | + Δb < 90°, the maximum interval in the longitude non-trivially depends on $\left\{b_0,\Delta b\right\}$ and is approximated by $$\Delta l=\arctan \left(\frac{\tan \left(\Delta b\right)}{\cos \left(b_0+\frac{b_0}{{90}^{°}}\Delta b\right)}\right).$$(34)

Because the position of the astrosphere model in the sky (or rather on the projected sphere) has no bearing for most purposes, b₀ = 0°, l₀ = 0° was used in most cases. In this case, the angular extent in both angular coordinates is identical, Δb = Δl. Even if the actual position of a star was used as (b₀, l₀), most stars appear close enough to the Galactic plane (b₀ ≈ 0°) that Δl ≈ Δb.

Because the outer edges of the 2D grid are now known as b₀ ± Δb, l₀ ± Δl, its axes were divided into dim(b, l) intervals. Here dim (b, l) is the number of model cells in ϑ, φ with additional resolution factors k_b,l. If the model is centred at the origin (d₀ = 0°) and k_b,l = 1, the 2D grid cells correspond to the columns of 3D model cells in the e_r -direction; the centres of all cells inside a column line up exactly. If the model is centred outside the origin (d₀ > R), resolution factors k_l = 2k_b must be used to obtain square grid cells because the models usually have twice as many cells in φ as in ϑ to obtain the same angular resolution.

The integration process was rather simple: when the angular coordinates (b, l) of the centre of a model cell appear inside one of the grid cells, its physical parameter H(κ) (e.g. the total energy loss rate of bremsstrahlung or the number density) were added to the grid cell value. Depending on the physical parameter H(κ), two types of weighting lead to two types of integration: LOS integration weights the cell by its length L(κ) along the LOS, whereas volumetric integration weights the cell by its volume V (κ) (see Fig. 1 bottom).

For most physical parameters, volumetric integration is the appropriate method. For all radiation types, the particles in the entire volume contribute to the total emission. Here, the parameter H(κ) is a measure of emission per volume, so the product of H(κ) and volume V (κ) is the total amount of emission. If H(κ) is the number density, H(κ)V (κ) gives the number of particles in the entire integration volume. Summing ∑ H(κ)V (κ) and dividing by the entire volume (∑ V (κ)) gives the volume-averaged number density. For the radiation parameters (bremsstrahlung, cyclotron, and Hα), the product H(κ)V (κ) divided by 4πd²(κ) gives the radiation flux of the cell, or more formalized, $${\text{Res}}_{\text{vol}}={\displaystyle \sum _{i=1}^N\frac{H_i{V}_i}{4\pi d_i^2}},$$(35)

with Res_vol the volumetric integration result, N cells inside the integrated volume, H_i, d_i, V_i the value of the physical parameter, the distance to the origin, and the volume of cell i ∈ [1, N]. The volume of cell (κ) is calculated through $$\begin{array}{ll}V(\kappa )=\hfill & {\displaystyle {\int }_{r_{\text{in}}(\kappa )}^{r_{\text{out}}(\kappa )}\text{d}}r{\displaystyle {\int }_{\phi (\kappa )-\frac{1}{2}\delta \phi }^{\phi (\kappa )+\frac{1}{2}\delta \phi }\text{d}}\phi {\displaystyle {\int }_{\vartheta -\frac{1}{2}\delta \vartheta }^{\vartheta +\frac{1}{2}\delta \vartheta }\text{d}}\vartheta r^2\cos \vartheta \hfill \\ =\hfill & \frac{1}{3}\left(r_{\text{out}}^3-r_{\text{in}}^3\right)\cdot \delta \phi \cdot \left(\sin \left(\vartheta +\frac{\delta \vartheta }{2}\right)-\sin \left(\vartheta -\frac{\delta \vartheta }{2}\right)\right),\hfill \end{array}$$(36)

where r(κ), φ(κ), ϑ(κ) are the spherical coordinates of cell (κ) in the simulation-centric coordinate system with ϑ = 0 the equatorial plane, δr, δφ, δϑ are the radial and angular increments of the simulation cells, and the shorthand notation r_out,in(κ) = r(κ) ± δr∕2.

For some physical parameters, volumetric integration is not the correct approach. If the observable depends on the length of the LOS instead of the emission volume, as is the case for the column density or the rotation measure, the cells must be weighted by their length L along the LOS instead of their volume V. Because L strongly depends on the position and orientation of the cell, an approximation was made to reduce the calculation runtime. Unlike the volumetric integration, where a model cell is treated as a finite volume with uniform physical values such as density or emissivity, a true LOS integration implies an infinitesimal extent in the directions perpendicular to the LOS. For a finite volume there is an infinite number of line segments parallel to the LOS through the volume in question; the length of these line segments along the LOS depends on the problem geometry, but varies widely in general.

In the special case when the centre of a model coincides with the origin (d₀ = 0), all LOS are equal to the radial rows of cells. The length of each cell along the LOS is the radial increment δr of the cell, which is equal for all cells. For other geometric arrangements, this is only true for the LOS that goes through the model centre. All other LOS will go through cells at an angle α(κ), calculated through $$\alpha (\kappa )=\sqrt{{\left(l(\kappa )-l(n)+\frac{1}{2}\delta l\right)}^2+{\left(b(\kappa )-b(m)+\frac{1}{2}\delta b\right)}^2},$$(37)

with (b(m), l(n)) the coordinate of the grid cell (m, n). The corresponding total angular increment is $$\delta \alpha (\kappa )=\min \left\{\left|\frac{\delta b}{\cos \left(x(\kappa )\right)}\right|,\left|\frac{\delta l}{\sin \left(x(\kappa )\right)}\right|\right\},$$(38)

with the shorthand $x=\arctan \left(l(\kappa )/b(\kappa )\right)$, so $$\delta \alpha (\kappa )=\min \left\{\delta b\sqrt{{\left(\frac{l(\kappa )}{b(\kappa )}\right)}^2+1},\delta l\frac{\sqrt{l{(\kappa )}^2+b{(\kappa )}^2}}{l(\kappa )}\right\}.$$(39)

To simplify the problem of a LOS going through a model cell, the latter was approximated by a conical frustum of height δr between base and top, an opening angle of δα∕2, and a central point r (so that the base and top are located at r ± δr along an r-axis). The exact model cell takes the shape of a truncated pyramid with a square base and top, where the length of the edges of base and top are equal to the diameters of the base and top of the approximated model cell. The angle α described above is now the angle between the central axis of the cone and the LOS direction. Because the integration method only selects cells with a central point inside the LOS, the length of the LOS through the conical frustrum can be selected as the length of the line segment parallel to the LOS through the cell centre. Thus, the problem is reduced to two dimensions by reducing the conical frustum to an isosceles trapezoid of height δr and bottom and top edges 2(r ± δr∕2)sin(δα∕2). The length through the model cell along the LOS now corresponds to the length of a line segment between two of the trapezoid edges through the trapezoid centre. This is further simplified by approximating the trapezoid by a rectangle with edges of the lengths δr and rδα, which are the lengths of the trapezoid bimedians. In total, the length of a model cell κ along the LOS is approximated by $$L(\kappa )=\min \left\{\frac{\delta r}{\cos (\alpha (\kappa ))},\frac{r\delta \alpha (\kappa )}{\sin (\alpha (\kappa ))}\right\}.$$(40)

The first approximation (square base to circular base of the truncated pyramid) underestimates the trapezoid base length by a factor of $\sqrt{2}$ at most, resulting in an approximation error of $\sqrt{(2{(\delta r)}^2+{(r\delta \alpha )}^2)/({(\delta r)}^2+{(r\delta \alpha )}^2)}$. The second approximation (trapezoid to rectangle) overestimates the maximum line segment length by less than $\sqrt{{(\delta r)}^2+{(r\delta \alpha )}^2}/(r+\delta r/2)$, giving a total maximum approximation error of $\sqrt{2{(\delta r)}^2+{(r\delta \alpha )}^2}/(r+\delta r/2)$. The error depends on the distance r of the model cell from its centre and is highest for cells close to the model centre and lowest for those on the outer edge. With the values of Table 1, this amounts to consistent 10.5% for λ Cephei, 5.2−8.6% for the heliosphere, 5.2−57% for Proxima Centauri, and 3.5−17% for V374 Pegasi. To reduce the influence of the approximation error, L(κ) was used purely for weighting the cells. The total length of the LOS itself was calculated by adding the radial increment δr to the distance between the first and the last cell along the LOS Δs. Therefore, the result of the LOS integration is $${\text{Res}}_{\text{LOS}}=\left({\displaystyle \sum _{i=1}^N{H}_i}{L}_i\right){\left({\displaystyle \sum _{i=1}^N{L}_i}\right)}^{-1}(\Delta s+\delta r).$$(41)

Both integration methods use the criterion that the angular coordinates (b₀, l₀) of a model cell are inside the respective grid cell. Because the coordinates of a model cell describe the centre of the model cell, the model cell may be positioned so close to the edge of the grid cell that a large part of its volume is outside the grid cell. Likewise, a large part of the simulation cell may be inside the grid cell; as long as the centre is outside the grid cell, the value of its physical parameter is not considered by the integration criterion. It can be assumed, however, that roughly the same number of model cells lie barely outside as inside the cell, so that these effects mostly cancel each other.

Because we need to transform between spherical and Cartesian coordinates, this effect is exacerbated by numerical uncertainties of the cell coordinates. To compensate for this, a grace interval was introduced: if the angular coordinates of the model cell are only slightly outside the respective grid cell, the model cell was still selected for integration. This grace interval was arbitrarily set to 10% of the grid cell extent in every direction. However, this led to multiple counting of model cells for different grid cells, resulting in too large volumes, fluxes, etc. In case of the volumetric integration, this error was corrected by calculating the true volume covered by each LOS, changing the results of Eq. (35) through $${\text{Res}}_{\text{vol}}^{\text{corr}.}=\frac{{\text{Res}}_{\text{vol}}}{{\sum }_{i=1}^N{V}_i}\cdot V_{\text{LOS}},$$(42)

where V_LOS is the calculated volume of each LOS, approximated by an octahedron. The vertices of this octahedron are the points where the LOS intersects a sphere of the simulation radius (R_ISM in Table 2). When the this grace interval was set to zero, then $V_{\text{LOS}}\approx {\sum }_{i=1}^N{V}_i$ with little error, validating this approach. When LOS integration was used instead of volumetric integration, the result was divided by the total number of integrated-over cells and multiplied with the total number of cells, thereby correcting for cells that were selected for multiple lines of sight.

In addition to the integration criterion, other parameters can be used to exclude specific cells: the number density n(κ), the temperature T(κ), and the distance to the origind(κ). This is useful for selecting or disregarding certain regions for integrations, for instance, the outer astrosheath.

The results of the integration are displayed in an orthogonal grid by default, with the displayed grid cells corresponding to the integration grid cells, and the value of the physical parameter represented by colour with a logarithmic colour scale. When the angular extent of the projection is large, the grid can be projected stereographically or in a Hammer projection. In practice, the stereographic projection is not necessary for most cases, while the Hammer projection is a much better depiction of all-sky projections (d₀ = 0) than the orthogonal grid.

Because the simulations only extend over a finite volume, the modelled astrospheres are cut off at the outer border of the simulated domain. Parts of the astrosphere structure were therefore not simulated and not taken into account for the projections. However, for all four examined models, all expected structures lie at least partially inside the simulated volume. Thus, a larger simulated domain is unlikely to yield any new fundamental results but would only quantitatively alter the projections.

Fig. 1 Upper left: sketch of the model being rotated; the homogeneous ISM inflow comes from the positive x-direction. Upper right: sketch of the model being shifted; in most cases, b0 = l0 = 0. Bottom: sketch of the LOS grid covering the projection of a model (red) with the weighting along a single LOS (blue): some cells (grey) are not taken into account because their centres are outside the LOS, while all cells inside the LOS with a common distance from the model centre share the same weight (coloured tiling patterns).

5 Results

The projection methods are applied to the models of λ Cephei in Sect. 5.1 and of the heliosphere in Sect. 5.2. Because the different types of radiation are the most interesting observables that are currently included in the projection methods, we focus on their fluxes. The rotation measure as an observable is also examined. A more limited investigation of the models of Proxima Centauri and V374 Pegasi is presented in Sect. 5.3.

5.1 λ Cephei

In Fig. 2 the mean number density (left panel) and Hα flux (right panel) of an integration area close to the central xz-plane of the λ Cephei model are displayed. The integration area has a thickness along the LOS of 0.6 pc (6% of the total LOS length). Not the entire computed domain is shown. This configuration was chosen to show the structure in a comparably thin box similar to an infinitesimally thin xz-plane. The ISM inflow comes from the left; the undisturbed ISM is the homogeneous area of high number density (10 cm⁻³) and high Hα flux (10⁻¹² erg cm⁻² s⁻¹) left of the parabolic shape, which is the BS. The domain of low number density (10⁻³−10⁻² cm⁻³) and medium Hα flux (10⁻¹⁶−10⁻¹⁴ erg cm⁻² s⁻¹) around (l, b) = (0°, 0°) is the undisturbed SW. The TS and the MD are visible as the outside border of this domain, where the number density jumps from roughly 10⁻³ to 10⁻² cm⁻³ and the Hα flux jumps from roughly 10⁻¹⁷ to 10⁻²³ erg cm⁻² s⁻¹. The AP separates the inner astrosheath (number density ~10⁻² cm⁻³, Hα flux ≤ 10⁻²² erg cm⁻² s⁻¹) from the outer astrosheath (number density (0.1−1) cm⁻³, Hα flux (10⁻¹⁸−10⁻¹⁵) erg cm⁻² s⁻¹), appearing as the colour transition between the TS and the BS. The domain in the inner astrosheath extending to the right from the MD (l > 0.2°) is also called the astrotail.

We display the calculated fluxes of Hα, bremsstrahlung, and cyclotron radiation for the λ Cephei model in Fig. 3. In the left panels, the model centre was placed at the origin (d₀ = 0) and rotated about the z-axis by 90° so that the ISM inflow comes from l = −90°. Different values of the respective physical parameter are mapped onto colour scales, with values lower than the minimum value of the colour scale plotted in that colour. The scaling is bounded to give the best possible representation of all observable features at the cost of cutting off multiple orders of magnitude of the minima, which all occur in the astrotail alone. In all three radiation types, this global minimum is clearly visible at around (b, l) = (0°, 90°), forming a near-perfect circle that appears distorted due to the projection. This corresponds to the astrotail, that is, the direction opposite to the stellar motion, or in the stellar rest frame the direction opposite to the ISM inflow. The astrotail is a region of low density (roughly 10⁻² cm⁻³); therefore all radiation processes are rather weak in this region (cf. Eqs. (7), (16), Sect. 3.2). The global maximum is expected in the opposite direction, where the number density jumps by a factor of 4 at the BS. This is true for the cyclotron radiation (third panel from the top), which is proportional not only to n but also to B² (cf. Sect. 3.2), which also jumps to higher values at the BS. In a roughly circular region with a radius of ~ 30° around the inflow direction (b, l) = (0°, −90°), the cyclotron flux reaches its almost constant maximum value of 7 × 10⁻¹³ erg cm⁻² s⁻¹. This angular extent coincides with the size of the HOA, an additional layer of high-temperature cells (T ≈ 7 × 10⁴ K) in the radial direction that is confined to the central 30° around the inflow axis (the axis through the star parallel to the homogeneous ISM inflow). The high temperature results in high thermal velocities and therefore also in high cyclotron emission. Outside this maximum, two concentric rings of high emission (5 × 10⁻¹³ erg cm⁻² s⁻¹) surround the central maximum, stemming from several regions with high density, high temperature, and/or strong magnetic fields close to the BS. Outside these rings, the cyclotron flux has values of (1− 3) × 10⁻¹³ erg cm⁻² s⁻¹ until it sharply declines due to the astrotail. Unlike the other emission types, the cyclotron flux shows an asymmetry in the astrotail minimum, being elongated to an elliptical shape whose semi-major axis lies at an angle of about 30° to the vertical axis. This is due to the orientation of the magnetic field (cf. Table 2).

While the bremsstrahlung flux (second panel from the top) also shows a central maximum (~ 2 × 10⁻⁷ erg cm⁻² s⁻¹) around the inflow direction, this maximum is local instead of global. A ring of low fluxes (~ 1.5 × 10⁻⁷ erg cm⁻² s⁻¹) with a radius of 30° around the inflow direction separates the central maximum from another high-emission ring ((1−2) × 10⁻⁷ erg cm⁻² s⁻¹) that itself is surrounded by a second low-emission ring. Outside these rings, the bremsstrahlung flux remains at values higher than 1 × 10⁻⁷ erg cm⁻² s⁻¹, reaching its maximum values of (2−3) × 10⁻⁷ erg cm⁻² s⁻¹ close to the astrotail minimum at high latitudes (b ≈±80°). These maxima come from long LOSs through the high-temperature outer astrosheath, while the inner rings correspond to those of the cyclotron flux. The bremsstrahlung flux is not symmetric to the equatorial plane of the model (b = 0°) but shows lower values close to the central maximum in the northern hemisphere ((0.8−1.0) × 10⁻⁷ erg cm⁻² s⁻¹). The high-emission rings described above likewise have lower values in the northern hemisphere than in the southern one.

The global maxima of the bremsstrahlung flux are also global for the Hα flux (top panel, 2 × 10⁻⁷ erg cm⁻² s⁻¹). Unlike the other emission types, there is no central maximum at the inflow direction (b, l) = (0°, −90°) but a local minimum (4 × 10⁻⁸ erg cm⁻² s⁻¹). This minimum has a homogeneous flux inside a circular region of 30° around the inflow direction and a gradient ring of (30−60)° with fluxes up to 1.3 × 10⁻⁷ erg cm⁻² s⁻¹ at the outeredge. This edge is well defined, separating the inner minimum from the higher-emission domain ((1.5−2) × 10⁻⁷ erg cm⁻² s⁻¹) in which the maxima lie. The minimum is caused by the decrease with higher temperatures of the effective radiation rate coefficient α_eff,3→2 (cf. Sect. 3.3), which is therefore much lower in the HOA than in the COA. This feature can also be observed in Fig. 2 (right panel), disrupting the arc of the BS at its vertex.

The various fluxes of the λ Cephei model at distance d₀ = 617 pc (the right three panels of Fig. 3) have a rather similar appearance for each geometric arrangement. In both the xz- (left) and the xy-plane (centre), corresponding to the side and top view of the astrosphere, the parabolic shape of the BS is clearly visible. The projected image is a result of the overlying parabolas of the outer astrosheath from different slices of the computational domain: in the central plane, the BS extends the farthest to the left (smaller l), while the BSs of the slices further from the central plane lie closer to the right (larger l). Because the shifting distance of the parabolas to the right decreases the closer the slice is to the central plane, the positions of the BSs in these slices are close to the position of the leftmost BS in the central slice. Therefore, the LOSs going through the leftmost edge of the parabolic shape encounter multiple BSs from different slices, giving the edge of the parabolic shape its high flux. Because the left and central images look similar, the structure of the astrosphere must be strongly symmetric around the inflow axis. Because the structure of a purely HD model is completely symmetric, it can be concluded that the influence of the magnetic field compared to that of thermal effects must be rather weak, or in other words, that the plasma β must be rather high. This has also been concluded by Scherer et al. (2019) with different arguments. For l > 0.3°, the LOSs of higher longitudes encounter only the astrotail, resulting in fluxes that are orders of magnitude below those from the outer astrosheath or the ISM, therefore causing the global minima. This unphysical effect is due to the finite extent of the model; in a larger model, the outer astrosheath would extend farther outward and would therefore appear in front of and behind the astrotail among these LOS, causing fluxes similar to the region inside the arc structure with 0.0° < l < 0.3°.

The difference between the fluxes from the outer astrosheath (the arc) and the undisturbed ISM (left from the parabolic structure) is strongest for cyclotron radiation (note the different scales for cyclotron radiation and bremsstrahlung/Hα). This is because cyclotron emission is proportional to nB² (cf. Sect. 3.2), with both n and B jumping to higher values at the BS. Bremsstrahlung and Hα emission both follow n² (cf. Eqs. (7) and (16)) and are therefore similar to each other. The Hα flux differs from the other two emission types by the behaviour of the BS at the parabola vertex: While the cyclotron and bremsstrahlung fluxes have their maxima here (taking the abovementioned geometric effect into account), the Hα flux is comparatively low. This is due to the high temperature of the HOA, analogous to the minimum in the Hα projection at d₀ = 0.

The right panels of Fig. 3 show the model projected on the yz-plane, the plane perpendicular to the inflow direction, corresponding to the front view of the astrosphere. In all three fluxes a distinct outer ring of low emission can be observed. This is created by the LOSs that do not encounter any cells from the outer astrosheath, corresponding to the top and bottom lines of the projections on the xy- and xz-planes (central columns of Fig. 3). For Hα, a central minimum is created by the HOA, in analogy to the central minima of the other Hα projections. For the bremsstrahlung and the cyclotron fluxes, faint ring-like structures can be observed, corresponding to the ring-like structures of the projections for d₀ = 0 (left column).

Projections at more angles are displayed in Fig. 4. Hα was chosen as the projected observable because the visible structures are similar for all three radiation types, and also because thisobservable is best suited for comparison with astronomical data. The rotational angle about the z-axis, henceforth labelled α, was varied for these projections in steps of 15° from 15° to 75°. The projection for α =90° corresponds to the projection on the xz- plane and is displayed in Fig. 3 (first row, second column); a projection for α =0° results in the projection on the yz-plane (Fig. 3, first row, fourth column), mirrored about the z-axis. Projections for 90° < α < 180° result in the same images as for 180° − α and for 180° < α < 360° in the mirror images of projections for 360° − α. Figure 4 shows that a crescent of high emission (≥ 3 × 10⁻¹¹ erg cm⁻² s⁻¹) arises to the left of the projection centre with increasing α, which becomes more narrow in width and extends into the right half of the projection until it obtains the distinctive parabolic shape of the bow shock at α = 90°. The area to the left of this crescent is free from any influence of the astrosphere and only consists of the undisturbedhomogeneous ISM. The central minimum from the HOA moves to the left with increasing α until it touches the bow shock crescent, interrupting the arc of high emission at its vertex. The astrotail effect decrease the flux in the area to the right of the crescent, reducing the otherwise mostly homogeneous high emission from the outer astrosheath to lower values at the inflow axis (b = 0) and resulting in the global minimum of the astrotail for α ≥ 60°. Except for the HOA-caused central minimum, the same holds true for the other radiation types.

The rotation measure (RM) of the λ Cephei model is displayed in the bottom row of Fig. 3; the red part (blue part) of the colour scale denominates positive (negative) RMs. The astrotail is clearly visible as the structure of low RMs in the all-sky projection (left panel, RM ≤ 1 rad m⁻²) and also on the far right in the projections on the xz- and xy-plane (central panels, RM ≤ 0.1 rad m⁻²). The alternating structure of positive and negative RM bands in the all-sky projection mostly comes from the homogeneous magnetic field of the undisturbed ISM. The BS can clearly be recognized in the xz- and xy-planes. In the yz-plane projection, only the radial gradient can be recognized as a structure.

Fig. 2 Mean number density (left) and Hα flux (right) of a flat box at the λ Cephei model central xz-plane with a thickness of 0.6 pc.

Fig. 3 Projections of the λ Cephei model in (from top to bottom) Hα flux, bremsstrahlung flux, cyclotron flux, and rotation measure at distances d0 = 0 (left column) and d0 = 617 pc (right columns), projected onto (from left to right) the sky and the xz-, xy-, and yz-planes. Values lower than the minimum value of the colour scales are plotted in the colour of that minimum value.

Fig. 4 Further projections of the λ Cephei model in the Hα flux, rotated about the x-axis by 0°, and about the z-axis by (from left to right) 15°, 30°, 45°, 60°, and 75°. A rotation of90° (180°) results in thesame projection on the xz- (yz-) plane as in Fig. 3. Scales are the same as in Fig. 3 (top row).

Fig. 5 Hα flux of the heliosphere model at distances d0 = 0 (left column) and d0 = 1 pc (right columns), projected onto (from left to right) the sky and the xz-, xy-, and yz-planes. Values lower than the minimum value of the colour scales are plotted in the colour of that minimum value.

5.2 Heliosphere

The heliosphere model was evaluated by the same method. The projections of all three radiation types look very similar, therefore we only display the Hα flux in Fig. 5. Unlike the λ Cephei model, the ISM fluid speed is too low for a BS to emerge, at least in single-fluid modelling (see, however, Scherer & Fichtner 2014). Instead of a sudden jump, the ISM parameters increase or decrease gradually in the entire region outside the heliopause (HP), forming a so-called bow wave (McComas et al. 2012; Zank et al. 2013).

The left panel shows the all-sky projection with the model centre at origin (d₀ = 0). As with the λ Cephei model, the heliotail is clearly visible as the global minimum. However, unlike the λ Cephei model, the tail of the heliosphere model is elongated in the b-direction. This stems from an asymmetric shape of the heliosphere: the average plasma β is much lower for the heliosphere model than for the λ Cephei model, therefore the effect of the magnetic field is much stronger. Because there is no BS, there also is no (hot) outer heliosheath, and the flux accordingly follows a simple gradient outside the heliotail instead of showing a ring-like structure around the bow direction. The lack of a BS is much more apparent in the projections at larger distances (d₀ = 1 pc, right panels), where a radial gradient corresponding to the bow wave is apparent in the xz-plane until the inflow hits the heliopause. The projection on the xy-plane shows an indentation on the right side, caused by the low emission of the heliotail. As in the λ Cephei model, the heliotail causes a ring-like structure in the projection on the yz-plane, although the ring is elongated in the b-direction due to the asymmetry of the heliotail.

The RM in the all-sky projection shows a similar structure of alternating bands of positive and negative RMs as the λ Cephei model, only with much lower values (RM < 10⁻⁷ rad m⁻²). Like the astrotail of the λ Cephei model, the heliotail causes a structure of lower RMs; no further structures can be seen. It is unlikely that such low values of the RM could be detected because fluctuations of the density and the magnetic field in the ISM outside the model range are assumed to cause RM values of similar magnitude.

Fig. 6 Left: Bremsstrahlung flux of the Proxima Centauri model at a distance of 1 pc, projected at the xz-plane. Right: cyclotron flux of the central region of the V374 Pegasi model at a distance of 9 pc, projected at the xy-plane. The cyclotron flux is scaled linearly.

5.3 Proxima Centauri and V374 Pegasi

The speed of Proxima Centauri relative to the ISM is high enough to produce a BS; the entire astrotail is contained in the model. The radiation fluxes are of the same order of magnitude as in the heliosphere model. Like in the λ Cephei model, the origin of the highest fluxes lies in the outer astrosheath, where the density and the magnetic induction are comparably high (n ≲ 10⁻² cm⁻³, B ≲ 10⁻⁵ G) and the temperature is moderate (T ≲ 10⁵ K). The RM is about one magnitude lower than the RM of the heliosphere model and shows a similar large-scale structure in the all-sky projection as the previous two models. Because the astrotail is completely contained in the model, there is only a faint minimum in the direction opposite to the inflow. As an example, the bremsstrahlung flux in the xz-plane is shown in Fig. 6 (left); the other fluxes are similar.

No BS appears to be produced in the V374 Pegasi model, although the structure outside the TS is complicated and requires further investigation. The projections show an elongated structure in the cyclotron flux of the xy-plane (see Fig. 6 right) that is also faintly visible in the other projection planes and observables. The fluxes are higher than those of the heliosphere model, but still far lower than the detection limits. It is noticeable, however, that the bremsstrahlung flux is higher by three orders of magnitude while the Hα flux is only higher by one order magnitude than that of the heliosphere model; this is due to the high temperatures of the V374 Pegasi model (see Table 2), which lead to low effective recombination rate coefficients. Because of the high magnetic induction (3 × 10⁻⁵ G outside the TS), the cyclotron flux is higher than that of the heliosphere model by six orders of magnitude. The critical frequency of bremsstrahlung is ν_c = 10¹⁸ Hz for the V374 Pegasi model, which is higher by two orders of magnitude than ν_c of the other three models; this is also caused by the high temperature (T ≳ 10⁶ K outside the TS) of the model.

6 Observability

To evaluate the model projections, the computed fluxes are first compared to current observational limits in Sect. 6.1. Subsequently, in Sect. 6.2 observational images of λ Cephei are compared to the corresponding projection.

6.1 Observational limits

The detectable minimum flux of Hα emission was assumed to be 10⁻¹⁹ erg cm⁻² s⁻¹ arcsec⁻² (Donahue et al. 1995). The angular area of the projection grid cells for the λ Cephei model at d₀ = 0 (left column) is 1.14 × 10⁸ arcsec^2 = 2.68 × 10⁻³ sr; at d₀ = 617 pc (right columns) it is 5.05 × 10⁴ arcsec^2. Multiplying the detection limit with the respective angular area for d₀ = 0 gives a minimum flux of 1.14 × 10⁻¹¹ erg cm⁻² s⁻¹; for the d₀ = 617 pc projections, the minimum flux amounts to 5.05 × 10⁻¹⁵ erg cm⁻² s⁻¹. The entire projections of the λ Cephei model would be clearly detectable.

The emission spectrum of bremsstrahlung was assumed to be flat up to the critical frequency ν_max = 10¹⁶ Hz, so that the spectral flux was calculated from the total flux by dividing by ν_max. Two estimates for the detectable minimum flux were made, one of which is given by the NLAO VLA Sky Survey (NVSS) with a rms brightness fluctuation of σ = 0.45 mJy∕(4° × 4°) ≈ 2.2 × 10⁻³⁵ erg s⁻¹ cm⁻² Hz⁻¹ arcsec⁻² at a frequency of 1.4 GHz (Condon et al. 1998). In combination with the solid angles described above, this amounts to detection limits of 2.5 × 10⁻²⁷ erg s⁻¹ cm⁻² Hz⁻¹ at d₀ = 0 and 1.1 × 10⁻³⁰ erg s⁻¹ cm⁻² Hz⁻¹ at d₀ = 617 pc. These values are below the calculated bremsstrahlung fluxes by some orders of magnitude. Spectrum-modifying effects such as self-absorption are most effective for low frequencies, therefore we assumed that the observed bremsstrahlung flux may be lower than the emitted flux by many ordersof magnitude. Therefore, an estimate of the detectable minimum flux at higher frequencies is of value. A detection limit of 10⁻¹⁵ erg s⁻¹ cm⁻² is given for a spectral range $[0.5,2]\text{keV}\widehat{=}[1.2,4.8]\times {10}^{17}\text{Hz}$ by Hasinger (2004). The corresponding spectral flux is 2.8 × 10⁻³³ erg s⁻¹ cm⁻² Hz⁻¹ and lower than the calculated values by several orders of magnitude. While the detection limit is given for a spectral range one order of magnitude higher than the calculated critical frequency, we assumed that the detection limit of the correct spectral range is in a similar order of magnitude. This means that the entire projection of the λ Cephei model would be detectable.

Becausecyclotron radiation is emitted at extremely low frequencies, the minimum detectable flux of the LOw Fequency ARray (LOFAR, van Haarlem et al. 2013) was used as a detection limit. LOFAR is sensitive in a spectral range [10, 240] MHz; a sensitivity of σ_{60 MHz} = 20 μJy is given for frequency ν = 60 MHz. The sensitivity is frequency-dependent with σ_ν ∝ ν^−0.7. Extrapolated to a frequency of 300 Hz, the detection limit would be σ_{300 Hz} = 103 mJy = 1.03 × 10⁻²⁴ erg s⁻¹ cm⁻² Hz⁻¹. While neither LOFAR nor any other astronomical instrument is sensitive in this spectral range, we still used σ_{300 Hz} as a minimum detectable flux for the sake of comparison. Compared to the calculated cyclotron fluxes, σ_{300 Hz} is lower by a few orders of magnitude. If an instrument with detection limits comparable to LOFAR were sensitive for the λ Cephei model gyration frequencies (ν_g < 500 Hz), it would be able to detect such astrospheres with ease. However, cyclotron (or synchrotron) radiation can be safely discarded as an origin to the detected emission of astrospheres because the frequencies are far too low.

For the heliosphere model, the fluxes are all much lower than those of the λ Cephei model, about eight to nine orders of magnitude for all three types of radiation. The gyration frequencies, at which cyclotron radiation is emitted, are lower than 20 Hz, while the critical bremsstrahlung frequency remains at ν = 10¹⁶ Hz. The angular area of the projection cells is 1.14 × 10⁸ arcsec² for d₀ = 0 and 1.90 × 10³ arcsec² for d₀ = 1 pc. The resulting detection limits of Hα are 1.14 × 10⁻¹¹ arcsec² (d₀ = 0) and 1.90 × 10⁻¹⁶ arcsec² (d₀ = 1 pc), which is far too high for any chance of detection. Similarly, the detection limits of bremsstrahlung and cyclotron fluxes are higher than the calculated fluxes by many orders of magnitude.

Because the fluxes for the model of Proxima Centauri are on the same order of magnitude as those of the heliosphere model, the model astrosphere would not be detectable with the assumed detection limits. While the fluxes for the model of V374 Pegasi are higher than those of the heliosphere model, they still remain far lower than the detection limits.

Fig. 7 Composite images of the shock structure around λ Cephei. The star is at the central position, marked by a magenta cross; north is up, east is left. The size of both plotted fields is 20′ × 20′. Left: WISE W3 band (12 μm, red), WISE W4 band (22 μm, green), IPHAS Hα (blue); centre: MIPS 24 μm (greyscale), IPHAS Hα (magenta). Right: composite image of the simulated shock structure, rotated about the z-axis by 25°, and the image from the centre panel. The simulation image was further rotated by 90° + 31.373° anticlockwise to align the BS directions.

6.2 Observational data

Composite images of the area around λ Cephei are displayed in Fig. 7 by combining observational data from three infrared bands, observed by the Wide-field Infrared Survey Explorer (WISE, left panel, Wright et al. 2010) and the Multiband Imaging Photometer for Spitzer (MIPS, centre panel, Werner et al. 2004; Rieke et al. 2004), and Hα, obtained as part of the INT Photometric Hα Survey of the Northern Galactic Plane (IPHAS, left and centre panels, Drew et al. 2005; Barentsen et al. 2014). The star moves in the south-southwest direction (the bottom right direction in the image plane; east is left and west is right), so that the bow shock is faintly visible as a wide arc southeast to west of the star in the long-wave infrared. A small part of the bow shock is clearly visible in both Hα and the infraredsoutheast of the star. While this structure is very apparent in Hα and the long-wave infrared, it is faint at 12 μm. The large disc-like structure in Hα centred on λ Cephei stems from ghosting effects and is purely instrumental. A much larger structure closer to the star in the southwest, visible at the edge of the ghosting disc, can be seen in the infrared bands but not in Hα.

To compare observed and simulated astrospheres, the model was projected with observed geometry of λ Cephei: Gaia measures a parallax angle of (1.6199 ± 0.1265) mas (Gaia Collaboration 2018), resulting in a distance of roughly d = 617 pc; Bouret et al. (2012) calculated a distance of 950 pc. Gaia measures a proper motion of μ_α cos(δ) = −6.328 mas yr⁻¹, μ_δ = −10.378 mas yr⁻¹ at right ascension α =103.83° and declination δ = 2.61°. The radial velocity is v_rad = −75.10 km s⁻¹ (Gontcharov 2006), pointing towards the observer. Assuming the distance derived from Gaia measurements, the angle of motion, that is, the angle between proper and radial motion, which in the simulation is the rotation about the z-axis, is $$\gamma =\arctan \left(\frac{d}{v_{\text{rad}}}\sqrt{{\mu }_{\alpha }^2{\cos }^2(\delta )+{\mu }_{\delta }^2}\right)\approx {25.35}^{°}.$$(43)

With an assumed distance of 950 pc, the angle of motion would be 36.11°. The proper motion angle, that is, the angle of the proper motion to the image coordinate system, is $\varphi =\arctan \left({\mu }_{\alpha }\cos (\delta )/{\mu }_{\delta }\right)\approx {31.373}^{°}$.

The projection of the simulated astrosphere has a distance between BS and star of about 0.2° = 12′, whereas the observed distance between BS and star is about 8′. Because the size of the simulated astrosphere scales with the boundary conditions, especially the inner number density and fluid velocity, the difference in scales is most likely caused by incorrect estimates of these values. The overall shape of the simulated astrosphere does not depend on the boundary conditions as strongly; scaling the images so that the bow shocks align is therefore a valid approach.

This was done for the right panel of Fig. 7: the MIPS composite image (central panel) was combined with a projection of the simulated astrosphere. This projection was generated with an angle about the z-axis of γ = 25° (see Eq. (43)); the image was further rotated by 90° + ϕ = 121.373° anticlockwise to align the inflow direction of the simulated astrosphere with the stellar proper motion. The projection image was scaled so that the BSs appear at the same distance to the star. The shapes of the BSs align well: the filament southeast of the star (bottom leftquadrant) can be identified as part of the BS. The unknown infrared emission structure inside the BSs roughly coincides with the astropause, but the inner structure of the astrosphere unfortunately disappears inside the Hα ghosting disc.

7 Conclusions

We presented and evaluated projections of 3D single-fluid MHD simulations of the astrospheres of λ Cephei and the heliosphere. The projections of two further astrosphere models, those of Proxima Centauri and V374 Pegasi, werediscussed. It was shown that the region between the astropause and the bow shock, the outer astrosheath,causes high fluxes in all examined radiation processes, resulting in a variety of shapes following the apparent shape of the outer astrosheath. This shape depends on the inclination angle between the LOS and the direction of the stellar motion, resulting in a ring with a central maximum when the angle is close to 0°, a filled parabola when the angle is close to 90°, and a crescent with a barely detached or attached central maximum when the angle is in between. Likewise, high values of the rotation measure follow the surface of the bow shock in the central plane, leading to parabolic shapes. This is consistent with the results of different approaches, such as the 3D HD models by Mohamed et al. (2012) or the 2D MHD models by Meyer et al. (2017), the latter of which also gives similar maps of their axisymmetric model at different inclination angles.

In astrospheres without a bow shock, for instance, a single-fluid model of the heliosphere, this region does not exist. No distinguishable shapes are visible in the fluxes or the rotation measure. It is unlikely that detailed observations of such astrospheres can be made. While the astrospheres of low-mass stars with supersonic motion relative to their ISM (such as Proxima Centauri) produce a bow shock and display comparably high fluxes in their outer astrosheaths, these fluxes are still much lower than those of the astrosheaths of runaway high-mass stars such as λ Cephei. Even though many more of these low-mass stars are far closer to us than high-mass stars, their fluxes are most likely too low to be detectable.

Unlike with other approaches (e.g. del Valle & Pohl 2018), no synchrotron radiation could be obtained with these models. This is due to the non-relativistic thermal speed of the electrons, which do not exceed temperatures of 10⁸ K. For relativistic electrons, which can be acquired by higher initial temperatures and (less likely by) different heating terms in the governing MHD equations, synchrotron radiation must be expected. While cyclotron radiation is emitted at very low and ultra-low frequencies, far below the spectral range of observational instruments, the bremsstrahlung spectrum reaches frequencies in the extreme ultraviolet and soft X-ray regime. Higher thermal speeds of electrons would result in even higher frequencies, as were observed by Sánchez-Ayaso et al. (2018), for example.

Acknowledgements

K.S. is grateful to the Deutsche Forschungsgemeinschaft (DFG), funding the project SCHE334/9-2. J.K. acknowledges financial support through the Ruhr Astroparticle and Plasma Physics (RAPP) Center, funded as MERCUR project St-2014-040. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This paper makes use of data obtained as part of the INT Photometric Hα Survey of the Northern Galactic Plane (IPHAS, www.iphas.org) carried out at the Isaac Newton Telescope (INT). The INT is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. All IPHAS data are processed by the Cambridge Astronomical Survey Unit, at the Institute of Astronomy in Cambridge. The bandmerged DR2 catalogue was assembled at the Centre for Astrophysics Research, University of Hertfordshire, supported by STFC grant ST/J001333/1.

References

All Tables

All Figures

Fig. 1 Upper left: sketch of the model being rotated; the homogeneous ISM inflow comes from the positive x-direction. Upper right: sketch of the model being shifted; in most cases, b0 = l0 = 0. Bottom: sketch of the LOS grid covering the projection of a model (red) with the weighting along a single LOS (blue): some cells (grey) are not taken into account because their centres are outside the LOS, while all cells inside the LOS with a common distance from the model centre share the same weight (coloured tiling patterns).

In the text

	Fig. 2 Mean number density (left) and Hα flux (right) of a flat box at the λ Cephei model central xz-plane with a thickness of 0.6 pc.
In the text

	Fig. 3 Projections of the λ Cephei model in (from top to bottom) Hα flux, bremsstrahlung flux, cyclotron flux, and rotation measure at distances d0 = 0 (left column) and d0 = 617 pc (right columns), projected onto (from left to right) the sky and the xz-, xy-, and yz-planes. Values lower than the minimum value of the colour scales are plotted in the colour of that minimum value.
In the text

	Fig. 4 Further projections of the λ Cephei model in the Hα flux, rotated about the x-axis by 0°, and about the z-axis by (from left to right) 15°, 30°, 45°, 60°, and 75°. A rotation of90° (180°) results in thesame projection on the xz- (yz-) plane as in Fig. 3. Scales are the same as in Fig. 3 (top row).
In the text

	Fig. 5 Hα flux of the heliosphere model at distances d0 = 0 (left column) and d0 = 1 pc (right columns), projected onto (from left to right) the sky and the xz-, xy-, and yz-planes. Values lower than the minimum value of the colour scales are plotted in the colour of that minimum value.
In the text

	Fig. 6 Left: Bremsstrahlung flux of the Proxima Centauri model at a distance of 1 pc, projected at the xz-plane. Right: cyclotron flux of the central region of the V374 Pegasi model at a distance of 9 pc, projected at the xy-plane. The cyclotron flux is scaled linearly.
In the text

Fig. 7 Composite images of the shock structure around λ Cephei. The star is at the central position, marked by a magenta cross; north is up, east is left. The size of both plotted fields is 20′ × 20′. Left: WISE W3 band (12 μm, red), WISE W4 band (22 μm, green), IPHAS Hα (blue); centre: MIPS 24 μm (greyscale), IPHAS Hα (magenta). Right: composite image of the simulated shock structure, rotated about the z-axis by 25°, and the image from the centre panel. The simulation image was further rotated by 90° + 31.373° anticlockwise to align the BS directions.

In the text