© ESO, 2012

1. Introduction

For decades, it has been understood that the mass lost by red giants and supergiants dominates the mass return to the Galaxy (Maeder 1992), yet the actual mechanisms of such stellar mass loss have remained a mystery. The generally accepted idea on the mass-loss mechanism for asymptotic giant branch (AGB) stars is based on pulsations and radiation pressure on newly formed dust grains. However, oxygen-rich AGB stars suffer from the so-called “acceleration deficit” dilemma, which states that mass-loss rates due to the formation of silicate dust alone are orders of magnitude less than observed ones (Woitke 2006b). The formation of both carbon and silicate grains (Höfner & Andersen 2007) or micron-sized Fe-free silicates (Höfner 2008; Norris et al. 2012) have been proposed as possible solutions to this dilemma. The proposed models for the mass loss of AGB stars are very unlikely to be applicable to red supergiants (RSGs) since they are irregular variables with small amplitudes.

Processes linked to convection, chromospheric activity, or rotation might play an important role as triggers for the mass loss in RSGs. Josselin & Plez (2007) proposes that turbulent pressure generated by convective motions, combined with radiative pressure on molecular lines, might initiate mass loss, but Alfvén waves generated by a magnetic field might also contribute (Hartmann & Avrett 1984). The last few years have provided more and more evidence that, at least for a fraction of evolved stars, the mass is not lost in a homogeneous way, but irregular and clumpy structures prevail (e.g. Weigelt et al. 2002; Kervella et al. 2009; Decin et al. 2011). This might have serious consequences for the total amount of mass lost during these late evolutionary stages. Solving the riddle of the mass loss of these (higher mass) targets is important in the framework of stellar evolution as a key to estimating when a supernova explosion may occur.

When these evolved stars with their surrounding winds move through the interstellar medium (ISM), a wind-ISM interface structure is expected in the form of a bow shock, a kind of cometary like structure pointing in the direction of motion. Photometric observations performed with the IRAS telescope (Noriega-Crespo et al. 1997) and later confirmed with the Spitzer Space Telescope (Ueta et al. 2006, 2008) and the Galaxy Evolution Explorer satellite (GALEX Martin et al. 2007; Sahai & Chronopoulos 2010) indeed showed the presence of bow shock structures around evolved stars. However, these early observations of the far-IR and UV-bright emission structures lacked any spectral resolution, and since they only had poor spatial resolution, allowed only basic morphological studies.

Thanks to the capabilities of the instruments on board the Herschel Space Observatory (Pilbratt et al. 2010), one is now able to image the circumstellar envelopes (CSEs) and the CSM-ISM interaction regions with unprecedented detail. In the MESS GTKP (Groenewegen et al. 2011) several known and new interaction objects have been imaged with the PACS and SPIRE instruments (e.g. CW Leo, Ladjal et al. 2010; Decin et al. 2011; Cox et al. 2012), among which several supergiants. Betelgeuse (α Ori, M2 Iab) is the closest supergiant, and it serves as the prototype of the red oxygen-rich supergiants. It has been an important astrophysical laboratory for many decades and its envelope has been studied with a variety of observational techniques. Betelgeuse is, therefore, a prime target for Herschel to study the still unknown mechanisms that lead to mass loss and eventually result in an energetic bow shock when the wind collides with the ISM.

1.1. Betelgeuse

Many detailed studies have already described some characteristics of the complex atmosphere, chromosphere, and dusty envelope of Betelgeuse. In this section, a summary is given of those properties that are important for this study.

Betelgeuse is a low-amplitude variable, whose variation was first described in 1840 by Sir John Herschel (Herschel 1840). Schwarzschild (1975) attributed these fluctuations to a changing granulations pattern formed by a few giant convection cells covering the surface of the star. Recently, a magnetic field of about 1 Gauss has been detected using the Zeeman effect (Aurière et al. 2010). It is very faint, but may play an important role in shaping the convection of Betelgeuse.

Using NRAO’s Very Large Array (VLA), combined with Hipparcos data, Harper et al. (2008) determined the parallax of Betelgeuse to be 5.07 ± 1.10 mas (or a distance of 197 ± 45 pc). The most likely star-formation scenario for Betelgeuse is that it is a runaway star from the Orion OB1 assocation and was originally a member of a high-mass multiple system within Ori OB1a (Harper et al. 2008). The galactic space motion in the local standard of rest (LSR) is U′ =  −12.7, V′ =  +0.5, W′ =  +28.2 km s^-1. Betelgeuse was probably formed about 10 to 12 million years ago from the molecular clouds observed in Orion, but has evolved rapidly owing to its high mass. No consensus has yet emerged regarding the star’s mass, but most studies report values between 10 M_⊙ and 20 M_⊙ (e.g., Dolan et al. 2008; Neilson et al. 2011).

Ultraviolet and Hα images prove there is a chromosphere with a spatial extent of  ~4 R _⋆  (Altenhoff et al. 1979; Goldberg 1981; Hebden et al. 1987), where the hot plasma has a temperature of around 5000 K (Gilliland & Dupree 1996). However, observations made with the VLA by Lim et al. (1998) suggest that much cooler (~1000−3000 K) gas co-exists at the same radial distances from the central target, and must be much more abundant because it dominates the radio emission. In addition, narrow-slit spectroscopy of Betelgeuse at 11.15 μm with the ISI by Verhoelst et al. (2006) reveals that silicate dust forms at distances beyond  ~20 R _⋆ , but that Al₂O₃ may form as close as  ~2 R _⋆ . This apparent dichotomy is solved when adopting a scenario in which a few inhomogeneously distributed large convective cells are responsible for lifting the cooler photospheric gas into the atmospheres and for producing shock waves, which could heat localized, but spatially distributed regions of the atmosphere to chromospheric temperatures (Lim et al. 1998).

Observations taken over the past few years have proven that the regions just above the stellar surface show some clear deviations from spherical symmetry. Kervella et al. (2009) obtained adaptive optics images of Betelgeuse in the 1.04 − 2.17 μm range in ten narrow-band filters. These images show that the CSE of Betelgeuse has a complex and irregular structure, with in particular a bright “plume” extending in the southwestern quadrant up to a radius of at least six times the photosphere. They proposed that the “plume” might be linked to either the presence of a convective hot spot on the photosphere or to polar mass loss possibly due to rotation of the star. Using AMBER on the VLTI, Ohnaka et al. (2009) found from CO lines that the atmosphere has an inhomogeneous velocity field, which could be explained by the gas moving outward and inward in a patchy pattern. This gives support to the idea of convective motion in the upper atmosphere or intermittent mass ejections in clumps or arcs. Recently, Kervella et al. (2011) observed Betelgeuse with VLT/VISIR in the N and Q band. The images showed a bright, extended and complex CSE within  ~2.5′′ of the central star with the appearance of a partial circular shell between 0.5 and 1.0′′, several knots, and filamentary structures. The ring-like structure might be related to the dust condensation zone.

Although some dust species, such as silicates and corundum, have been identified in the CSE around Betelgeuse, the envelope shows a remarkable deficiency in dust (Skinner & Whitmore 1987; Verhoelst et al. 2006). In addition, the molecular content is unusually low. There is evidence that only a fraction of carbon is associated into CO (Huggins 1987). The extended, warm chromosphere is clearly unfavourable for the survival of molecules and the formation of dust species. Detailed studies by Glassgold & Huggins (1986) show that the most common form of all species are neutral atoms and first ions.

A parsec-size shell around Betelgeuse was detected with IRAS by Noriega-Crespo et al. (1997), and later on with AKARI by Ueta et al. (2008) (see Appendix A). The shell is asymmetric with a mean radius of 6′ and the outermost structure at  ~7′. Noriega-Crespo et al. (1997) proposed that the shell is confined by the ram pressure of the ISM, hence represents the bow shock around Betelgeuse. To the north-east of Betelgeuse, a ridge or filament was detected, which was said to belong to the surrounding ISM. Neither IRAS nor AKARI could clearly resolve the bow shock and filament structure. The derived dust temperatures in the bow shock and filament vary between  ~10 K and 40 K (Ueta et al. 2008).

Recently, Le Bertre et al. (2012) has detected a detached H i gas shell of  ~2′ in radius surrounding Betelgeuse, and reported the detection of atomic hydrogen associated with the far-infrared arc located 6′ north-east of Betelgeuse.

In this paper, we focus on the morphological appearance of the inner envelope and bow shock surrounding Betelgeuse. The thermodynamical structure and chemical content of the envelope and bow shock will be discussed in Decin et al. (in prep., from here on called Paper II). In Sect. 2, the different observations with their respective calibrations are described. The observational properties of the inner envelope and bow shock structure are analysed in Sect. 3. With the aid of hydrodynamical simulations, the origin of the multiple arcs is discussed in Sect. 4. The origin of the linear bar is described in Sect. 5. The conclusions are summarized in Sect. 6.

2. Observations

In this section, an overview is given of the different observations toward Betelgeuse and its extended envelope, the data-reduction for each data-set is described, and some first results are deduced from the data.

2.1. Herschel PACS and SPIRE images

Infrared images were obtained using the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) and the Spectral and Photometric Imaging Receiver (Griffin et al. 2010, SPIRE) on board the Herschel satellite (Pilbratt et al. 2010), and they are part of the MESS guaranteed-time key programme (Groenewegen et al. 2011). For both instruments, the “scan map” observing mode was used with medium (20′′/s) scan speed in the PACS 70, 100, and 160 μm filter settings and 30′′/s for the SPIRE 250, 350, and 500 μm filters. To create a uniform coverage, avoid striping artefacts, and increase redundancy, two observations at orthogonal scan directions were concatenated. The PACS images were obtained on September 13, 2010 (OBSIDS 1342204435 and 1342204436) and on March 29, 2012 (OBSIDS 1342242656 and 1342242657), the SPIRE images on March 11, 2010 (OBSID 1342192099).

