Issue 
A&A
Volume 576, April 2015



Article Number  A97  
Number of page(s)  7  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201425091  
Published online  09 April 2015 
Cosmic rays in astrospheres
^{1}
Institut für Theoretische Physik IV: Weltraum und Astrophysik,
RuhrUniversität Bochum, 44780
Bochum,
Germany
email:
kls@tp4.rub.de; hf@tp4.rub.de; jk@tp4.rub.de; tow@tp4.rub.de
^{2}
Research Department, Plasmas with Complex Interactions,
RuhrUniversität Bochum, 44780
Bochum,
Germany
^{3}
Center for Space Research, NorthWest University,
2520
Potchefstroom, South
Africa
email: 12834858@nwu.ac.za; Stefan.Ferreira@nwu.ac.za; DuToit.Strauss@nwu.ac.za
^{4}
Astronomisches Institut, RuhrUniversität Bochum,
44780
Bochum,
Germany
email: bomans@astro.rub.de; kweis@astro.rub.de; thomas.wodzinski@astro.rub.de
Received: 1 October 2014
Accepted: 10 February 2015
Context. Cosmic rays passing through large astrospheres can be efficiently cooled inside these “cavities” in the interstellar medium. Moreover, the energy spectra of these energetic particles are already modulated in front of the astrospherical bow shocks.
Aims. We study the cosmic ray flux in and around λ Cephei as an example for an astrosphere. The largescale plasma flow is modeled hydrodynamically with radiative cooling.
Methods. We study the cosmic ray flux in a stellar wind cavity using a transport model based on stochastic differential equations. The required parameters, most importantly, the elements of the diffusion tensor, are based on the heliospheric parameters. The magnetic field required for the diffusion coefficients is calculated kinematically. We discuss the transport in an astrospheric scenario with varying parameters for the transport coefficients.
Results. We show that large stellar wind cavities can act as sinks for the Galactic cosmic ray flux and thus can give rise to smallscale anisotropies in the direction to the observer.
Conclusions. Smallscale cosmic ray anisotropies can naturally be explained by the modulation of cosmic ray spectra in huge stellar wind cavities.
Key words: stars: winds, outflows / cosmic rays / hydrodynamics
© ESO, 2015
1. Introduction
Recently, simulations of astrospheres around hot stars have gained new interest, see for example Decin et al. (2012), Cox et al. (2012), Arthur (2012), van Marle et al. (2014). These authors modeled astrospheres using a (magneto)hydrodynamic approach, either in 1D or 2D. In this work, such astrospheric models are used for the first time to estimate the cosmic ray flux (CRF) through it. Because of the large spatial extent of O star astrospheres (wind bubbles), these objects can efficiently cool the spectrum of Galactic cosmic rays (GCR). Especially λ Cephei is an interesting example, being the brightest runaway Of star in the sky (type O6If(n)p). We estimate the CRF at different energies for λ Cephei as an example of an O star astrosphere using stochastic differential equations (SDE; Strauss et al. 2013) to solve the GCR transport equation. Runaway O and B stars are common and part of a sizable population in the Galaxy, a significant number of which show bowshock nebulae (e.g., Huthoff & Kaper 2002). In Sect. 3 we discuss the radiative cooling functions, while in Sect. 4 we show the astrosphere model results. In Sect. 5 we estimate the CRF.
2. Largescale structure of astrospheres
Winds around runaway stars, or in general, stars with a nonzero relative speed with respect to the surrounding interstellar medium (ISM), develop bulletshaped astrospheres.
The hydrodynamic largescale structure is sketched in Fig. 1, the notation of which is described below.
The hypersonic stellar wind (Mach numbers Ma ≫ 1) undergoes a shock transition to subsonic velocities at the termination shock (TS) in the inflow direction. Then a tangential discontinuity, the astropause (AP), is formed between the ISM and the stellar wind, where the velocity normal to it vanishes: there is no mass transport through the AP. Other quantities such as the tangential velocity, temperature, and density are discontinuous, while the thermal pressure is the same on both sides. If the relative speed, or the interstellar wind speed as seen upwind in the rest frame of the star, is supersonic in the ISM, a bow shock (BS) exists. If the relative speed is subsonic, there will be no BS, see Table 1 for the stellar parameters and for the stellarcentric model distances of the TS, AP, and BS. The region between the BS and AP is called outer astrosheath, the region between the AP and TS the inner astrosheath. The AP around the inflow direction at the stagnation line is sometimes called the nose, while the region beyond the downwind TS is called the astrotail. The latter can extend deep into the ISM. The region inside the TS is called the inner astrosphere.
Fig. 1 Largescale structure of an astrosphere. For details see text. 
Parameters for λ Cephei (van Leeuwen 2007).
Stellar and interstellar boundary conditions.
In the downwind direction, the termination shock forms a triple point (T), from which the Mach disk (MD) extends down to the stagnation line; this is the line through the stagnation point and the star. A tangential discontinuity (TD) emerges down into the tail. A reflected shock (not shown here) also extends from the triple point toward the TD.
This is the standard shape of an astrosphere using a single hydrodynamic fluid (Baranov et al. 1971; Pauls et al. 1995). For the large dimensions of O star astrospheres, cooling operates inside the outer astrosheath, but usually not in the inner astrosheath. Cooling is also present beyond the model boundary (see Table 2) of the inner astrosphere, which leads to relative sizes of these regions different from that of pure hydrodynamic flow.
The astrosphere is always bulletshaped, which can be seen by the conservation of momentum: ${\mathit{\rho}}_{\mathrm{sw}}{\mathit{v}}_{\mathrm{sw}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{sw}}\mathrm{=}{\mathit{\rho}}_{\mathrm{ISM}}{\mathit{v}}_{\mathrm{ISM}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{ISM}}$, where ρ,v,P are the density, velocity, and the thermal pressure of the stellar wind (subscript sw) and the ISM (subscript ISM). The stellar wind pressure is usually neglible inside the TS. For hypersonic flows this holds true for the thermal ISM pressure in inflow direction, while for subsonic inflows the thermal pressure dominates. In the tail direction, only the thermal ISM presssure is present beyond the TS. This means that even if the inflow is subsonic, there is a total pressure asymmetry between upwind and downwind as long as the ISM velocity does not vanish. Thus on long time and large spatial scales, a bulletshaped astrosphere will develop.
3. Cooling and heating
As a result of the large dimensions of O star astrospheres and their outer astrosheath, the plasma can effectively be cooled by Coulomb collisions. In this process, an electron in an atom (molecule) will be excited, and after returning to a lower energy state, the reemitted photon will carry away the energy (Sutherland & Dopita 1993,and references therein). This can lead to an effective cooling of the shocked ISM (ISM_{2} in Fig. 1).
Several cooling functions are discussed in the literature (e.g., Rosner et al. 1978; Mellema & Lundqvist 2002; Townsend 2009; Schure et al. 2009; Reitberger et al. 2014). The differences in the cooling functions are caused by different abundances and different levels of approximation. Because we do not know the abundances in the ISM surrounding λ Cephei, we use here the analytic representation by Siewert et al. (2004), which lies in between all the other cooling functions mentioned above, and can therefore be considered as a useful principal representation.
Mainly the hot shocked ISM is affected, while the stellar wind gas is not. The reason for this is that the number densities in the shocked ISM are higher by a factor of 3.6, which is the compression ratio between the shocked and unshocked ISM, and that it is by a factor 10^{2} hotter, that is, T_{ISM,shocked} ≈ 2.5 × 10^{5} K. In this region most of the mentioned cooling functions yield similar values. Inside the inner astrosheath, that is, the region between the AP and the TS, the number density is on the order of n = 10^{3} cm^{3} and the velocities are on the order of v = 600 km s^{1}, and thus the cooling length scale, depending on v and n, becomes huge and is neither important there nor inside the TS, see below.
The heating length scale only depends on the number density n and also increases to scales much larger than the distances inside the AP. Therefore, as a result of the huge ram pressure of the stellar wind inside the TS, we can, from a dynamical point of view, safely neglect heating and cooling, because it does not influence the adiabatic expansion of the stellar wind. Moreover, the thermal pressure is negligible compared with the ram pressure in the inner astrosphere. With the help of the momentum equation, we can uniquely determine the shocked thermal pressure, which dominates in the inner astrosheath.
In the following we always assume quasineutrality n_{e} = n_{p} ≡ n, where n_{e},n_{p} are the electron and proton number densities, respectively. Furthermore, we neglect the contribution of heavy ions. In Table 1 we summarize the parameters of λ Cephei derived from observations and the characteristic distances of its model astrosphere. The parallax translates into a tangential velocity of 41.1 km s^{1}, which together with the radial velocity of 75.1 km s^{1} gives a total speed of 85.5 km s^{1}. Our hydrodynamic model is threedimensional (for later use), and in view of the uncertainties of the ISM state, we have chosen a set of parameters as given in Table 2. The standard procedure in modeling is to choose one axis as the inflow direction (here the y axis). As long as the ISM is homogeneous, with a vanishing magnetic field, and the stellar wind is spherically symmetric, the astrosphere is symmetric around the inflow direction.
The derived inner boundary conditions for the model are estimated from the massloss rate and the terminal velocity and taken at 0.03 pc. They are presented in Table 3 together with those of the ISM.
Some characteristic numbers using the cooling function from Siewert et al. (2004).
As stated above, the cooling acts differently in the subsonic regions from the way it does in the supersonic regions because in the former the thermal pressure P is dominant, while in the latter it is the ram pressure ρv^{2}/ 2, where ρ and v are the mass density (mainly protons, but helium or other elements may contribute to the mass density) and bulk speed, respectively. While in the following we discuss only protons, we can add, for example, helium, which leads to partial densities, temperatures, and pressures for which the approximations made below can be separated. We are interested in the shortest characteristic length, which is in the subsonic case inversely proportional to the number densities, and thus the protons dominate.
In the hypersonic case, the estimation below can differ by up to 40% when including helium. This would then also require including helium as a new species in the Euler energy equations (including different species, see Scherer et al. 2014), which would violate our assumption of a single fluid. The interaction of other species concerning the energy loss by Coulomb collisions is, however, already included in the cooling functions, so that we can continue with the singlefluid equations consisting of protons including the cooling term for a first analysis.
From the stationary energy equation we obtain in the subsonic region with the assumption nm_{p}v^{2}/ 2 ≪ P: $\begin{array}{ccc}{\nabla}\mathrm{\xb7}\left(\frac{\mathit{\gamma}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\mathit{P}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{n}{\mathit{m}}_{\mathrm{p}}{\mathit{v}}^{\mathrm{2}}\right){v}& \mathrm{\approx}& \frac{\mathit{\gamma Pv}}{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{L}}_{\mathrm{cool}\mathit{,s}}}\\ & \mathrm{=}& \frac{\mathrm{5}\mathit{nkTv}}{{\mathit{L}}_{\mathrm{cool}\mathit{,s}}}\mathrm{=}\mathrm{}{\mathit{n}}^{\mathrm{2}}\mathrm{\Lambda}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathit{,}\end{array}$(1)whereγ = 5 / 3 is the adiabatic index, m_{p} the proton mass, P = 2nkT, k is the Boltzmann constant, v,T the bulk velocity and temperature of the plasma flow, and Λ the cooling function. Taking the absolute values in Eq. (1) and replacing ∇ by the inverse subsonic cooling length, L_{cool,s} can be estimated as $\begin{array}{ccc}{\mathit{L}}_{\mathrm{cool}\mathit{,s}}\mathrm{\approx}\frac{\mathit{\gamma Pv}}{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{n}}^{\mathrm{2}}\mathrm{\Lambda}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{5}\mathit{kTv}}{\mathit{n}\mathrm{\Lambda}\mathrm{\left(}\mathit{T}\mathrm{\right)}}& & \end{array}$(2)and the subsonic cooling time $\begin{array}{ccc}{\mathit{\tau}}_{\mathrm{cool}\mathit{,s}}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{cool}\mathit{,s}}}{\mathit{v}}\mathrm{=}\frac{\mathrm{5}\mathit{kT}}{\mathit{n}\mathrm{\Lambda}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{\xb7}& & \end{array}$(3)This cooling time is up to a factor 3 the same as that given in Sutherland & Dopita (1993). For supersonic interstellar flows, that is, ${\mathit{\rho}}_{\mathrm{ISM}}{\mathit{v}}_{\mathrm{ISM}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{\gg}\mathit{\gamma}\mathit{/}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}\mathit{P}$, it follows: $\begin{array}{ccc}{\nabla}\mathrm{\xb7}\left(\frac{\mathit{\gamma}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\mathit{P}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{n}{\mathit{m}}_{\mathrm{p}}{\mathit{v}}^{\mathrm{2}}\right){v}\mathrm{\approx}\mathrm{\mathcal{O}}\left(\frac{{\mathit{m}}_{\mathrm{p}}\mathit{n}{\mathit{v}}^{\mathrm{3}}}{\mathrm{2}{\mathit{L}}_{\mathrm{cool}\mathit{,h}}}\right)\mathit{,}& & \end{array}$(4)and we derive the supersonic (index h for “hypersonic” to distinguish it from the subsonic one) cooling length L_{cool,h} and time τ_{cool,h}$\begin{array}{ccc}{\mathit{L}}_{\mathrm{cool}\mathit{,h}}\mathrm{\approx}\frac{{\mathit{m}}_{\mathrm{p}}{\mathit{v}}^{\mathrm{3}}}{\mathrm{2}\mathit{n}\mathrm{\Lambda}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\u2001\mathrm{and}\u2001{\mathit{\tau}}_{\mathrm{cool}\mathit{,h}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{p}}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}\mathit{n}\mathrm{\Lambda}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{\xb7}& & \end{array}$(5)For the heating function by photoionization and some supplemental heat source we use the approach by Reynolds et al. (1999), see also Kosiński & Hanasz (2006). The heating rate for photoionization depending on electron collisions is limited by recombination and is thus proportional to ${\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}$, while additional heating terms that can include photoelectric heating by dust, dissipation of turbulence, interactions with cosmic rays are proportional to n_{e}: $\begin{array}{ccc}\mathrm{\Gamma}\mathrm{=}{\mathit{n}}^{\mathrm{2}}{\mathit{G}}_{\mathrm{0}}\mathrm{+}\mathit{n}{\mathit{G}}_{\mathrm{1}}\mathit{,}& & \end{array}$(6)with the constants G_{0} = 10^{24} erg cm^{3} s^{1} and G_{1} = 10^{25} erg s^{1} (Kosiński & Hanasz 2006). Replacing the righthand side of Eq. (1) by Γ, we obtain $\begin{array}{ccc}{\mathit{L}}_{\mathrm{heat}\mathit{,s}}\mathrm{=}\frac{\mathrm{5}\mathit{kTv}}{\mathit{n}{\mathit{G}}_{\mathrm{0}}\mathrm{+}{\mathit{G}}_{\mathrm{1}}}\u2001\mathrm{and}\u2001{\mathit{L}}_{\mathrm{heat}\mathit{,h}}\mathrm{=}\frac{{\mathit{m}}_{\mathit{p}}{\mathit{v}}^{\mathrm{3}}}{\mathrm{2}\mathrm{(}\mathit{n}{\mathit{G}}_{\mathrm{0}}\mathrm{+}{\mathit{G}}_{\mathrm{1}}\mathrm{)}}\mathit{,}& & \end{array}$(7)and the heating times τ_{heat,s},τ_{heat,h} by dividing the respective heating lengths by the speed v.
We can estimate the cooling and heating lengths and times for the shocked ISM (ISM_{2}), and analogously for the shocked stellar wind SW_{2} (see Fig. 1). The results are displayed in Table 3 together with the characteristic lengths and times for the hypersonic ISM (ISM_{1}), because it is also cooled. For the ISM the dependence of the cooling and heating lengths is shown in Fig. 2. The vertical line denotes the temperature at which the ram pressure equals the thermal pressure. Left of this line the hypersonic length scales from Eqs. (5) and (7) are displayed, while on the right the subsonic scales are shown (Eqs. (2) and (7)).
The supersonic parameters are number density n = 11 cm^{3} and a speed of v = 80 km s^{1} while for the subsonic parameters, the density was multiplied by the compression ratio s = 3.62 and the speed was divided by s using the RankineHugoniot relations.
Figure 2 shows that in the subsonic case for temperatures above ≈5 × 10^{5} K the heating scale length is always longer than that for the cooling. Thus cooling is more efficient than heating for the discussed functions and parameters. In the supersonic case (T< 5 × 10^{5} K), the heating scale lengths are only longer down to temperatures of ≈10^{4} K, while for temperatures below ≈4 × 10^{2} K the cooling scale length is more important. Thus, above ≈4 × 10^{2} K cooling is more efficient than heating, and below this, heating is important. We can read from this figure that the length scale for cooling of the shocked ISM temperature T ≈ 2.5 × 10^{5} K is ≈0.01 pc. From the model we see that the BS to AP distance is 0.17 pc, which is more than ten times the cooling length L_{cool,s}. This strong cooling is balanced by the fact that at least at the stagnation line the shocked ram pressure $\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}{\mathit{\rho}}_{\mathrm{2}\mathit{,}\mathrm{ISM}}{\mathit{v}}_{\mathrm{2}\mathit{,n,}\mathrm{ISM}}^{\mathrm{2}}$ has to be converted into thermal pressure toward the astropause, through which mass flux is zero.
It is also evident from Fig. 2 that if the temperature falls below ≈4 × 10^{4} K, the flow again becomes supersonic.
Different values of the number density n or speed v or different cooling or heating functions change the characteristic lengths, but the principle estimates for the dynamics of the flow field remain the same.
In the stationary case, where no relative speed between the star and the ISM exists, one should take for the speed v that of the outwardmoving shock front. Moreover, the criterion given by Schwarz et al. (1972), namely that dlog (Λ(T))/dlog (T) > 2 to guarantee thermal stability, is generally not satisfied, and thus one expects clumping in between the AP and the BS and possibly in the ISM. The cooling lengths can be scaled by the speed and the density: L_{cool,s} does not change as long as v/n = const. This is not true for the heating length L_{heat,s}.
Fig. 2 Number density and speed are fixed at n = 11 cm^{3} and v = 80 km s^{1}. All scales are logarithmic. On the horizontal axis the temperature is displayed, while on the left vertical axis the characteristic length is shown in AU and in parsec on the right y axis. The black vertical line denotes the temperature where the ram pressure equals the thermal pressure. 
A similar result was also derived by van Marle & Keppens (2011). In the supersonic ISM the ram pressure dominates the cooling scales. It is expected that a magnetic field will increase the characteristic length scale due to its additional pressure and thus helps in stabilizing the shock structures (van Marle et al. 2014). The above estimate of the characteristic cooling length is crucial in determining the resolution required for numerical models of the largescale structure of astrospheres. From Fig. 2 it is clear that this should be higher than the shortest characteristic cooling lengths, which are on the order of >0.01 pc for the numbers given there. Moreover, the cooling of the surrounding ISM can also easily be inferred from the above estimate, and thus it can be approximately determined whether it is cooled during the computation.
4. Astrosphere model
The hydrodynamic model for the star λ Cephei is made in 3D; the third dimension is needed for future comparison with models including bipolar winds or magnetic fields. The boundary conditions are given in Table 2. While the stellar mass loss and the terminal speed can be determined by observations, the stellar wind temperature as well as the interstellar parameters are sophisticated guesses. The model solves the Euler equations for λ Cephei using the Cronos MHD code as described in Kissmann et al. (2008), Kleimann et al. (2009), and Wiengarten et al. (2013). The results are displayed in Fig. 3 for the proton number density of λ Cephei. In Fig. 3 the wiggles along the AP caused by the thermal instabilities can be clearly recognized. They are due to the cooling functions.
This model provides the underlying plasma structure needed as input for the transport equation discussed below. The model is solved on a spherical grid with a resolution of 0.005 pc in radial dimension, and 2^{◦} and 3^{◦} in ϑ and ϕ dimension, respectively. The large distance of the outer boundary is needed because during the evolution of the astrosphere it becomes much broader than shown in Fig. 3 and finally shrinks to the state shown here after ca. 170 kyr.
Fig. 3 Logarithmic number density for λ Cephei (including heating and cooling). The axes are given in pc, while the color bar is a logarithmic scale for the proton number density. 
Fig. 4 Continuumcorrected and slightly spatially filtered Hα image of the bow shock around λ Cephei. North is up, east is left, the size of the plotted field is 22′ × 22′. The imaging data were taken as part of the IPHAS (Drew et al. 2005) survey. A bright and welldefined bow shock nebula is visible toward the SE of the star. The roughly spherical Hα structure centered on λ Cephei is an outoffocus ghost of the bright star, the dark inner rings are the equivalent structure from the rfilter image used for continuum subtraction. In addition to these purely instrumental effects, several Hα bright features are visible close to the star. A more detailed discussion of the multiwavelength properties of the bow shock and the circumstellar environment of this runaway O6 star will be presented in an upcoming paper. 
5. Cosmic ray fluxes
To model the flux of GCRs, the Parker transport equation $\frac{\mathit{\partial f}}{\mathit{\partial t}}\mathrm{=}\mathrm{}{v}\mathrm{\xb7}\mathrm{\nabla}\mathit{f}\mathrm{+}\mathrm{\nabla}\mathrm{\xb7}\left({K}\mathrm{\xb7}\mathrm{\nabla}\mathit{f}\right)\mathrm{+}\frac{\mathrm{\nabla}\mathrm{\xb7}{v}}{\mathrm{3}}\frac{\mathit{\partial f}}{\mathit{\partial}\mathrm{ln}\mathit{P}}$(8)has been frequently used in the literature (see Potgieter et al. 2001, and reference therein). In this equation, f is the (nearly) isotropic GCR distribution function, v the bulk plasma flow, K the diffusion tensor (including, in a 3D geometry, separate components directed parallel and perpendicular to the mean magnetic field), and P is the particle rigidity. As an initial approach, spherical symmetry is assumed, so that the Parker equation reduces to $\frac{\mathit{\partial f}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathrm{}\mathit{v}\frac{\mathit{\partial f}}{\mathit{\partial r}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial}}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}{\mathit{\kappa}}_{\mathit{rr}}\frac{\mathit{\partial f}}{\mathit{\partial r}}\right)\mathrm{+}\frac{\mathit{P}}{\mathrm{3}{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial}}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{)}}\frac{\mathit{\partial f}}{\mathit{\partial}\mathrm{ln}\mathit{P}}\mathit{,}$(9)where r is radial distance and κ_{rr} is the effective radial diffusion coefficient. In this work, Eq. (9) is solved by transforming it into the equivalent set of SDEs
$\begin{array}{ccc}& & \mathrm{d}\mathit{r}\mathrm{=}\left[\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial}}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}{\mathit{\kappa}}_{\mathit{rr}}\mathrm{)}\mathrm{}\mathit{v}\right]\mathrm{d}\mathit{t}\mathrm{+}\sqrt{\mathrm{2}{\mathit{\kappa}}_{\mathit{rr}}}\mathrm{d}\mathit{W,}\\ & & \mathrm{d}\mathit{P}\mathrm{=}\left[\frac{\mathit{P}}{\mathrm{3}{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial}}{\mathit{\partial r}}\mathrm{(}{\mathit{r}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{)}}\right]\mathrm{d}\mathit{t,}\end{array}$which is then integrated numerically (Strauss et al. 2011), where the solution of the 1D equation along r can be found in Strauss et al. (2011).
These equations are coupled to the simulated HD geometry by reading in the modeled values of v and ∇·v (governing energy changes) directly from the HD simulations (along the stagnation line) and solving for the GCR flux using essentially a test particle approach. The magnetic field enters the computations via the diffusion coefficient κ_{rr} = κ_{rr}(B). In the 1D scenario, the geometry of B does not enter the computations, although for an azimuthal field (the case inside the TS), κ_{rr} ≈ κ_{⊥}, thus reducing to a diffusion coefficient perpendicular to the mean field. As a boundary condition for the GCR flux, a local interstellar spectrum is specified at the edge of the computational domain. The GCR differential intensity is related to the distribution function by j = P^{2}f.
Based on experience gained from modulation studies inside the heliosphere, κ_{rr} can be decomposed into a radial and energy dependence, κ_{rr} = κ_{1}(r)κ_{2}(P), where, for this study ${\mathit{\kappa}}_{\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ {\mathit{\kappa}}_{\mathrm{sw}}\mathrm{if}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{BS}}\\ {\mathit{\kappa}}_{\mathrm{ISM}}\mathrm{if}\mathit{r}\mathrm{\ge}{\mathit{r}}_{\mathrm{BS}}\end{array}\mathit{,}$(12)with r_{BS} the radial position of the BS and where κ_{ISM} is independent of position and κ_{sw} = κ_{0}B_{0}/B. κ_{0} is a normalization constant, usually specified near the inner boundary, and B_{0} is a constant included for dimensional consistency. Assuming that the astrospheric magnetic field is about 80 times higher than the heliospheric case (Naze 2014), a value of 0.027 kpc^{2} Myr^{1} is used, scaling up the heliospheric diffusion coefficient by about two orders of magnitude for this astrospheric case. The B^{1} radial dependence is based on the results of Engelbrecht & Burger (2013), but approximated in such a fashion because in situ observations are, of course, not available. The energy dependence of κ_{rr} is taken from Büsching & Potgieter (2008)${\mathit{\kappa}}_{\mathrm{2}}\mathrm{\left(}\mathit{P}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ {\left(\frac{\mathit{P}}{{\mathit{P}}_{\mathrm{0}}}\right)}^{\mathrm{+}\mathrm{0.6}}\mathrm{if}\mathit{P}\mathit{>}{\mathit{P}}_{\mathrm{0}}\\ {\left(\frac{\mathit{P}}{{\mathit{P}}_{\mathrm{0}}}\right)}^{0.48}\mathrm{if}\mathit{P}\mathrm{\le}{\mathit{P}}_{\mathrm{0}}\end{array}\mathit{,}$(13)with P_{0} = 4 GV. Instead of showing κ_{rr}, we follow the convention of showing the corresponding mean free path (mfp), $\mathit{\lambda}\mathrm{=}\frac{\mathrm{3}{\mathit{\kappa}}_{\mathit{rr}}}{{\mathit{v}}_{\mathrm{p}}}\mathit{,}$(14)where v_{p} is particle speed. The magnitude of κ_{rr} is changed in the next section to illustrate its effect on the resulting particle intensities.
Because the mfp used is, in general, not known, we study the behavior of three different values inside λ Cephei and three different values of the interstellar mfp.
5.1. Mean free path
The structure of the astrosphere can be recognized in the lower part of Fig. 5 by the respective jumps in the mfps (from the right, which is the inflow direction): the BS (1.1 pc), AP (0.93 pc) and TS (0.66 pc), represented by the vertical dotted lines. The TS is marked by a sharp decrease of the mfp, the AP by a change of the slope, and finally the BS by the sharp drop at the ISM side. The resulting differential flux (DI) for 1 GeV particles along the stagnation line is displayed in Fig. 6. In all cases the modulation, that is, the decrease in the differential intensity DI, starts far away in front of the BS, a feature that was discussed by Scherer et al. (2011), Herbst et al. (2012), and Strauss et al. (2013) for the case of the heliosphere.
Fig. 5 Modulation of 1 GeV particles in λ Cephei along the stagnation line for different stellar and interstellar mfp. The mfp are shown in the lower panel and the colors correspond to those in the upper panel. The magenta, orange, and green curves in the upper panel are identical when normalized to their respective values at infinity. The dotted lines show the position of the TS, AP, and BS. 
This outer astrospheric modulation depends on the ratio between the parallel and perpendicular diffusion coefficient and vanishes for high ratios. Thus, for an appropriate choice of the transport parameters, the modulation of GCRs starts in front of the astrosphere and then rapidly decreases at the astrospheric BS. It is also evident that in contrast to the heliosphere, the GCRs are modulated directly behind the BS: this is due to the cooling effects, which shrink the region between the AP and BS to 0.17 pc and are a barrierlike feature for the cosmic rays (Potgieter & Langner 2004).
Fig. 6 Modulation of particles with different energies. Upper panel: normalized differential particle intensities in part./ (m^{2}s sr MeV) for energies of 1 and 10 GeV are shown (in black and blue, respectively). The y axes are in a linear scale with different values for all three panels. Middle panel: DI for particles at 100 GeV and 1 TeV are displayed (black and blue). Lower panel: DIs for 10 and 100 TeV particles are shown. All DIs are normalized to their corresponding interstellar DI and are decrease towards the star. 
In the lower part of Fig. 5 we plot the mfp used in the modeling. The inner mfp (λ_{sw}) obtained from the model are then divided or multiplied by a factor two and three to demonstrate its influence. A small λ_{sw} strongly increases the modulation, and the flux of 1 GeV particles almost vanishes inside the AP. Increasing λ_{sw} by the same factor leads to a nearly vanishing modulation. Changing the interstellar λ_{ISM} from 50 pc to 100 pc and 500 pc does not change the DI remarkably. These variations show the effects expected from theory. Especially for short mfp inside the astrosphere the GCR spectra are efficiently cooled. The 20% modulation in front of the astrosphere is presumably caused by the sharp drop of the mfp between the AP and the BS.
From Fig. 5 we can also see that most of the modulation – up to a factor 5 or even more – occurs between the BS and the AP (depending on the transport parameters). Thus, because astrospheres are large volumes surrounding the star, the GCR spectrum can be cooled and, especially, smallscale anisotropies in the Galactic CRF can be observed far away from large astrospheres. This prediction will be improved by further modeling and studying the distribution of hot stars in the Galaxy.
At the BS, the mfp falls to quite low values, rises sharply at the AP, and finally jumps again at the TS. In the inner astrosheath the mfp increase continuously toward the shock because we assumed that the magnetic field is inversely proportional to the speed v as for the Parker spiral field in the heliosphere (Potgieter et al. 2001). Thus at the stagnation line, the speed must vanish at the AP because this is a tangential discontinuity with no mass flow through it. Inside the TS, the magnetic field increases with r in the direction to the star, and the mfp decreases. Because of the strong compression of the outer astrosheath, the mfp diminishes to an even lower value, which then forms a barrier for the cosmic ray transport.
This qualitative behavior is theoretically well understood and confirmed by observations in the heliosphere, except for the additional modulation in the outer astrosheath. The latter differs because of the cooling function. Thus, the cooling function has a direct influence on the cosmic ray transport.
5.2. Energy dependence
In Fig. 6 we study the behavior of the GCR flux at different energies. The 1 GeV particles are modulated by a factor of almost five, while for higher energies this factor becomes smaller, as expected. But even for the highest energies of 10 to 100 TeV, a small modulation of a few per mille (1 per mille = 0.1%) outside of the BS is visible.
In Fig. 6 we plot the DIs for different energies, from 1 GeV to 100 TeV in logarithmic steps of ten. In the upper panel the lower energies (form 1 to 10 GeV) are strongly modulated, while for the 10 GeV and 1 TeV particles the DI is only diminished by a few percent (middle panel). The lower panel shows the modulation of the 10 and 100 TeV particles, which is in the per mille range or even much lower for the 100 TeV particles, which are hardly affected. All energies are modulated ahead of the BS, and thus the modulation volume is larger than that of the astrosphere. The DI of the TeV particles can be increased by diminishing the mfp inside the astrosphere or can become less by increasing the mfp; this is not shown in Fig. 6. The mfp is based on that in the heliosphere, and the modulations displayed in Fig. 6 are based on it.
5.3. Observational evidence
Smallscale modulations of a few per mille in the TeV range are observed with, for example, the IceCube experiment (e.g., Abbasi et al. 2012,and references therein). In particular, their Fig. 5 shows a resolution of 3^{◦}, which cannot resolve the astrosphere of λ Cephei with angular diameter 0.2^{◦}. Nevertheless, the pixels show variations that may be attributed to local GCR fluctuations, provided the statistics in the pixels is large enough and the astrospheres are close enough. Because runaway O/B stars or OB associations are quite frequent in the Galaxy (Huthoff & Kaper 2002), their local disturbance of the GCR spectrum can explain smallscale variations in the GCR flux. λ Cephei was only used as an example, but the above holds true in general for all large astrospheres.
Our simulations may help to understand these variations. The variation in the DIs can be increased or decreased by choosing a smaller or larger mfp inside the astrosphere, therefore a better estimate of their magnetic fields is needed, which can help to understand their turbulence levels.
For the higher energies a modulation outside the astrosphere can also be observed, where the “astrosphere of influence” is roughly twice the BS distance for bulletshaped astrospheres like that of λ Cephei. This holds true in the direction of the inflowing ISM, while it is more complicated in other directions. These details will be explored in future. We studied here only protons, but because the transport coefficients depend on rigidity, the behavior of other species can roughly be estimated. For helium, for instance, a similar behavior is expected as for the protons shifted by a factor of two in Fig. 6. At these high energies the rigidity has roughly the same value as the particle energy. This means that multiplying the rigidity by a factor of 2 is the same as doing it in the particle energy. Therefore, helium is slightly less modulated when passing through astrospheres than the protons.
If there were a few large astrospheres (or stellar wind bubbles) in the direction toward an observer, the GCR spectrum would be slightly cooler than in other directions. This can contribute to the understanding of the smallscale anisotropies present in the IceCube data (Abbasi et al. 2012).
The parameter set for the supersonic wind and ISM are chosen to be close to the observed values, but for modeling purposes we worked with rounded values. In addition, for the transport model we took the turbulence levels based on that of the heliosphere. Both of these aspects need further attention, but the principal effect remains: large stellar astrospheres can affect the local interstellar cosmic ray spectra.
6. Conclusion
Based on our models, we studied the transport of GCRs in an astrosphere. We have shown that even ahead of an astrosphere, where the mfp is still undisturbed, a cooling of the GCR spectra occurs because particles are trapped by scattering into the astrosphere, in which they can be effectively cooled. The “astrosphere of influence” around the modeled astrospheres is roughly twice its hydrodynamic dimension. Thus, stellar wind bubbles (astrospheres) can cool the Galactic spectrum, and small anisotropies are expected in the direction to an observer.
As a result of its effect in compressing the outer astrosheath, the cooling function directly influences the GCR modulation in this region. Additionally, the cooling and heating functions require high resolutions for global models because of their characteristic length scales. They also give a rough estimate whether the surrounding ISM is also affected during the calculations.
The modulation in astrospheres affects particles up to 100 TeV. This can help to understand the anisotropies on small angular scales, for example, by the IceCube experiment, among others. The angular extent of λ Cephei’s astrosphere is too small to be resolved by these experiments because of its large distance (≈650 pc). Nevertheless, nearby hot stars or stellar associations can have an larger angular extent and may possibly be observed. In our simulations we saw an effective cooling of the GCR spectrum. Thus, the conclusion is that smallscale cosmic ray anisotropies may be explained by the modulation in such huge cavities.
Because the model is fully 3D, the modulation or other parameters along a line of sight toward Earth can in principle be theoretically determined. Here we demonstrated that large astrospheres can modulate the Galactic cosmic ray spectrum quite significantly, and no homogeneous spectrum over all directions can be expected. To calculate the modulation of GCRs along a line of sight or at Earth requires more sophisticated methods than discussed here, but this is being developed.
Acknowledgments
K.S., H.F., J.K., and T.W. are grateful to the Deutsche Forschungsgemeinschaft, DFG funding the projects FI706/151 and SCHE334/101. D.B. is supported by the DFG Research Unit FOR 1254. This work is also based on the research supported in part by the South African NRF. Any opinion, finding, and conclusion or recommendation expressed in this material is that of the authors, and the NRF does not accept any liability in this regard. K.S., H.F. and R.D.S. appreciate discussions at the team meeting “Heliosheath Processes and Structure of the Heliopause: Modeling Energetic Particles, Cosmic Rays, and Magnetic Fields” supported by the International Space Science Institute in Bern, Switzerland. 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
 Abbasi, R.,Abdou, Y.,AbuZayyad, T., et al. 2012, ApJ, 746, 33 [Google Scholar]
 Arthur, S. J. 2012, MNRAS, 421, 1283 [NASA ADS] [CrossRef] [Google Scholar]
 Baranov, V. B.,Krasnobaev, K. V., &Kulikovskii, A. G. 1971, Sov. Phys. Doklady, 15, 791 [Google Scholar]
 Büsching, I., &Potgieter, M. S. 2008, Adv. Space Res., 42, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Cox, N. L. J.,Kerschbaum, F., van Marle, A. J., et al. 2012, A&A, 543, C1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Decin, L.,Cox, N. L. J.,Royer, P., et al. 2012, A&A, 548, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drew, J. E.,Greimel, R.,Irwin, M. J., et al. 2005, MNRAS, 362, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Engelbrecht, N. E., &Burger, R. A. 2013, ApJ, 772, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Henrichs, H., & Sudnik, N. P. 2013, in Massive Stars: From alpha to Omega, 71 [Google Scholar]
 Herbst, K.,Heber, B.,Kopp, A.,Sternal, O., &Steinhilber, F. 2012, ApJ, 761, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Huthoff, F., &Kaper, L. 2002, A&A, 383, 999 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kissmann, R.,Kleimann, J.,Fichtner, H., &Grauer, R. 2008, MNRAS, 391, 1577 [Google Scholar]
 Kleimann, J.,Kopp, A.,Fichtner, H., &Grauer, R. 2009, Annales Geophysicae, 27, 989 [Google Scholar]
 Kosiński, R., &Hanasz, M. 2006, MNRAS, 368, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Mellema, G., &Lundqvist, P. 2002, A&A, 394, 901 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Naze, Y. 2014, in Putting A Stars into Context: Evolution, Environment, and Related Stars, Proc. Int. Conf. held on June 3–7, 2013 at Moscow M.V. Lomonosov State University, Russia, eds. G. Mathys, E. Griffin, O. Kochukhov, R. Monier, & G. Wahlgren (Moscow: Publishing house, Pero), 340 [Google Scholar]
 Pauls, H. L.,Zank, G. P., &Williams, L. L. 1995, J. Geophys. Res., 100, 21595 [NASA ADS] [CrossRef] [Google Scholar]
 Potgieter, M. S., &Langner, U. W. 2004, ApJ, 602, 993 [NASA ADS] [CrossRef] [Google Scholar]
 Potgieter, M. S.,Burger, R. A., &Ferreira, S. E. S. 2001, Space Sci. Rev., 97, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Reitberger, K.,Kissmann, R.,Reimer, A., &Reimer, O. 2014, ApJ, 789, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, R. J.,Haffner, L. M., &Tufte, S. L. 1999, ApJ, 525, L21 [Google Scholar]
 Rosner, R.,Tucker, W. H., &Vaiana, G. S. 1978, ApJ, 220, 643 [NASA ADS] [CrossRef] [Google Scholar]
 Scherer, K.,Fichtner, H.,Strauss, R. D., et al. 2011, ApJ, 735, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Scherer, K.,Fichtner, H.,Fahr, H.J.,Bzowski, M., &Ferreira, S. 2014, A&A, 563, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schure, K. M.,Kosenko, D.,Kaastra, J. S.,Keppens, R., &Vink, J. 2009, A&A, 508, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, J.,McCray, R., &Stein, R. F. 1972, ApJ, 175, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Siewert, M.,Pohl, M., &Schlickeiser, R. 2004, A&A, 425, 405 [Google Scholar]
 Strauss, R. D.,Potgieter, M. S.,Kopp, A., &Büsching, I. 2011, J. Geophys. Res. Space Phys., 116, 12105 [Google Scholar]
 Strauss, R. D.,Potgieter, M. S.,Ferreira, S. E. S.,Fichtner, H., &Scherer, K. 2013, ApJ, 765, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Sutherland, R. S., &Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, A. R.,Gibson, S. J.,Peracaula, M., et al. 2003, AJ, 125, 3145 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D. 2009, ApJS, 181, 391 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Marle, A. J., &Keppens, R. 2011, Computers and Fluids, 42, 44 [Google Scholar]
 van Marle, A. J.,Decin, L., &Meliani, Z. 2014, A&A, 561, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wiengarten, T.,Kleimann, J.,Fichtner, H., et al. 2013, J. Geophys. Res. Space Phys., 118, 29 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Some characteristic numbers using the cooling function from Siewert et al. (2004).
All Figures
Fig. 1 Largescale structure of an astrosphere. For details see text. 

In the text 
Fig. 2 Number density and speed are fixed at n = 11 cm^{3} and v = 80 km s^{1}. All scales are logarithmic. On the horizontal axis the temperature is displayed, while on the left vertical axis the characteristic length is shown in AU and in parsec on the right y axis. The black vertical line denotes the temperature where the ram pressure equals the thermal pressure. 

In the text 
Fig. 3 Logarithmic number density for λ Cephei (including heating and cooling). The axes are given in pc, while the color bar is a logarithmic scale for the proton number density. 

In the text 
Fig. 4 Continuumcorrected and slightly spatially filtered Hα image of the bow shock around λ Cephei. North is up, east is left, the size of the plotted field is 22′ × 22′. The imaging data were taken as part of the IPHAS (Drew et al. 2005) survey. A bright and welldefined bow shock nebula is visible toward the SE of the star. The roughly spherical Hα structure centered on λ Cephei is an outoffocus ghost of the bright star, the dark inner rings are the equivalent structure from the rfilter image used for continuum subtraction. In addition to these purely instrumental effects, several Hα bright features are visible close to the star. A more detailed discussion of the multiwavelength properties of the bow shock and the circumstellar environment of this runaway O6 star will be presented in an upcoming paper. 

In the text 
Fig. 5 Modulation of 1 GeV particles in λ Cephei along the stagnation line for different stellar and interstellar mfp. The mfp are shown in the lower panel and the colors correspond to those in the upper panel. The magenta, orange, and green curves in the upper panel are identical when normalized to their respective values at infinity. The dotted lines show the position of the TS, AP, and BS. 

In the text 
Fig. 6 Modulation of particles with different energies. Upper panel: normalized differential particle intensities in part./ (m^{2}s sr MeV) for energies of 1 and 10 GeV are shown (in black and blue, respectively). The y axes are in a linear scale with different values for all three panels. Middle panel: DI for particles at 100 GeV and 1 TeV are displayed (black and blue). Lower panel: DIs for 10 and 100 TeV particles are shown. All DIs are normalized to their corresponding interstellar DI and are decrease towards the star. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.