The data reduction was performed using the Herschel Interactive Processing Environment (HIPE, see Groenewegen et al. 2011, for more details). To remove low-frequency noise, which causes additive brightness drifts, the scanamorphos routine (Roussel 2012) was applied. Deconvolution of the data was done using the method as described in Ottensamer et al. (2011).

The PACS instrument offers a pixel size of 3.2′′ in the 70 μm and 100 μm bands and of 6.4′′in the 160 μm band, but the images presented in this paper are oversampled by a factor of 3.2, resulting in a sampling of 1′′ and 2′′ per pixel, respectively. The PACS FWHM at 70, 100, and 160 μm is 5.7′′, 6.7′′, and 11.4′′, respectively. The pixel size of the three SPIRE bands is 6′′ at 250 μm, 10′′ at 350 μm, and 14′′ at 500 μm, with an FWHM of the beam size being 18.1′′, 25.2′′, and 36.6′′, respectively.

The individual scan maps are plotted in Fig. 1, and a false-colour image displaying all the PACS data is shown in Fig. 2. Remarkably, a multiple arc-like structure is detected at  ~6−7′ away from the central star, together with an enigmatic linear bar-like structure at  ~9′. A closer zoom into the arc-like structure is displayed in Fig. 3. The inner radius of the closest arc is at 280′′ (as measured from the central target). The other most prominent arcs have an inner radius at 280′′, 310′′, 350′′ and 375′′. The brightest parts of the arcs are reasonably well fitted by concentric ellipses (see Fig. B.1 in the Appendix B). The width of the arcs is typically  ~20′′. Clearly, substructures of the size of several arcseconds, probably related to density differences, are seen in the arcs. The arc-like structure is probably related to the CSM-ISM interaction phase, when the circumstellar material collides at a wind speed of  ~15 km s^-1 with the ISM, creating a bow shock (see Sect. 4 for more details). The inner ridge of the linear bar is at 535′′ from the central star. The length is  ~1600′′. The width is  ~30′′ in the northern part; moving down towards the south-east, the bar widens and splits into two parts. The angle between the direction of space velocity (drawn in Fig. 2) and the linear bar is  ~80°.

In the whole sample of 78 AGB stars and red supergiants surveyed in the framework of the MESS GTKP, Betelgeuse is the most spectacular source, because it is the only one showing this multiple-arc structure and the presence of a linear bar. The linear bar and the full arc were already detected with IRAS (Noriega-Crespo et al. 1997) and AKARI (Ueta et al. 2008), but their individual structures could not be resolved (see Appendix A).

To emphasize the shell morphology in the inner envelope, we removed the extended envelope halo by subtracting a smooth, azimuthally averaged profile represented by a power law r ^− α for r < 240′′ (see Decin et al. 2011). The method has been tested in detail to determine if artefacts from the point-spread function (PSF) still affect the data (Royer et al., in prep.). Only the region within 15−20′′ from the central target should be interpreted with care owing to the complex PSF and the bright contribution from the central target at these infrared wavelengths. The resulting image is shown in Fig. 4. The inner envelope structure is clearly non-homogeneous and testifies to a turbulent mass-loss history. A detailed discussion of the inner envelope structure is presented in Sect. 3.1.

2.2. WISE data

The complex environment around Betelgeuse has also been imaged by the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010). WISE has four photometric bands centred at 3.4, 4.6, 12, and 22 μm, with an angular resolution of 6.1′′, 6.4′′, 6.5′′, and 12′′, respectively. The 5σ point source sensitivity in unconfused regions on the ecliptic is better than 0.08, 0.11, 1, and 6 mJy, respectively. The WISE All-Sky Data were released on March 14, 2012. The WISE data images are shown in Fig. 5. The extreme brightness of Betelgeuse at these infrared wavelengths makes it very difficult to find any nearby associated emission structures. The linear bar is visible only at 12 and 22 μm, with a (background-subtracted) surface brightness of  ~1.2 MJy/sr and  ~1.9 MJy/sr, respectively.

2.3. GALEX observations

Betelgeuse was also observed with GALEX on November 29, 2008. The pipeline-calibrated FUV and NUV images were retrieved from the GALEX archive. The bandpass is 1344 − 1786 Å for an angular resolution of 4.5′′ for the FUV filter, and 1771−2831 Å for a 6′′ angular resolution for the NUV filter (Morrissey et al. 2005). The pixel size is 1.5′′ × 1.5′′. The integration time for Betelgeuse was 48 916 s. Some very faint extra extinction accompanying the outermost arc around Betelgeuse was detected in the FUV filter (see Fig. 6). The FUV flux of the faint extra extinction is  ~7% lower than in the adjacent pixels. The NUV filter is unfortunately satured by a ghost of the central star.

Until now, bright UV emission in the wind-ISM interaction region has only been detected in two AGB stars: o Cet (Martin et al. 2007) and CW Leo (Sahai & Chronopoulos 2010); α Ori is the first target for which extra UV extinction in the wind-ISM interaction zone is reported. No UV emission or extinction is seen in the region of the linear bar. The origin of the FUV/NUV emission of the bow shock around o Cet has been attributed to H₂ molecules in the cold gas that are collisionally excited by hot electrons from the post-shock gas (Martin et al. 2007). The origin of the FUV/NUV emission in the interaction zone between the ISM and the stellar wind of CW Leo has not been modelled in detail, although it was suggested by Sahai & Chronopoulos (2010) that also here collisionally excited H₂ emission might be the origin. The FUV/NUV emission around CW Leo peaks at slightly larger distances from the central target than the PACS/SPIRE infrared flux excess (see Fig. 2 in Ladjal et al. 2010). We note that for the planetary nebula NGC 6720, dust and H₂ are co-spatial, and it is argued that H₂ has been formed on grain surfaces (van Hoof et al. 2010). If the FUV extinction toward the bow shock around Betelgeuse is the FUV counterpart of the outer arcs as detected with PACS, this suggests that background UV radiation is absorbed by dust in the bow shock.

2.4. GALFA H i 21 cm observations

H i 21 cm observations in the vicinity of Betelgeuse were extracted from the results of the Galactic Arecibo L-Band Feed Array H i (GALFA-H i) survey (Peek et al. 2011). GALFA-H i is a survey of the Galactic interstellar medium with a spatial resolution of  ~4′, covering a large area (13 000 deg²) with a high spectral resolution (0.18 km s^-1) from  −700 < v_LSR <  + 700 km s^-1 in the 21-cm line hyperfine transition of neutral hydrogen conducted at the Arecibo Observatory. ALFA is a seven-beam feed array, with the beams arranged on the sky in a hexagon. Each beam is slightly elliptical along the zenith-angle direction, approximately 3.3′  ×  3.8′FWHM, with a gain of  ~11 K/Jy for the central beam and  ~8.5 K/Jy for the six other beams. Details on the data reduction are given in Peek et al. (2011). Typical noise levels are 80 mK rms in an integrated 1 km s^-1 channel, corresponding to an integration time of 30 s.

Several selected data cubes are displayed in Fig. 7. At the CO v_LSR velocity of  ~4 km s^-1 (De Beck et al. 2010), H i emission is clearly detected. In Fig. 8, we present the H I spectrum obtained on the position of Betelgeuse. The spectrum is in good agreement with the Nançay RadioTelescope (NRT) and Leiden/Argentina/Bonn (LAB) data presented in Fig. 3 by Le Bertre et al. (2012). The emission peaks at 6.992 km s^-1, and no distinctive feature is seen in the velocity range of the circumstellar CO emission. It illustrates that in the direction of Betelgeuse, the H i emission is dominated by galactic emission due to interstellar atomic hydrogen along the line of sight (Le Bertre et al. 2012). The H i spectrum was also inspected at different offset positions. Taking the central star in the “on-position” and the “off-positions” at different offsets gave the H i spectrum of Betelgeuse (see Fig. 9). Although affected by interstellar confusion, emission at the CO v_LSR is clearly detected in the spectra with the offsets taken in the east-west direction. Figure 10 shows the H i spectrum at different offset positions for which the off-position was taken in a region outside the envelope, arcs or linear bar (i.e. at 11′ to the west and 10′ to the north of the central target). The interstellar confusion is seen for instance in the emission feature at  −3 km s^-1, and the absorption at +12 km s^-1. The Betelgeuse H i line profile has an FWHM of  ~3.5 km s^-1. The integrated line intensity is  ~2−7 Jy km s^-1. Assuming a smooth outflow and using the standard relation, M_H   i = 2.36 × 10^-7D^2∫S_ν   dv, with M_H   i in M_⊙ , the distance D in pc, and ^∫S_ν   dv in Jy km s^-1, this translates to an envelope mass of  ~0.02–0.07 M_⊙ in atomic hydrogen at 197 pc. As discussed already, the density in the inner envelope is not homogeneous, for which reason this mass value should be considered as very approximate.

Significant emission is also seen at other velocities, of which a few examples are shown in Fig. 7. Using VLA data, Le Bertre et al. (2012) note a seemingly spatial association between the H i emission integrated over the range  − 1.5 to 8.9 km s^-1 and the IRAS 60 μm image (see their Fig. 9). This coincidence is not seen between the GALFA-H i and Herschel 70 μm data, probably owing to the low spatial resolution of the GALFA-H i data (see upper panel in Fig. 11).

The bottom panels in Fig. 7, showing the emission at a velocity of 18 and 21.17 km s^-1, notably display an alignment of the H i emission with the bar and arc-like structures detected in the Herschel PACS images (see also bottom panel in Fig. 11). This H i emission probably has a galactic ISM origin, and might be part of cold, low-density, atomic gas structures (see also Sect. 5). To gain further insight into the origin of the bar-like structure seen in the PACS observations, we compared the galactic H i emission at different velocities in a larger region covering a right ascension from 05^h39^m56^s to 05^h59^m33^s and a declination from 07${\mathrm{◦}}^{}$21′10.23′′ to 9${\mathrm{◦}}^{}$42′30.47′′. From these data, it is clear that the ISM contains a lot of small-scale structure and that most of the emission seen in the direction of Betelgeuse is not related in any way to the mass loss of Betelgeuse. Because it is situated beyond the Local Bubble, Betelgeuse probably lies near the edge of a cavity of low neutral hydrogen density and close to a steep gradient in ISM column density between l = 180° − 200° (see also Harper et al. 2008). The ISM small-scale structure in the direction of Betelgeuse has also been detected by Knapp & Bowers (1988) via the observation of discrete molecular CO clouds, probably associated with the λ Ori molecular cloud complex.

3. Observational properties of the infrared (extended) envelope around Betelgeuse

In this section, some observational properties as deduced from the infrared data for the inner envelope structure, and the arcs and linear bar are described in Sects. 3.1 and 3.2, respectively.

3.1. Inner envelope structure

Figure 4 shows that the inner envelope structure of Betelgeuse is clearly non-homogeneous. The detection of non-spherical structures has already been reported on a sub-arcsecond-scale (e.g. Lim et al. 1998; Kervella et al. 2009, 2011) (see Sect. 1.1). The Herschel images show the first evidence of a high degree of non-homogeneity of the material lost by the star beyond 15′′, which even persists until the material collides with the ISM. As described in Sect. 2.1 and shown in Fig. 4, very pronounced asymmetries are visible within 110′′ from the central star, although some weaker flux enhancements are visible until  ~300′′ away. The typical angular extent for these arcs in the inner envelope ranges from  ~10° to 90°. Considering that asymmetries are seen on a subarcsecond-scale close to the stellar atmosphere and in the dust formation region, the Herschel images suggest that the arc-like structures seen in the free-flowing envelope of Betelgeuse (i.e. before the material collides with the ISM) are the relics of a non-homogeneous mass-loss process.

The radial distance of the most pronounced asymmetries (within 2′) agrees with the size of a ring of H i emission recently detected by Le Bertre et al. (2012), although no counterpart to the H i emission plateau further away from the star in the south-western region is seen in the Herschel images. The Herschel and H i observations seem to suggest that the mean gas and dust density (averaged over 360°) in the inner envelope drastically changed at  ~2′ (or 32 000 yrs ago for an expansion velocity of 3.5 km s^-1). This might be due to, e.g., a change in mass-loss rate or the creation of an inner bow shock inside a fragmented, filamentary halo (see Sect. 4.1).

The angular extent of the arc-like structures in the inner envelope of Betelgeuse is considerably smaller than for the AGB star CW Leo, the only other target in the MESS sample whose inner envelope has already been studied in detail (Decin et al. 2011). In CW Leo, almost spherical, ring-like, non-concentric dust arcs where detected until 320′′, having an angular extent between  ~40° and  ~200°. The shells of CW Leo have a typical width of 5′′–8′′, and the shell separation varies in the range of  ~10′′–35′′, corresponding to  ~500−1700 yrs. This (again) suggests that the mass-loss mechanism of the AGB star CW Leo, based on pulsations and radiation pressure on dust grains that formed non-isotropically, does not hold for Betelgeuse, which is an irregular variable with small amplitude variations. As described in Sect. 1.1, the mass-loss process might by convection-induced, which would yield locally strong variations in gas and dust density.

Assuming that dust is the main contributor to the Herschel/PACS images (see Paper II), one can derive the dust temperature in the arcs from a modified blackbody of the form B_ν(T)·λ ^− β, as expected for a grain emissivity Q_abs ~ λ ^− β with the emissivity index β ranging between one¹ (typical of, e.g., layered amorphous silicate grains, Knapp et al. 1993) and two (typical for, e.g., crystalline silicate grains, Tielens & Allamandola 1987; Mennella et al. 1998). When β is equal to 2.0, the dust temperature for the clumps beyond 60′′ ranges from  ~25 to 65 K. For β equal to 1.0, the dust temperature increases to values between 35 K and 140 K (see Fig. 13 in Sect. 4.1).

The Herschel images prove that clumpy structures are prevalent over the full envelope and might eventually have an impact on the shape of the bow shock structure.

3.2. Multiple arcs and linear bar

3.2.1. Orientation of the bow shock

Using the Herschel PACS images, Cox et al. (2012) derived a de-projected stand-off distance for the outer arcs of 4.98′ and a position angle of 47.7°. Based on the distance of 131 pc and proper motion values of van Leeuwen (2007), a space motion vector inclination i_LSR of 8.9° for a v_rad,LSR of 3.4 km s^-1 was deduced. From the AKARI data, Ueta et al. (2008) deduced a de-projected stand-off distance for the outer arcs of 4.8′  ±  0.1′, a position angle of 55°  ±  2° and an inclination angle for the outer arcs of 56°  ±  4°, where the bow shock cone is oriented to the plane of the sky. From the difference between the galactic space-velocity components of the bow shock cone and the heliocentric galactic space-velocity components of Betelgeuse, Ueta et al. (2008) find that the ISM around Betelgeuse flows at  ~11 km s^-1 into the position angle of  ~ 95° out of the plane of the sky (towards us). Differences between the Herschel (Cox et al. 2012) and AKARI (Ueta et al. 2008) results are mainly because of using different values for the radial velocity, proper motion, and distance. Ueta et al. (2008) and Cox et al. (2012) have deduced that the peculiar velocity of the star with respect to the ISM is between 24 and 33 km s^-1 for an interstellar hydrogen nucleus density between 1.5 and 1.9 cm^-3.

3.2.2. Instabilities in the arcs

The arc-like structure around Betelgeuse does not show clear evidence of large-scale instabilities, created by e.g. Kelvin-Helmholtz (KH), Rayleigh-Taylor (RT), Vishniac non-linear thin shell (NTSI Vishniac 1994), or transverse acceleration (TAI Dgani et al. 1996) instabilities. This is in contrast to other targets in the MESS sample as, e.g., R Leo for which the Herschel image shows clear signatures of RT instabilities that are slightly deformed by to the action of the KH instabilities (see Fig. C.1 in the Appendix C). The upper limit on the maximum length for instabilities in the arcs around Betelgeuse as traced from the Herschel images is  ~30′′. Although suggested by Ueta et al. (2008), vortex shedding (Wareing et al. 2007b) downstream in the tail of the bow shock is not seen in the Herschel images. With a wind velocity of  ~14.5 km s^-1 (as deduced from low rotational CO transitions, see Paper II), the ratio between the peculiar velocity of the star, v _⋆  and the wind velocity, v_w, is more than 1 and instabilities in the shell of the bow shock are expected (Dgani et al. 1996). This is also clear from the hydrodynamical simulations presented by van Marle et al. (2011) for parameters typical of a red supergiant resembling Betelgeuse. The non-occurrence of large-scale instabilities can put constraints on the physical parameters determining the morphology of the bow shock (see Sect. 4.2).

3.2.3. Dust temperature in the arcs and linear bar

Both the multiple arcs and the linear bar are best visible in the PACS images. The surface brightness in the multiple-arc like structure and linear bar has been measured for the different PACS filters using aperture photometry (see Table 1). Uncertainties are  ~15%, mainly dominated by the absolute flux uncertainty. The PACS and AKARI fluxes are consistent with each other for the arc-like region, with the exception of the AKARI 160 μm point which is a factor  ~3 higher. However, the flux measurements in the bar-like structure are quite different (see Fig. 12). We note that both the arcs and the bar are extremely faint in the AKARI LW images (see Appendix A), and contamination cannot be excluded. These data points will therefore be neglected.

Neither [O i] line emission at 63 μm nor [C ii] at 157 μm is detected with Herschel in the bow shock region or linear bar (see Paper II for more details), which is the reason we assume that dust emission is the main contributor to the Herschel/PACS images. Using a modified blackbody with a spectral index β varying between one (typical for, e.g., layered amorphous silicate grains or amorphous carbon grains, Knapp et al. 1993; Dupac et al. 2003) and two (typical of, e.g., crystalline silicate grains or graphitic grains, Tielens & Allamandola 1987; Mennella et al. 1998; Dupac et al. 2003), we derived the (mean) dust temperature from the Herschel data using a Monte-Carlo uncertainty estimation with normally distributed noise. For the arcs, a (mean) dust temperature of 87 ± 7 K (β = 1.04 ± 0.12) is derived and for the bar 64 ± 2 K (β = 1.0 ± 0.05) (see Fig. 12)². For the arcs, 91% of the β-values have values  < 1.1, and the maximum value for β is 1.88 with only 0.4% having values in the interval between 1.75 and 1.90. For the bar, the maximum β-value is 1.02. Low values for the spectral index have already been found in circumstellar environments (e.g., Knapp et al. 1993). Our result of a low spectral index (around 1) for quite high dust temperatures around 65 K (bar) and 85 K (arcs) is in good agreement with the inverse temperature dependence of the dust spectral index as derived by Dupac et al. (2003). According to Dupac et al. (2003), a low spectral index value might be explained by (1) the occurrence of large grains; (2) the fact that warm regions could harbor aggregates of silicates, porous graphite, or amorphous carbon, having a spectral index around 1; or (3) an intrinsic dependence of the spectral index on the temperature due to quantum processes, where in the case of high temperatures processes such as thermally-activated relaxation processes and temperature-dependent absorption could dominate. The chemical composition of the arcs might contain layered amorphous silicate grains (cf. Verhoelst et al. 2006) or larger grains (as recently detected around the AGB star W Hya, Norris et al. 2012) explaining the low β-index. In the case where the linear bar has an interstellar origin (see Sect. 5), amorphous silicates – an important component of the interstellar dust budget – might explain the low β-index as well. However, very small grains whose nature is thought to be carbonaceous (Desert et al. 1986) might also be the origin of the low β-index. Observations in the near-IR are prerequisites to constraining the detailed chemical content in the arcs and in the bar.

A modified blackbody with a temperature around 63 K can explain the WISE 22 μm data of the linear bar, but not the WISE 12 μm data. Admittedly, the flux in the WISE 12 μm band is quite uncertain and should be interpreted as an upper limit due to the large contribution of the central target (see Fig. 5). It might be that the Herschel and AKARI images show the emission by larger dust grains compared to the WISE data, which can be fitted with a modified blackbody of  ~145 K. Smaller grains indeed attain higher temperatures, hence emit at shorter wavelengths, than large grains, since the former absorb more efficiently per unit mass.

The dust temperature map (for β = 1) and hydrogen column density³ in the arcs and linear bar as derived from the Herschel PACS data are shown in Fig. 13. A lower limit for the flux values of 9, 8, and 4 MJy/sr was used for the 70, 100, and 160 μm data, respectively. Mainly in the eastern region, the flux values in the 160 μm image are too low to derive the dust temperature and hydrogen column density properly. More physical information on the eastern region can be obtained from either the slope of a linear fit to the flux values at 70, 100, and 160 μm or from the ratio of the flux at 70 μm to the flux at 100 μm (see bottom panels in Fig. 13).The arcs and bar have comparable dust colour temperatures. The dust temperature is higher both in the multiple arcs and in the linear bar in the direction of space motion of Betelgeuse. A gradient in temperature is visible, with the temperature decreasing (and column density increasing) at larger distances from the central target.

Considering energy conservation, one can show that the energy flux of the gas entering the shock front per unit area and per unit time is around five orders of magnitude lower than the stellar radiation energy flux at the distance of the arcs and bar. As a result, the additional heating of the dust grain from the hot post-shock gas is negligible, and the temperature of the dust grain is dictated by radiative equilibrium conditions when only taking the stellar radiation into account (see also Sect. 4.2.1). Assuming a smooth continuous outflow, the temperature of an iron or amorphous silicate grain only heated by the stellar radiation is in the order of ranges from about 45 K to 60 K at a distance of 280 − 530′′ away from the central star (see Appendix D). The correspondence between the mean dust temperatures derived from the data and these theoretical calculations is satisfactory, taken the assumptions of the theoretical model into account. Specifically, one can imagine that the dust temperatures in a clumpy envelope might indeed be somewhat higher than the values shown in Fig. D.1, since stellar photons can probably reach the outer envelope regions more easily. In addition, smaller grains will attain a higher temperature since they absorb more efficiently per unit mass.

3.2.4. Mass of the arcs and linear bar

For an opacity, κ_λ, of 185 cm²/g at 70 μm, 90 cm²/g at 100 μm, and of 35 cm²/g at 160 μm (Mennella et al. 1998) and assuming a dust-to-gas ratio of 0.002 (Verhoelst et al. 2006), we deduce a total gas and dust mass of  ~2.4 × 10^-3 M_⊙ for the arcs and  ~2.1 × 10^-3 M_⊙ for the linear bar. The derived (gas+dust) mass is very sensitive to the assumed dust temperature and dust-to-gas ratio. Lowering the dust temperature in both the arcs and bar to 30 K would yield a mass of 0.07 M_⊙ and 0.029 M_⊙, respectively. Assuming that the linear bar has an ISM origin (see Sect. 5) with a dust-to-gas ratio around 0.01, the mass in the bar would be around 0.01 M_⊙.

Assuming a smooth outflow, a first theoretical estimate for the mass in the bow shock shell can be derived from the equation (Mohamed et al. 2012) ${\mathit{M}}_{\mathrm{shell}}\mathrm{=}\frac{\mathit{Ṁ}\hspace{0.17em}{\mathit{R}}_{\mathrm{SO}}}{{\mathit{v}}_{\mathrm{w}}}\mathit{,}$(1)with R_SO the stand-off distance of the bow shock and v_w the wind velocity. This expression assumes the mass is distributed over 4π steradian in a spherical shell and only takes the contribution from the stellar wind into account. Thus for θ ≲ 90°, the expected shell mass would be half of M_shell. A bow shock shell mass of 2.4 × 10^-3 M_⊙ for the arcs would imply a mass-loss rate of 2 × 10^-7 M_⊙/yr (for a wind velocity of 14.5 km s^-1). However, CO and H i observations tracing the inner envelope regions yield a mass-loss rate that is one order of magnitude higher,  ~2 × 10^-6 M_⊙/yr (e.g., Paper II; Rodgers & Glassgold 1991; Le Bertre et al. 2012). This hints that the mean mass-loss rate might have varied substantially, as also suggested in Sect. 3.1.

The hydrodynamical models presented in Sect. 4.2.1 yield a mass in the arcs of  ~0.1 M_⊙ (including 90% CSM and 10% ISM material). According to Mohamed et al. (2012), the observed low mass in the arcs might indicate that the bow shock created by the RSG wind has a very young age, of the order of 20 000 yrs, and may not yet have reached a steady state. But, the existence of (non-homogeneous) clumps in the inner wind travelling all the way to the bow shock region might also explain the low mass in the arcs as derived from the Herschel images.

4. Discussion: origin of multiple arcs

In this section, we focus on the origin of the multiple arcs seen in the bow shock around Betelgeuse. Section 4.1 gives details on some observational properties. In Sect. 4.2, the observational data are compared with hydrodynamical simulations. In Sect. 4.3, we discuss the origin of the multiple arcs based on the observational and theoretical constraints.

4.1. Constraints from the observations

From the observations described in previous sections, two concerns might be important for constraining the origin of the multiple arcs seen in the bow shock region of Betelgeuse. (1) The different arcs are reasonably well fitted with concentric ellipses (see Fig. B.1 in the Appendix B). (2) The inner envelope shows clear evidence of a pronounced asymmetric clumpy structure. This leads to the hypothesis that the origin of the different arcs is a projection effect on the plane of the sky of the contact discontinuity created from the collision between different inhomogeneous mass-loss events and the ISM. As shown in Fig. 6 of Cox et al. (2012), different inclination angles lead to projected outlines of the contact discontinuity, where concentric ellipses with a smaller inclination angle yields a smaller projected size of the bow shock outline. This would mean that the arcs are shaped by an interaction of the inhomogeneous stellar mass loss at past and present epochs.

Another explanation for the different arcs (and also linear bar) might be inferred from the recent results of van Marle et al. (2011), who shows that small grains follow the movement of the gas mass elements, while larger CSM grains tend to penetrate deeper into the ISM (see also Sect. 4.2). This analysis takes only the inertia of the grains into account. However, since the grains are charged, they will gyrate around the magnetic field lines. Assuming a gas temperature in the shock of 10 000 K and a shock velocity equal to the space velocity, the typical Larmor radius for a grain with a size of 5 nm is  ~1 × 10¹⁴ cm, while larger grains with a size of 1 μm have a Larmor radius in the order of 4 × 10¹⁸ cm. With a typical width of  ~20′′ for the arcs, grains with a size  ≲0.1 μm are position-coupled to the hot gas. One then could postulate that arcs further away from the central target could contain larger dust grains (≳0.1 μm), while the inner arcs contain the smallest grains coupled to the gas. In the case of randomized mass-loss variations (as suggested by the data tracing the inner envelope), one could consider the possibility that each mass-loss event has a favourable grain-size spectrum (i.e., the grain size distribution function reaches a maximum for some specific grain sizes). An analogous situation is obtained by Fleischer (1994), who modelled the dust formation in the case of a dynamical atmosphere and showed that far out in the wind, where the grain growth is definitely stopped, the grain-size spectrum shows distinct grain-size peaks. A larger grain size implies a higher drift velocity and a larger Larmor radius, which eventually might lead to a bow shock arc further away from the central target. However, it remains remarkable that only Betelgeuse in the whole MESS sample shows the appearance of multiple arcs in the bow shock region, while the same argument as given above might hold for other targets in the sample, especially the Mira-type variables with high-amplitude variations being favourable for a grain-size spectrum with strong distinct peaks. In addition, the dust temperature maps (Fig. 13) show no significant temperature difference between the different arcs, while it is well known that smaller dust grains are generally hotter than larger grains since the former absorb more efficiently per unit mass.

The observed fragmentation of the bow shock might also be understood in terms of the effects predicted by Dgani & Soker (1998). They postulate that RT instabilities might fragment a bow shock in the direction of motion, enabling the ISM to flow into the inner parts of the envelope and creating a bow shock well inside the almost spherical, but very filamentary, haloes. For targets close to the Galactic plane (as is the case of Betelgeuse, with b =  −9°), the ISM magnetic field can make RT instabilties very efficient. This might explain the recent detection by Le Bertre et al. (2012) of a detached H i shell with a radius of 2′. Magnetic fields would suppress certain modes but accentuate others, changing the appearance of the envelope, and potentially leading to “RT rolls” (or stripes). An example resembling the image of Betelgeuse in shown in Fig. 3a in Dgani & Soker (1998). We should remark, however, that the calculations as performed by Dgani & Soker (1998) were for the case of a planetary nebula, which have much higher stellar wind velocities than red supergiants or AGB stars. In the sample of 78 AGB stars and supergiants analysed by Cox et al. (2012), Betelgeuse is the target closest to the Galactic plane (z =  −5 pc), and it is the only target showing this multiple arc-like structure. Four other targets in the MESS sample are at similar distances from the Galactic plane (z < 20 pc), one being a supergiant (μ Cep) which resembles Betelgeuse in different aspects, including a similar dust mass-loss rate to Betelgeuse (Josselin & Plez 2007). The main difference for parameters influencing the bow shock morphology is the wind velocity, which is a factor 2.4 higher in μ Cep (De Beck et al. 2010). While for Betelgeuse v _⋆ /v_w > 1, this ratio is  <1 for μ Cep in which case it has been predicted by Dgani et al. (1996) that the bow shock is stabler. The Herschel images presented by Cox et al. (2012) show, that (possible) instabilities in the bow shock of μ Cep have comparable sizes to Betelgeuse, while different arcs are indeed not seen.

If present, the strength of the ISM magnetic field can be predicted from the angular separation between the arcs. An angular separation of 30′′ at a distance of 300′′ yields an Alfvén speed in the pre-shocked ISM of  ~4 km s^-1 (Eq. (4) in Dgani 1998). For an ISM density of 1.9 cm^-3 (Cox et al. 2012) and assuming that mainly H i is present in the bow shock, this yields a magnetic field of 3 μG. The average strength of the total magnetic field in the Milky Way is about 6 μG near the Sun (Biermann 2001) and increases to 20–40 μG in the Galactic centre region. Dense clouds of cold molecular gas can host fields of up to several mG strength (Heiles & Crutcher 2005).

4.2. Constraints from hydrodynamical modelling

For insight into the origin of the complex CSM-ISM interface and the non-occurrence of large-scale instabilities in the arcs, we have simulated the hydrodynamical evolution of the envelope of Betelgeuse. The simulations also show if dust grains are good tracers for the gas-driven dynamics in the bow shock. The hydrodynamical computations were done using MPI-AMRVAC code (Keppens et al. 2012; van Marle et al. 2011). More details on the numerical method and initial conditions are given in Appendix E.

The simulation was run for five different situations, using variations of a basic model to explore several aspects of the bow shock morphology and its dependence on the properties of the wind and the ambient medium. The physical parameters for the basic model are given in Table 2. The gas density in the ambient medium is set at 3 × 10^-24 g cm^-3 (or n_H ≈ 2 cm^-3) and a temperature of 100 K (reflecting the temperature of a cold neutral medium). The five different simulations have the following specifications:

• The first simulation (A) uses the basic wind parameters that arespecified in Table 2 and the ambient mediumspecified above. We assume that all dust grains have the same size,5 nm, and density, 3.3 g cm^-3.

• The second simulation (B) is the same as simulation A, except for the size of the dust grains, which is increased to 100 nm to demonstrate the behaviour of large dust grains.

• The third simulation (C) has a warm ambient medium of 8000 K, reflecting the temperature of a warm neutral or partially ionized medium. The ambient medium is assumed to be heated by an outside radiation source. Therefore, whereas the ambient medium is kept at a minimum temperature of 8000 K, the stellar wind, which is protected from UV-radiation by the ionized gas of the bow shock, has a minimum temperature of 100 K.

• The fourth simulation (D) has the same ambient medium temperature as simulation C, but with lower ambient medium density (5 × 10^-25 g cm^-3 ) and higher space velocity (v_⋆ = 72.5 km s^-1). This simulation was added to compare our results with the recent hydrodynamical models of Mohamed et al. (2012).

• For the fifth simulation (E), as for the first and second, we set the minimum temperature to 100 K throughout the entire domain. However, instead of a smooth wind, we assume a periodic variation in the mass loss rate. For 1000  years the mass-loss rate is assumed to be high (Ṁ = 2.5 × 10^-5 M_⊙/yr), then, for the next 9000 years, it is low (Ṁ = 2.5 × 10^-8 M_⊙/yr). This “picket-fence” type variation is then repeated periodically. The two mass-loss rates are related as Ṁ_high    =    1000   Ṁ_low. The total amount of mass lost over the 10 000 year period is the same as for the first two simulations. The period of 10 000  years reflects the typical timescale for AGB thermal pulses.

The effects of changing one or more input parameters in the simulations are discussed in Sects. 4.2.1–4.2.3. A link to the online movies is provided in the online Appendix F. A reflection on the origin of the multiple arcs based on the outcome of the simulations is given in Sect. 4.3.

4.2.1. Cold ambient medium

The result for simulation A with an ambient medium temperature of 100 K is shown in Fig. 14. After a simulation time of 300 000 yrs, the place of the bow shock interaction has stabilized: the termination shock⁴ occurs at 0.25 pc (or 262′′) and the bowshock⁵ at 0.34 pc (or  ~ 356′′), in good agreement with the observations (see Sect. 3.2, neglecting the inclination angle). The turbulent astropause/astrosheath⁶ thus has a width of  ~90′′. Although both the wind termination shock and the forward shock are smooth, the contact discontinuity, which separates the shocked wind from the shocked ambient medium, shows instabilities. These do not penetrate the shocks on either side, but cause deviations in the shape of the forward shock.

Figure 14 shows three types of instabilities in the gas. (1) A thin layer of shocked wind material, right on the contact discontinuity, has a higher density and lower temperature than the surrounding gas. This is the result of the radiative cooling instability, which causes dense gas to cool faster. As a result, dense gas loses thermal pressure and is compressed, which increases the density and therefore the cooling rate. (2) The high density gas of the shocked wind extends “fingers” into the shocked ambient medium, which form in the forward shock and increase in size as they travel downstream. These are Rayleigh-Taylor instabilities, resulting from low density, but high-pressure, material (the shocked ambient medium) pushing back denser shocked wind. (3) Due to the shear-force between the shocked wind and the shocked ambient medium, the Rayleigh-Taylor instabilities deform and show a circular motion. These are Kelvin-Helmholtz instabilities and are primarily visible downstream from the star. All three features were also found by van Marle et al. (2011).

The temperatures in the shocked gas are lower than the kinetic energy of the free-streaming wind and the ambient medium. For a fully adiabatic shock, the temperatures at the front of the bow shock would be about 27 300 K for the shocked interstellar medium and around 7200 K for the shocked wind. Instead, the maximum temperature for the shocked interstellar medium is 14 000 K and the shocked wind does not reach more than 4000 K. van Marle et al. (2011) find similar temperatures. This is mainly due to the effective radiative cooling. Also, as the thermalized gas moves downstream, it decompresses, trading temperature for an increasing velocity, which causes the tail of the shock to be cooler than the front. This can be seen by comparing the absolute velocity and temperature of the gas (right panel of Fig. 14).

Comparing the gas and dust density (left panel of Fig. 14) shows that the dust particles penetrate the shocked wind layer, but are brought to a stop at the contact discontinuity, which they tend to follow. In that sense, small dust grains are good tracers of the gas-driven dynamics at the contact discontinuity. This confirms the earlier results found by van Marle et al. (2011). As a result, the instabilities at the contact discontinuity should also be visible in the infrared (assuming that the observations have sufficient resolution, which is the case for the Herschel instruments). Since the dust is concentrated at the contact discontinuity, other features that are typical of the shocks would most likely not show up in infrared images. For example, in the region behind the star, the wind termination shock curves back toward the axis of motion as its thermal pressure balances against the ram pressure of the wind. This effect is completely invisible in the dust. The behaviour of the circumstellar environment over the full 200 000 year period, starting from the moment when radiative cooling is introduced, is presented online in animated form (Appendix F). This movie shows that the number of instabilities increases with the introduction of radiative cooling. These instabilities form in the front of the bow shock and then travel downstream.

The total mass of the gas trapped inside the bow shock region, found by integrating the density between the termination shock and bow shock for 0 ≤ θ ≤ 90°, is approximately 0.1   M_⊙ of which 90% is shocked wind material. The shell mass is considerably higher than the value obtained from the Herschel data (see Sect. 3.2). As already suggested, a possible explanation may be that the envelope around Betelgeuse is very clumpy or that Betelgeuse is a very young RSG star so that the bow shock is not yet fully formed (Mohamed et al. 2012). This could also explain the smoothness of the bow shock, since instabilities need time to form. Most of the mass is situated at the sides of the star, where the volume of the shock region is much larger; e.g., when integrated only over 0 ≤ θ ≤ 45° the total mass becomes 1.8 × 10^-2 M_⊙, of which about one third is shocked ISM.

Assuming radiative-equilibrium balance, the heating of the dust grain due to energy transfer from the hot post-shock gas can be estimated. An approximation is given by Eq. (5.40) from Tielens (2005), which assumes the absorption efficiency of the dust as a function of wavelength to be a simple power law: Q(λ) = Q₀(λ/λ₀)^β, with Q₀ = 1, λ₀ = 2π   a and β = 1. As a result, the dust temperature only depends on the mean intensity of the radiation field J and the size of the dust grains a, ${\mathit{T}}_{\mathrm{d}} \mathrm{\propto } (\frac{\mathit{J}}{\mathit{a}}{)}^{\mathrm{0.2}}\mathrm{·}$(2)The total energy flux entering the shock front is $\begin{array}{ccc}{\mathit{F}}_{\mathrm{total}}& \mathrm{=}& {\mathit{F}}_{\mathrm{wind}}\mathrm{+}{\mathit{F}}_{\mathrm{ISM}}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathit{Ṁ}\frac{{\mathit{v}}_{\mathrm{w}}}{\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\rho }}_{\mathrm{ISM}}\hspace{0.17em}{\mathit{v}}_{\mathit{\star }}^{\mathrm{3}}\\ & \mathrm{\approx }& \mathrm{5.5}\mathrm{\times }{\mathrm{10}}^{-5}\hspace{0.17em}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}\mathit{.}\end{array}$(3)Assuming that only 20% of the energy flux is converted to radiation, the additional heating from the dust grain due to the hot post-shock gas is only  ~6 K, while the heating due to the stellar radiation is one order of magnitude higher (see Fig. D.1 in the Appendix D).

The effect of grain size on the distribution of dust in the bow shock region is demonstrated in Fig. 15. This shows the result for simulation B for a larger grain with size 100 nm. Instead of following the contact discontinuity, the larger dust grains are concentrated at the forward shock and, in the tail of the bow shock, even penetrates the undisturbed ambient medium. This behaviour, which was also shown in van Marle et al. (2011), is the result of the larger momentum of large grains, compared to the drag force. The large grains are less tightly coupled to the gas and can therefore penetrate deeper into the bow shock region and the ambient medium beyond. However, since larger grains are less numerous than smaller ones, the infrared morphology of the bow shock will mainly be determined by the smaller dust species.

4.2.2. Warm ambient medium

When the ambient medium is kept at a minimum temperature of 8000 K (simulation C), the results look different, as shown in Fig. 16. The Rayleigh-Taylor instabilities are much smaller. This is the result of radiative cooling. The instabilities, which are made up of stellar wind material, are allowed to cool to 100 K, whereas the surrounding shocked ambient medium has a minimum temperature of 8 000 K. Therefore, the gas in the Rayleigh-Taylor fingers has a lower thermal pressure than the surrounding gas and is compressed, inhibiting their growth. This is different from the results shown in Cox et al. (2012), which also showed large instabilities for the warm ambient medium. However, these may have been partially due to numerical effects of modelling a spherical wind expansion on a cylindrical grid, which lead to a large instability close to the symmetry axis.

The layer of shocked interstellar matter extends further from the axis of motion and has an almost triangular shape. The combination of the smaller instabilities and the larger distance between contact discontinuity and forward shock also reduces the effect of the instabilities on the shape of the shock, which is completely unaffected.

Because of the high thermal pressure of the ambient medium, the wind region is more constrained; in fact, the free-streaming wind has not yet reached the lower boundary, but is still constrained by the shocked wind. This was also seen in the results of Cox et al. (2012). As in the first simulation, the small dust grains follow the contact discontinuity.

The result of simulation D with the same high ambient medium temperature, but lower ambient medium density and higher space velocity, is shown in Fig. 17. The high stellar velocity creates a stronger shear force between the shocked wind and the shocked ambient medium. As a result, the Kelvin-Helmholtz instabilities are much greater, leading to a large-scale instability that changes the shape of the entire bow shock. This is contrary to the results of Mohamed et al. (2012), who obtained a smooth bow shock for similar input parameters. This difference may be the result of a different numerical treatment (grid versus SPH), but may also result from the total timescale of the simulation, because Kelvin-Helmholtz instabilities need time to grow (see accompanying movie in Appendix F). Owing to the larger collision speed at the bow shock, the shock temperatures are much higher (~10⁵ K) than for the other simulations. This largely negates the effect of the high ambient medium temperature, since its thermal energy is much lower than the kinetic energy. While the general structure of the bow shock is unstable, the local dust shell directly in front of the star is quite smooth. Also, the instability of the bow shock leads to a structure that curves in the opposite direction from the wind termination shock. Similar features have been found by Wareing et al. (2007a) and Cox et al. (2012).

4.2.3. Clumpy mass-loss episodes

The hypothesis that the different arcs might be the result of the collision of different clumpy mass-loss episodes having filling factors  < 1 with the surrounding ISM (see Sect. 4.1) can be checked using the hydrodynamical simulations. In the first instance, this can be done in a simplified 2D scheme, where the collision between the ISM and a CSE with a variable mass-loss rate has been modelled. A full 3D model is beyond the scope of this paper. As outlined above, we simulated this effect by varying the mass-loss rate by a factor 1000 over a period of 10 000 yrs. Since the change in dust density at the edge of a clump is very steep, the mass-loss rate variations were implemented in a “picket-fence” way.

Figure 18 shows the result of simulation E with variable mass loss. The mass leaves the star concentrated in shells that are separated by voids with a density that is three orders of magnitude lower. As they travel outward, the shells expand owing to their own internal pressure. In the wake of the star the shells eventually merge. To the front of the star this does not occur because the shells hit the wind termination shock before they can expand far enough.

The variation in mass-loss rate (and therefore in ram pressure) causes a widespread disturbance in the morphology of the bow shock region. The wind termination shock, where the isotropic thermal pressure of the shocked gas is equal to the ram pressure of the wind, shows a wave-like pattern in the downwind region as it conforms to the variable ram pressure of the wind. At the front of the bow shock the variation in the free-streaming wind increases the instability at the contact discontinuity. The resulting Rayleigh-Taylor and Kelvin-Helmholtz instabilities are more numerous than in the case of a smooth wind (Fig. 14) and show more small-scale structure.

Finally, the impact of subsequent shell collisions means that multiple arc-like structures appear between the termination shock and the contact discontinuity (for a detailed view see Fig. 19). This may be an indication of the origin of the multiple arcs observed in the bow shock of Betelgeuse.

The animation of the simulation (Appendix F) shows the influence of the variable mass-loss rate, in terms of both the formation and growth of local features, as well as of the large-scale morphology of the bow shock. Despite the variation in ram pressure, the location of the bow shock mostly remains unchanged. This is the result of the period of the variation relative to the dynamical timescale of the bow shock region. The bow shock location can only change if the dynamical timescale is significantly shorter than the period of the variation. This is not the case in this simulation. The bow shock region has a cross-section of the bow shock of approximately 0.1 pc at the front and up to a parsec to the sides of the star. With a sound speed of about 10 km s^-1 (for a temperature of 10 000 K), the timescale varies between 10 000 and 100 000 years. Since the period of the variations is only about 10 000 years, the bow shock does not have enough time to adjust to the variations.

4.3. Discussion

A. Multiple arcs:

the hydrodynamical simulations show some good agreement between the observations and the theoretical predictions, which might point us toward the origin of the arcs; e.g., the observed ratio between the projected – minimal – distance between the star and the bow shock outline in the direction of relative motion and the distance perpendicular to that, i.e. at θ = 90°, is  ~0.66 for the different arcs. The same ratio is obtained in simulations A, C, and E, while the ratio is significantly lower for simulations B and D.

The hydrodynamical simulations show that the formation of dust clumps in the inner wind regions, simulated in all its generality by implying spherically symmetric mass-loss variations, can have an impact on the observed morphology of the bow shock region, and potentially might lead to the formation of multiple arcs in the bow shock. Dust clumps with higher grain size have a larger drift velocity, large Larmor radius, and longer stopping length when passing through the bow shock, and eventually might create a bow shock arc further away from the central target. However, as said above, we would then expect that other targets in the MESS sample would show multiple arcs in the bow shock region as well.

An effect not included in the MPI-AMRVAC-code is photo-ionization. In the case of the CO waves (or ripples) detected on the surface of the Orion molecular cloud, Berné et al. (2010) argues that ultraviolet radiation has created a (small) insulating photo-ablative layer, allowing the development of a KH instability with a wavelength at least one order of magnitude longer than the insulating layer. A speculative idea is that the multiple arcs found around Betelgeuse are undulations created by photo-chemical effects stimulating the growth of some instability with wavelength around the observed width of the arcs, which are then traced by the smaller dust grains that tend to follow the gas instabilities. However, also for this effect, it remains surprising that Betelgeuse is the only target in the MESS sample showing multiple arcs.

In our opinion, the most probably reason for the origin of the multiple arcs is the combined effect of (1) a turbulent non-homogeneous mass loss, which might enable the ISM to flow into the envelope and create multiple bow shocks, together with (2) the effect of the interstellar magnetic field, which can make some instabilities very efficient while suppressing other modes (see Sect. 4.1). The angular separation between the arcs points toward a magnetic field of  ~3 μG.

B. Suppression of large-scale instabilities:

remarkably, the size of (local) density clumps tracing potential instabilities inside the arcs as seen in the Herschel images (with an upper limit to the maximum length of  ~30′′, see Sect. 3.2) is often smaller than predicted by the hydrodynamical simulations. Most of the instabilities seen in Figs. 14 and 19 are Rayleigh-Taylor (RT) instabilities, which are slightly deformed by the action of Kelvin-Helmholtz (KH) instabilities. Thin-shell instabilities such as the NTSI (Vishniac 1994) or TAI (Dgani et al. 1996) do not occur, since the shell between both shocks is too thick: with a typical thickness of the shell in the order of 0.1 pc, the crossing time is in the order of 10⁴ yrs (for a velocity of  ~10 km s^-1) by which time an individual mass element is already displaced by a few times 10¹⁷ cm, preventing the growth of an instability. Thermal instabilities are also ruled out, since radiative shocks with a power-law cooling function (Λ ∝ T^γ) are stable against thermal instabilities if γ ≳ 1 (Gaetz et al. 1988).

The only example in our simulations showing only small-scale instabilities is the case of a warm ambient medium around Betelgeuse (with a temperature of  ~8000 K) for which a strong growth of the instabilities is prevented by compressing the Rayleigh-Taylor fingers. As shown in Fig. 16, instabilities in the direction of the space motion of the target then have typical lengths  ≲ 1 × 10¹⁷ cm (or 34′′).

We should also realize that the size and growth rate of the instabilities are directly related to the density of the gas (Chandrasekhar 1961). In our simulations, compression due to radiative cooling increases the density of the shocked wind. However, it is possible that we over-estimate the cooling. The energy loss due to radiative cooling is proportional to the ion density times the electron density. If the gas is only partially ionized (as in the 3000 − 8000 K regime), the cooling may be reduced, which in turn would reduce the compression and therefore the density contrast at the contact discontinuity, damping the instabilities.

Not considered in our simulations is the effect of a magnetic field. Depending on the direction of the magnetic field with respect to the shear flow, a magnetic field might (de)stabilize KH instabilities (Miura & Pritchett 1982; Keppens et al. 1999). Miura & Pritchett (1982) found that in a sheared magnetohydrodynamical flow in a compressible plasma, modes with kΔ < 2 are unstable, with Δ the scale length of the shear layer and k the wavenumber. Compressibility and a magnetic field component parallel to the flow (B₀ ∥ v₀) are found to be stabilizing effects. For the transverse case (B₀ ⊥ v₀), only the fast magnetosonic mode is destabilized, but if k·B₀ ≠ 0, the instability contains Alfvén-mode and slow-mode components as well. Keppens et al. (1999) studied the case of an initial magnetic field aligned with the shear flow. In the case the initial magnetic field is unidirectional everywhere (uniform case), the growth of KH instabilities is compressed, while a reversed field (i.e. a field that changes sign at the interphase) acts to destabilize the linear phase of the KH instability. As a result, the non-occurrence of large-scale instabilities might potentially be explained by the presence of a magnetic field. As suggested in Sect. 4.1, the detection of multiple arcs in the bow shock might also point towards the presence of an ISM magnetic field. With a strength of  ~3 μG, the magnetic pressure is  ~3.6  ×  10^-13 dyn cm^-2, being two orders of magnitude lower than the ram pressure at the contact discontinuity. Without dedicated simulations, it is not possible to predict how the presence of a weak ISM magnetic field might influence the growth rate of the instabilities.

One might also argue that the bow shock shell is still very young so that instabilities might not have had enough time to grow. Using the prescriptions of Brighenti & D’Ercole (1995) the typical growth time for RT instabilities to reach a size of 30′′ at a distance of  ~200 pc is only  ~ 7000 yrs, while the KH growth timescale is at least an order of magnitude smaller⁷. However, according to our simulations, the shell would not have had time yet to reach the observed stand-off distance. Changing the mass-loss rate by a factor of  ~10 only reduces the current stabilization time from  ~100 000 yrs to  ~30 000 yrs.

The clumpy structures found in the inner envelope region might also have an imprint on the appearance of the bow shock. With shock velocities below 50 km s^-1, the percentage of dust mass destroyed (i.e. returned to the gas phase) is negligible (Tielens 2005). But, one should also take into account that dust clumps tend to dissipate. One possible depletion mechanism is heat conduction, which is thought to be unimportant in the case of cool supergiant (and AGB) winds owing to the small difference in temperature between the clump and its environment. According to Hartquist et al. (1986), clouds are, however, ablated when they are in relative motion with respect to the surrounding medium. Following Smith et al. (1988), we estimate that all dust clumps with a size smaller than  ~10¹⁵ cm (or 0.3′′ at a distance of 200 pc) will be ablated by the time the clump reaches the bow shock region; larger dust clumps might still be visible in the bow shock region. The effect of clumpiness on the development of instabilities has not yet been studied. In the case where most of the dust grains formed in the inner wind had already been ablated by the time they arrived at the bow shock region, the clumps seen in the outer arcs might also reflect ISM inhomogeneities progressively engulfed by the expanding shell.

5. Discussion: Origin of the linear bar

It is striking that the arcs and the linear bar show the same far-infrared colour temperature and density (as deduced from the PACS data). This might suggest that they are illuminated by the same source, probably Betelgeuse itself, and belong to the same physical phenomenon (i.e. the bow shock). As already argued in previous section, larger grains are capable of leaving the shocked wind layer, and can penetrate the shocked or even unshocked ISM. The large grains almost completely ignore the morphology of the shocked gas shell and penetrate the unshocked ISM a considerable distance (~0.1–0.2 pc, van Marle et al. 2011). Since they are very low in number, they will not influence the gas morphology. However, even for very large grains, the shape of the bow shock still shows some curvature (see Fig. 4 in van Marle et al. 2011), which is not seen in the Herschel images. Recently, Mackey et al. (2012) have suggested that the linear bar is the relic of a blue supergiant (BSG) wind interaction with the ISM, where Betelgeuse has recently entered the red supergiant (RSG) phase. The BSG reverse shock collapses in on the RSG wind producing an inner shocked shell (the arcs in the Herschel image). The r^-2 decrease in dust emission makes the contact discontinuity look a little more like a bar (compared to the projected density in a 2D simulation) because the strongly curved parts at larger radii are fainter. If so, Mackey et al. (2012) predict that the bar should eventually begin to curve back into a bow shock shape, with a radius of curvature that can be larger than the current stand-off distance of the bar suggests. The large FOV Herschel images cannot give extra support to this proposed scenario since no curvature is seen.

We thus should consider the possibility that the linear bar is not related to Betelgeuse, but is an interstellar structure, which is by accident co-spaced with Betelgeuse. If so, the motion of the star and angular separation between the star and bar imply a collision between the two in about 17 500 yrs, while the arcs and linear bar would collide in about 5000 yr.

The linear bar might be the edge of an interstellar cloud. The results obtained by Dickey & Garwood (1989) prove that it is not unreasonable that Betelgeuse is near an interstellar cloud. Dickey & Garwood (1989) showed that a few hundred clouds are present within 1 kpc having a column density around 10¹⁹ cm^-2 (cf. Fig. 13). With an estimated mass for the linear bar around  ~0.01 M_⊙, the molecular cloud mass spectrum has a density of  ~10 pc^-3 (see Fig. 2 in Dickey & Garwood 1989). Using H i 21 cm observations, Matthews et al. (2011) has recently detected an extended circumstellar wake associated with the CSE around the AGB star X Her. Mapping shows that the star lies (in projection) near the periphery of a much larger H i cloud that also exhibits signatures of interaction with the ISM. While the radial velocity of the cloud near X Her overlaps that of the CSE, this is not the case for α Ori (see Fig. 7). However, the probability of a mere chance superposition in position, velocity, and direction of space motion between a cloud and target is extremely low. Betelgeuse might thus be close to an interstellar cloud, but have a different radial velocity.

The bar might also be a linear filament in the interstellar cirrus, as also seen in the Galactic centre, whose possible origin is linked to the Galactic magnetic field (Yusef-Zadeh et al. 1984; Yusef-Zadeh & Morris 1987). Linear filaments might be caused by particle acceleration by shocks related to preexisting filamentary magnetic structure in the ISM (Heyvaerts et al. 1988) or by particle energization at reconnection sites at the ionized surfaces of molecular clouds and subsequent particle loading of ISM magnetic field lines (Serabyn & Morris 1994). Filamentation can also be the result of synchrotron thermal instabilites (Bodo et al. 1990): both the thermal condensation mode and slow magnetosonic wave can create structures that are aligned with the background magnetic field. Rosner & Bodo (1996) show that a particle acceleration process, akin to the acceleration of the anomalous cosmic-ray component associated with the solar wind termination shock, assisted by radiative instabilities driven by synchrotron emission can provide an explanation for the filamentation process. Some fraction of accelerated particles will be line-tied to the external (ISM) magnetic field. The local field magnetic field is parallel to the Galactic plane and follows the local spiral arms (Andreasyan & Makarov 1989), having a pitch angle of about  − 8° (Han 2008). Optical polarization measurements in the direction of Betelgeuse on a scale of several degrees show that the position angle of the local ISM galactic field (Bingham & Shakeshaft 1967; Heiles 2000; Heiles & Crutcher 2005) is consistent with the position of the linear bar. However, typical temperatures of interstellar cirrus clouds are around 15 K (e.g. Miville-Deschênes et al. 2010), much lower than the temperatures we deduced from the Herschel maps. Only when the filamentary cirrus is by chance close to Betelgeuse can the radiation of the target illuminate and heat the filamentary structure.

6. Conclusions

This paper presented new Herschel PACS and SPIRE images of the red supergiant Betelgeuse. The Herschel data show that the collision between the wind of the red supergiant and the interstellar medium has created a very complex structure. Multiple arcs at a distance of 6–7′ are detected with typical dust temperatures in the range of 40–140 K. At a distance of  ~9′ from the central target, an intriguing linear bar is detected, with the same colour temperature as the arcs. Betelgeuse is the only target in the MESS sample showing these two phenomena.

The inner envelope structure is very inhomogeneous and testifies to a turbulent mass-loss history. The Herschel images show the first evidence of a high degree of non-homogeneous distribution for the material beyond 15′′, which even persists until the material collides with the ISM. This inhomogeneity is probably linked to the giant convection cells recently detected in the outer atmosphere of Betelgeuse, which might enable localized dust creation and ejection.

The Herschel data have been complemented with data from the UV (from GALEX), near-IR (from WISE), and radio (from the GALFA-HI survey) to disentangle the different components that might result in this complex bow-shock structure and in the linear bar. The average temperature in the arcs is  ~85 K and in the bar is  ~63 K. A low dust emissivity index, β, around 1 is deduced from fitting the spectral energy distribution. This might point toward (1) the occurrence of large grains; (2) a grain composition resembling layered amorphous silicates grains, porous graphite, or carbonaceous grains; or (3) high temperature processes influencing the quantum effects. The dust temperature in the arcs and bar is regulated by the stellar radiation field of Betelgeuse, while heating due to the hot post-shock gas is negligible.

The Herschel PACS images of the inner envelope and the H i results by Le Bertre et al. (2012) suggest that the mean dust and gas density might have changed considerably  ~32 000 yrs ago (radius of  ~2′). The idea of a strong variation in the mean density is also supported by the low bow shock shell mass deduced from the Herschel data. This change in mean density might be from a change in mean mass-loss rate or the formation of an inner bow shock well inside a filamentary halo, in which case the ISM is able to flow into the inner parts of the envelope creating a bow shock.

Detailed hydrodynamical simulations were performed to understand the morphology of the bow shock around Betelgeuse. The stand-off distance and the ratio between the project distance between the star and the bow shock outline in the direction of relative motion and the distance perpendicular to that are predicted well by the hydrodynamical simulations, except in the case of large dust grains (around 100 nm; simulation B) or when a low ambient medium density and high space velocity are considered (simulation D). It is clear that small dust grains are good tracers for the gas-driven dynamics at the contact discontinuity. Most of the mass in the bow shock is situated to the sides of the star, and only  ~20% is situated in a cone with angular extent  <45°. The dust temperature in the bow shock is mainly determined by heating by the stellar radiation, while energy transfer from the hot post-shock gas can be neglected. Increasing the temperature in either the ambient interstellar medium or in the shocked wind prevents the growth of strong instabilities. Allowing for strong variations in the mass-loss rate in the inner envelope has an impact on the appearance of the bow shock, and eventually multiple arc-like structures might appear between the termination shock and the contact discontinuity.

Based on the observations and on hydrodynamical simulations, different hypotheses are formulated to explain the origin of the multiple arcs and to understand why no large-scale instabilities are seen in the bow shock region. In our opinion, the two main ingredients to explain both features are (1) a clumpy mass-loss process and (2) the influence of the Galactic magnetic field. The occurrence of clumps might also explain the low total (gas and dust) mass in the arcs as deduced from the data (~2.5 × 10^-3 M_⊙). The Galactic magnetic field might serve as an extra trigger to fragment the bow shock region. The angular separation between the arcs seen in the Herschel data is compatible with a magnetic field of  ~3 μG, in agreement with the average strength of the magnetic field in the Milky Way.

The linear bar might be the edge of an interstellar cloud illuminated by Betelgeuse or a linear filament whose a possible origin is linked to the Galactic magnetic field. Since no curvature is present in the bar, we believe that the bar is not directly linked to a previous blue supergiant wind, as recently discussed by Mackey et al. (2012).

Online material

Appendix A: Previous images of the bow shock around Betelgeuse

For completeness, we show in this section the bow shock around Betelgeuse as observed with IRAS by Noriega-Crespo et al. (1997) (Fig. A.1) and AKARI by Ueta et al. (2008) (Fig. A.2).

Appendix B: Additional figures

Figure B.1 PACS 70 μm image on which concentric ellipses were fitted to the different arcs.

Appendix C: Herschel PACS image of R Leo

The Herschel PACS 70 μm image of the AGB star R Leo are compared to a hydrodynamical simulation in Fig. C.1. Rayleigh-Taylor instabilities, slightly deformed owing to the action of Kelvin-Helmholtz instabilities, are visible both in the data and in the simulations.

Appendix D: Dust grain temperature in the envelope around Betelgeuse

Under the assumption of a smooth continuous outflow and using the parameters as specified in Table 2, the dust grain temperatures for Fe and amorphous silicates in the envelope around Betelgeuse have been calculated using the mcmax code (Min et al. 2009). Only the stellar radiation is taken into account as the energy source to compute the dust temperatures assuming radiative equilibrium. The particle shape is represented by CDE (continuous distribution of ellipsoids Bohren & Huffman 1983) for particles of which the equivalent spherical grain would have a grain size of 0.1 μm. The composition of the amorphous silicates is taken from de Vries et al. (2010), who derived that the infrared spectrum of the AGB star Mira is best fitted with 65.7% MgFeSiO₄ and 34.3% Mg₂SiO₄, implying an average composition of Mg_1.36Fe_0.64SiO₄, having a β-index of  ~1.8. The resulting dust temperatures are shown in Fig. D.1. The temperature for the amorphous silicate grains is lower than for the Fe-grains, owing to the higher opacity of Fe at the near-IR wavelengths, where the star emits most of its photons.

Appendix E: Numerical simulations of the bow shock

E.1. Numerical method

For the simulations the MPI-AMRVAC hydrodynamics code (Keppens et al. 2012) was used. For a comparison to similar work, we refer to van Marle et al. (2011). The MPI-AMRVAC code solves the conservation equations for mass, momentum, and energy on an adaptive mesh grid. It includes optically thin radiative cooling, using the exact integration algorithm (Townsend 2009), and a cooling curve for gas at solar metallicity (Wang Ye, priv. comm.). For these simulations, a two-dimensional approach was taken, assuming an inclination angle of 0°. As deduced by Cox et al. (2012), the inclination angle of the space motion with respect to the plane of the sky is small (~8°).

Dust was included in the simulations using a two-fluid approach, where the dust is treated as a pressure-less fluid (Paardekooper & Mellema 2006; Woitke 2006a; van Marle et al. 2011; Cox et al. 2012). The dust is linked to the gas through the drag force (Kwok 1975): ${\mathit{f}}_{\mathrm{d}} \mathrm{=} \mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathit{T}}\mathrm{)}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{n}}_{\mathrm{d}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\hspace{0.17em}{\mathit{a}}^{\mathrm{2}}\hspace{0.17em}\mathrm{\Delta }\mathit{v}\hspace{0.17em}\sqrt{\mathrm{(}\mathrm{\Delta }\mathit{v}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{T}}^{\mathrm{2}}}\mathit{,}$(E.1)with n_d the dust particle density, ρ the gas density, a the radius of the dust grain, ${\mathit{v}}_{\mathit{T}}\mathrm{=}\sqrt{\mathrm{3}\mathit{kT}\mathit{/}{\mathit{m}}_{\mathrm{h}}}$ the thermal speed, and Δv the velocity difference between dust and gas. For a given mass, the particle density decreases with a^-3. Therefore, the drag force is weaker for larger dust grains. The sticking coefficient for gas-dust collisions, α_T, is estimated as ${\mathit{\alpha }}_{\mathit{T}} \mathrm{=} \mathrm{0.1} \mathrm{+} \mathrm{0.35}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\sqrt{\mathrm{(}{\mathit{T}}_{\mathrm{g}}\mathrm{+}{\mathit{T}}_{\mathrm{d}}\mathrm{)}\mathit{/}\mathrm{500}}}$(E.2)(Decin et al. 2006), with T_d the dust temperature and T_g the gas temperature. This is an improvement over our earlier simulations (van Marle et al. 2011; Cox et al. 2012), where we kept the sticking coefficient constant. Because the dust temperature is typically much lower than the gas temperature in the shocked CSM-ISM region, we approximate the sticking coefficient by setting T_d    =    0⁸.

Unlike our previous models (van Marle et al. 2011; Cox et al. 2012), which used a cylindrical grid, these simulations are run on a two-dimensional spherical grid with a physical domain of R = 3 pc by θ = 180°. The grid has a minimal resolution of 160 radial by 80 azimuthal grid cells and five additional levels of refinement are allowed, giving the simulations an effective resolution of 2560 ×    1280 grid cells.

E.2. Initial conditions

At the inner radial boundary, the wind parameters are set according to the values in Table 2. Since the star is moving relative to the ambient medium, the outer radial boundary is divided into two zones: (1) for θ < 90° an inflow boundary is set with the ambient gas entering the grid as a parallel stream with v = 28.3 km s^-1; (2) for θ > 90° the gas is allowed to flow out of the grid. This is similar to the approach of Villaver et al. (2012). The gas density in the ambient medium is set at 3 × 10^-24 g cm^-3 (or n_H ≈ 2 cm^-3) and a temperature of 100 K (reflecting the temperature of a cold neutral medium). This temperature also serves as a lower limit for the gas temperature throughout the physical domain. We assume that there is no dust in the ambient medium.

To avoid numerical instabilities we start the simulation without radiative cooling and let it run for 100 000 years. By that time the bow shock has reached its equilibrium position (where the ram pressure of the wind is balanced by the ram pressure of motion of the star through the ambient medium). From this point in time we include radiative cooling and let the simulation run for 200 000 years.

Appendix F: Animations of bow shock simulations

Fig. F.1. Movie of Aori_TISM100K showing gas density and dust particle density for the 200 000  year period from the introduction of radiative cooling to the end of the simulation for the model with a cold ambient medium. (186.8 MB) Access here

Fig. F.2. Similar to Fig. F.1, Aori_TISM8000K showing the 200 000 year period from the introduction of radiative cooling to the end of the simulation for the model with a warm ambient medium. (188.9 MB) Access here

Fig. F.3. Similar to Fig. F.1, the third movie (Aori_vardM - ) shows the 200 000 year period from the introduction of radiative cooling to the end of the simulation for the model with a cold ambient medium and a varying mass-loss rate. (189.8 MB) Access here

Fig. F.4. Similar to Fig. F.2, the fourth movie (Aori_TISM8000K_highvel) shows the 200 000 year period from the introduction of radiative cooling to the end of the simulation for the model with a warm ambient medium and a high stellar velocity. (186.0 MB) Access here

Fig. F.5. Similar to Fig. F.1, the fifth movie (Aori_largegrain) shows the 200 000 year period from the introduction of radiative cooling to the end of the simulation for the model with a cool ambient medium, but for large (100 nm) dust grains. (189.8 MB) Access here

---

Acknowledgments

This research has made use of NASA’s Astrophysics Data System Bibliographic Service and the SIMBAD database, operated at the CDS, Strasbourg, France. PACS was developed by a consortium of institutes led by the MPE (Germany) and including UVIE (Austria); KUL, CSL, IMEC (Belgium); CEA, OAMP (France); MPIA (Germany); IFSI, OAP/AOT, OAA/CAISMI, LENS, SISSA (Italy); IAC (Spain). The development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI (Italy), and CICT/MCT (Spain). SPIRE was developed by a consortium of institutes led by Cardiff Univ. (UK) and include: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). The development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). The development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC (UK); and NASA (USA). NLJC, PR, BV, JB, and MG acknowledge support from the Belgian Federal Science Policy Office via the PRODEX Programme of ESA. FK acknowledges funding by the Austrian Science Fund FWF under project number P23586-N16, RO under project number I163-N16. The Galactic Arecibo L-Band Feed Array H i (GALFA-H i) Survey data set was obtained with the Arecibo L-band Feed Array (ALFA) on the Arecibo 305m telescope. Arecibo Observatory is part of the National Astronomy and Ionosphere Center, which is operated by Cornell University under Cooperative Agreement with the National Science Foundation of the United States of America.

References

All Tables

All Figures