© ESO, 2012

1. Introduction

The chemical enrichment of the Universe plays an important role in explaining stellar and galaxy evolution. Apart from supernovae, asymptotic giant branch (AGB) stars and supergiants have crucial roles in supplying heavy elements into the interstellar medium (ISM). It has been well established that evolved stars lose material during their AGB and supergiant phase via a dusty wind, and this wind material will eventually be dissipated into the ISM (e.g. Tielens et al. 2005). The presence of extended envelopes around AGB stars and supergiants is a direct result of their stellar mass loss. The properties of these (detached) shells are directly affected by the mass-loss history of low- to intermediate-mass stars and thus contain information on this late stage of stellar evolution.

Both the chemical composition and physical conditions of the stellar-wind material returned to the ISM can be altered by several processes, such as mixing of nucleosynthesis products into the envelope, dust formation (Ferrarotti & Gail 2006), shock-induced chemistry in the inner wind envelope (Duari et al. 1999), as well as photo-dissociation due to interstellar cosmic rays and UV photons penetrating the outer envelope (Willacy & Millar 1997; Decin et al. 2010). There is a marked difference between the dust and gas formed in circumstellar environments and present in the ISM (e.g. Jones 2001). To reconcile these differences in the composition of circumstellar and interstellar dust, additional processes are expected to play an important role. One such mechanism involves processing and alteration of dust in shocks that can occur when the stellar wind material interacts with the local ISM. Shocks are not only important in that they can alter the composition of the dusty wind but they also generate turbulence and acoustic noise, which affect the evolution of the ISM (Cox & Smith 1974; McKee & Ostriker 1977) and are believed to be important to the initial phases of star formation (McKee & Cowie 1975; Spitzer 1982).

IRAS observations already hinted at the existence of low-temperature circumstellar dust (van der Veen & Habing 1988). Clear observational evidence of detached shells around evolved stars was found by, e.g., Olofsson et al. (1988), Stencel et al. (1988), and Hawkins (1990). Following these discoveries, infrared observations with IRAS and ISO have revealed a growing number of AGB stars with extended and sometimes clearly detached thermal dust emission (e.g. Young et al., 1993b,a; Waters et al., 1994; Izumiura et al., 1996; Hashimoto & Izumiura, 1998). At the same time, CO radio line emission revealed a number of objects with detached gas shells (Olofsson, 1996, and references therein). Several detached-shell objects were also imaged in the optical (e.g. González Delgado et al., 2001; Maercker et al., 2010; Olofsson et al., 2010). Interferometric radio maps uncovered their detailed structure, which turned out to be remarkable spherical thick CO-line-emitting, geometrically-thin shells (Lindqvist et al. 1999; Olofsson et al. 2000). These observations showed that possibly short phases of intense mass-loss rate are required to form such large detached (spherical) molecular shells of swept-up wind and interstellar material (Rowan-Robinson et al. 1986; Zijlstra & Weinberger 2002; Libert et al. 2007). The enhanced mass-loss rate events could be induced by a thermal pulse (van der Veen & Habing 1988; Vassiliadis & Wood 1993) or, alternatively, the He core-flash at the end of the RGB phase (Dominy 1984).

IRAS also observed for the first time the infrared signatures of bow shocks associated to e.g. ultra-compact H ii regions (van Buren et al. 1990; Mac Low et al. 1991) and runaway OB stars (Noriega-Crespo et al. 1997). Stencel et al. (1988) presented the first IRAS detection of a bow shock around a red supergiant, α Ori. These data were reprocessed by Noriega-Crespo et al. (1997) which revealed a detached bow shock at  ~6′, as well as a linear bar at 9′. Only recently have infrared space telescopes, such as Spitzer and AKARI, revealed much more detail on the infrared signatures of bow shocks around evolved stars, such as R Hya (Ueta et al. 2006), α Ori (Ueta et al. 2008), R Cas (Hashimoto & Izumiura 1998; Ueta et al. 2010), and U Hya (Izumiura et al. 2011). The Herschel Space Observatory now reveals these infrared shells and bow shock regions at the best spatial resolution ever. Detached infrared shells around TT Cyg, U Ant, and AQ And are discussed in Kerschbaum et al. (2010) and on more objects in Kerschbaum et al. (2011). Bow shocks are reported for CW Leo (Ladjal et al. 2010), X Her and TX Psc (Jorissen et al. 2011), and α Ori (Decin et al. in preparation). For o Cet the inner part of the stellar wind bubble bounded and formed by the termination shock is seen (Mayer et al. 2011). For some of these infrared bow shocks, counterparts are detected in UV emission (o Cet: Martin et al. 2007; CW Leo: Sahai & Chronopoulos 2010). These observations indicate that the wind material undergoes additional processing in wind-ISM shocks and that such processing is perhaps more common than previously envisioned.

This paper presents the properties – such as morphology, size, and brightness – of detached or spatially extended far-infrared emission associated to bow shock interaction regions and detached shells around a large sample of AGB stars and red supergiants. The paper is structured as follows. First we discuss the observations, the data processing, and map-reconstruction scheme, and present the infrared images in Sect. 2. In Sect. 3.1 we introduce a simple morphological classification system differentiating observed shapes; in particular, we make a distinction between bow shocks and detached rings. In Sects. 3.2 and 3.3 we present radial profiles, as well as infrared flux measurements of the extended emission. In order to make a qualitative and quantitative comparison of the observed bow shocks and detached rings we review first the basic physics and relevant stellar and interstellar parameters underlying the formation of a bow shock or detached shell (Sect. 4). In Sect. 5 we discuss and collect the relevant parameters in order to predict, for example, the size of the bow shock. Next, in Sect. 6, we present hydrodynamical simulations to elaborate on the origin of the observed morphologies, in particular, the features related to (turbulent) instabilities. The observations are compared to the theoretical predictions and the simulations in Sect. 7. Furthermore, the observations are interpreted in the context of stellar mass-loss and ISM properties as well as with respect to the interaction between these two material flows. Then we discuss the requirements to observe wind-ISM interaction, and thus address the frequent occurrence of bow shocks for our entire sample of AGB stars and red supergiants. The main results and conclusions are summarised in Sect. 8.

2. AGB stars and supergiants in the MESS Herschel Key Programme

The MESS (Mass-loss of Evolved StarS) programme (Groenewegen et al. 2011) is a large far-infrared Herschel (Pilbratt et al. 2010) survey aimed at studying the mass-loss history of evolved stellar objects, from AGB stars to post-AGB, planetary nebulae, luminous blue variables, and supernovae. One key goal of the MESS programme is to resolve the circumstellar envelopes around a representative sample of evolved stars (from AGB stars to PNe), thereby studying the global evolution of mass-loss and circumstellar envelope structure. The details of the programme, the complete target list, data processing approach, and first science results are summarised in Groenewegen et al. (2011).

Briefly, infrared scan maps of a large sample of AGB stars and supergiants are obtained with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) using the “scan map” observing mode with the “medium” (20′′ s^-1) scan speed in the blue (70 μm) and red (160 μm) filter setting. To create a uniform coverage, avoid striping artefacts, and increase redundancy, two observations at orthogonal scan directions are systematically concatenated. Scan lengths range from 6 to 34′, scan leg steps are 77.5′ or 155′, and repetition factors ranged from 2 to 8.

Data processing was performed by applying the standard pipeline steps. In particular, a careful correction of the detector signal drifts is critical to reveal faint extended emission structures with PACS. After deglitching, we applied two different map-making algorithms to assess the quality of the final map, and to verify the faint emission associated to the wind-ISM interaction; these are “PhotProject” with high-pass filter (HIPE v.7.0.) and “Scanamorphos” (Roussel, in prep.). The maps shown here are those obtained with the latter method as these are less susceptible to the filtering of true extended emission. All frames are projected onto an image with a pixel size of 1′′ and 2′′ for the blue (70 μm) and red (160 μm) bands, respectively. Thus, the final image over-samples the instrumental point-spread functions, which have full-width at half-maxima of 5.6′′ and 11.4′′ at these wavelengths, respectively.

Aperture photometry of the central point-like sources observed with PACS at 70 and 160 μm shows that both the “PhotProject” and “Scanamorphos” maps give flux densities that agree within 5%, well within the 15% calibration uncertainties. However, the processing with the updated calibration files in HIPE 7.0 results in 70 and 160 μm fluxes that are 11% and 15% lower, respectively, compared to the values presented in Groenewegen et al. (2011) using HIPE 4.4.

Figures 1 to 4 present the objects in our sample which reveal signs of extended and/or detached dust emission around the central star in either or both the PACS 70 and 160 μm bands. Colour versions of these figures are available in the electronic edition providing enhanced visibility. The different observed morphologies are divided into several classes (see Sect. 3). Out of 78 AGB stars and supergiants imaged with PACS as part of the MESS survey, 32 stars show clear evidence of wind-ISM interaction, 15 stars show detached rings, and 6 show extended irregular emission (Tables 1 and 2). The remaining 30 stars do not show evidence of wind-ISM interaction (Table 3). Tables 1 to 3 provide the basic properties of the observed AGB stars and supergiants: IRAS identifier (Col. 1), Name (Col. 2), Distance (Col. 3), Mass-loss rate (Col. 4), z (Col. 5); height above the Galactic plane, n_H (Col. 6); local ISM density, μ (Col. 7); proper motion, v_⋆ (Col. 8); space velocity, PA (Col. 9); the proper motion position angle measured from north to east, i (Col. 10); the inclination angle of the space motion with respect to the plane of the sky, and v_w (Col. 11) the terminal wind velocity. In addition, we also provide the predicted stand-off distance, R₀ in arcminutes (Col. 12) and parsec (Col. 13), as well as, for Tables 1 and 2, the measured stand-off distance R₀ in arcminutes (Col. 14), parsec (Col. 15), the position angle θ (Col. 16) and the inferred local ISM density, n_H (Col. 17). Columns 18 and 19 give information on the spectral type/circumstellar chemistry and binarity, respectively.

After deconvolution, three stars (R Scl, TX Psc, U Cam) reveal evidence of both a wind-ISM interaction and a detached ring. To measure the ring of TX Psc, a de-convolved 2D-image was used. Detailed studies of some individual (classes of) objects included in the survey here are being presented in companion papers. The complex multiple wind-ISM interaction shells associated to α Ori are presented, together with complementary observations, in Decin et al. (in prep.). The interaction around X Her and TX Psc is described by Jorissen et al. (2011) and that of o Cet by Mayer et al. (2011). First results on the detached ring-like shells around AQ And, TT Cyg and U Ant were already discussed by Kerschbaum et al. (2010), while CW Leo (IRC+10 216) was discussed previously by Ladjal et al. (2010) and Decin et al. (2011). Also, the objects with both inner shells and large detached bow shock regions will be discussed in more detail in a forthcoming paper. This notwithstanding, these sources are included here to complete the sample of observed AGB stars and red supergiants.

3. Observational diagnostics of wind-ISM interaction regions

3.1. Morphological classification of the detached far-infrared emission

From the observations (Figs. 1 to 4) one can immediately recognise different overall shapes of the detected (detached) extended emission. The two most obvious cases are arcs or “fermata” and “rings”. Two additional distinct morphologies are the double opposing arcs or “eye” and the “irregular” emission. The morphological classification is summarised in Table 4 and indicated for all sources in Tables 1 and 2.

“Fermata” (Class I) interaction is characterised by a relatively smooth elliptical arc spanning an azimuthal opening angle  ≥ 120°. This class resembles closest the wind-ISM bow shock shape predicted by theoretical models (Sect. 4). We note that there are different cases within this class reminiscent of the different shapes that occur due to variations in the stellar and interstellar parameters. This will be discussed in detail in Sect. 6. For example, X Pav, EP Aqr and X Her show a peculiar bullet-like shape with back flow emission in the wake of the bow shock. Others show clear signatures of Kelvin-Helmholtz (KH) and/or Rayleigh-Taylor (RT) instabilities (see Sect. 6). Class II (“eyes”) includes objects with two elliptical non-concentric arcs observed at opposing sides of the central source, both have a covering angle of  ≤ 180°. In two cases, VY UMa and U Cam, the two arcs are connected and there is even tentative evidence for a jet structure in the mid plane. Class III consists of spherical structures, i.e. circular “rings”. This category includes typical detached spherical-shell objects, but could also include true rings (i.e. not spheres). From the seven, well-studied CO-detected detached shells around carbon stars (e.g. U Ant, U Cam, TT Cyg, R Scl, S Sct, V644 Sco, DR Ser; Maercker et al. 2010), only the last two are not included in the MESS survey. All known and observed detached CO shell objects show a spatially resolved dust ring, co-spatial with the gas emission. Larger, and thus presumably older, rings are also detected around Y CVn and AQ And. For these objects no corresponding CO shell has been detected, possibly due to the photo-dissociation of CO by the interstellar radiation field (Libert et al. 2007; Kerschbaum et al. 2010). The ring around W Hya is incomplete and deviates slightly from sphericity, and a jet-like structure is also visible in the north-east direction. Sources with diffuse irregular extended emission are classified “irregular” (Class IV). All other targets, for which no evidence of diffuse extended emission has been observed, are assigned to Class X, “non-detection”.

Interestingly, for several objects, such as R Scl, U Cam, and TX Psc, we detect small detached rings in addition to the classical bow shock region further away from the central source. These objects are subsequently assigned to both Class I and Class III (Tables 1 and 2). For these objects it may be the case that we observe a young spherical shell (originating from a wind-wind interaction due to a recent thermal pulse) expanding in a relatively low density medium. In this scenario the local environment of the star has been blown out by the earlier wind/mass-loss which swept out ISM material and created a wind-ISM bow shock. At a later stage a thermal pulse produced a density/temperature enhancement traced by the infrared emission. Alternatively, the small inner ring represents a structure delineating the interface between the (back flowing) termination shock and the free expansion zone. For other objects, such as CW Leo and α Ori, the Herschel/PACS observations reveal irregular multiple incomplete shells in the inner regions of the stellar wind envelope (e.g. Decin et al. 2011).

3.2. Radial distance of arcs and detached shells

For some objects the extended emission is very faint with a low contrast with respect to the sky background. In order to improve the measurement of the angular distance/size of arcs and shells we have constructed radial profiles for Class II and Class III objects. Table 5 gives the radii of the “eyes” (II) and “rings” (III) as identified in the azimuthally averaged radial profiles shown next to the image panels in Figs. 2 and 3 (adopted azimuth angles are given at the top of each radial profile panel; for the “rings” objects it is always 360°. For the “eyes” the radii and radial profiles of both arcs are given. For a few cases with very faint extended emission, the assignment is ambiguous. V Pav possibly belongs to Class III instead of Class II as the dust emission faintly traces a full circle. RT Cap is tentatively assigned to Class III but its ring is not complete and rather consists of two distinct arcs. However, the arcs are concentric, which is not the case for true “eye” objects. For S Cep the ring is also not complete but resembles more a full circle with gaps. The detached shell around Y CVn is very faint and barely visible at a radius of  ~190′′, but in this case the detection is supported by the radial profile and its earlier detection with ISO/ISOPHOT by Izumiura et al. (1996). Tentatively, a faint arc bow shock structure (“fermata”) extending over  ~100° (from position angle  ~30° to 130°) is present at 9′ east of the central star. X TrA is also listed in Class III although also here the ring is faint, with only a brighter arc to the east. Further observations are necessary to resolve these ambiguities.

3.3. Far-infrared dust emission

In Table 5 we give the aperture flux measured in both PACS bands for the elliptical bow shocks and circular rings. To determine the appropriate aperture we fit the extended dust emission with an ellipse. The inner and outer annulus of the aperture are subsequently derived from the azimuthally averaged radial profile. For detached rings the fitted ellipse is circular and the aperture is taken over the entire annulus. Bow shocks have a nearly elliptical shape, and their aperture flux is determined over a limited azimuthal angular range of the annulus. An elliptical annulus was applied in order to measure flux in the non-spherical arcs. The shape of the elliptical annulus was identified from the image with the most significant detection and then directly applied to the other band. The annuli of the detached rings were nearly spherical. For the fermata the flux is only integrated over a limited azimuthal angle covering the observed emission. An illustrative example of a bow shock aperture measurement is shown for θ Aps in Fig. 5.

Cold dust grains emit strongly at mid- to far-infrared wavelengths. Dust emission can, in a simple approximation, be described as a blackbody modified by the grain emissivity; F_ν ∝ B_ν·Q_ν, with Q_ν ≈ ν^β, where β depends on the type of dust considered. For example, β = 1.1 for amorphous carbon in the ISM (Rouleau & Martin 1991) as well as for carbon grains in carbon stars (Jura 1986), while astronomical silicates have β = 2 (Volk & Kwok 1988).

Only for fast shocks (v_shock > 50 km s^-1) are high shock temperatures reached that give rise to UV radiation (e.g. observed for the bow shock of CW Leo) and possibly a strong [O i] 63 μm line. But, as the grains are only weakly coupled to the gas, due to the low densities, T_dust ≪ T_gas. Low velocity shocks on the other hand give primarily rise to dust emission in the infrared. Despite observational efforts with Spitzer and Herschel, no infrared line emission has been detected in AGB star bow shocks yet (Ueta 2011; Decin et al. 2011), and we assume here that the observed emission is entirely due to thermal dust emission (see also Sect. 4).

The observed infrared flux – for optically thin emission – gives a direct measure of the emitting dust mass, and thus – given a gas-to-dust ratio – the total mass. Assuming specific dust properties one can write (Li 2005): ${\mathit{M}}_{\mathrm{dust}}\mathrm{=}\frac{{\mathit{d}}^{\mathrm{2}}{\mathit{F}}_{\mathit{\lambda }}}{\mathit{\kappa }\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}{\mathit{B}}_{\mathit{\lambda }}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{d}}^{\mathrm{2}}{\mathit{F}}_{\mathit{\nu }}}{\mathit{\kappa }\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}}$(1)with F_ν the observed flux (in erg s^-1 cm^-2 Hz^-1), λ the wavelength (in cm), d the distance (in cm), T_dust the dust temperature (in K), B_ν the Planck black body curve in frequency units, and and κ_ν the dust opacity at the observed wavelengths. We adopt a dust opacity of 60 at 70 μm and 10 cm² g^-1 at 160 μm (from Li & Draine 2001, who give κ(λ) = 2.9210⁵(λ/μm)^-2, i.e. β = 2). Table 5 shows that the 70 over 160 μm flux density (F_ν) ratio varies between 1.2 and 3.5, corresponding to a dust temperature of  ~30 K to 40 K (for β = 2). Here we adopt 35 K, which means the dust is slightly warmer than the cirrus dust (T_dust ~ 15 to 20 K; Boulanger 2001). The estimated total dust and gas masses (adopting an average gas-to-dust ratio of 200) emitting in the bow shock region or detached shell are summarised in Table 5. The total observed masses range from  ~1 × 10^-5 to  ~6 × 10^-4   M_⊙. We note that the dust mass inferred from the measured infrared emission is sensitive to the dust temperature, which is not well constrained by the available data, the dust emissivity law (κ(λ)), which is uncertain by order of magnitude and strongly dependent on the chemical composition of the dust, as well as the gas-to-dust ratio.

The total mass of the ambient medium that could be swept up by the stellar wind is roughly the volume set by the stand-off distance of the bow shock or the radius of the detached ring, r, times the mass volume density: M_shell,ISM = 4/3π  r³  ρ_ISM. In Table 5 we give M_ISM using both the ISM densities, n_H, derived using Eq. (5) as well as using n_H inferred from the measured stand-off distance and Eq. (2). For many cases the observed gas and dust mass is less than the potential swept-up ISM mass, but the derived masses are similar in a few cases. The potentially swept-up ISM mass, M_ISM, ranges from 10^-5 to 10^-1 M_⊙. In comparison, assuming constant mass-loss rates, the stellar mass loss after 1000 years amounts to 1 − 50 × 10^-4   M_⊙, again similar to the observed total mass. This consistency is not surprising, since to first order, the total mass accumulated in the bow shock depends on its age (i.e. mass of stellar material piled-up) and the ISM density (mass of swept-up ambient ISM material). However, for bow shocks (spanning a limited azimuthal angle), only a fraction – typically less than half – of the stellar wind and ISM material is entrained into the bow shock region. Furthermore, both ISM and wind material flow from the bow shock apex along the contact discontinuity to be finally shed in the tails. For detached shells, both the ISM and stellar wind mass are presumably preserved, unless either formation or destruction of dust grains and molecules alters the gas-to-dust ratio and its chemical composition (e.g. photo-destruction of CO by the ISRF or processing of grains in the shock).

4. Interaction between stellar winds and the ISM

Bow shocks are common in astrophysical contexts and can occur where two material flows – of different density, velocity, or viscosity – collide. Examples are bow shocks around compact H ii regions due to stellar winds (van Buren et al. 1990; Mac Low et al. 1991; Raga et al. 1995), but also due to winds in binary systems (Stevens et al. 1992; Parkin & Pittard 2008; Pittard 2009; van Marle et al. 2011a) or to movement of an object with a magnetic field through a medium (Baranov et al. 1971; Cordes et al. 1993), between slow and fast winds (Steffen et al. 1998; Steffen & Schönberner 2000), and between stellar winds and the surrounding ISM (Matsuda et al. 1989; Brighenti & D’Ercole 1995; Comeron & Kaper 1998).

In this section we discuss the formation of bow shocks due to interaction of a stellar wind with a low-density medium. First the general case of a star with a stellar wind moving through the ISM is considered. We discuss also the special case of a stationary star with a stellar wind expanding into either an older stellar wind or the interstellar medium.

4.1. The size and shape of a bow shock due to “wind-ISM” interaction

If a supersonic stellar wind and the ambient medium collide, a bow shock interface can be created at a distance from the moving star where the ram pressures and momentum fluxes of the wind and the ISM balance each other. This is called the contact discontinuity. It is a result of interaction between two fluids, but shocks are not necessary. In the absence of shocks, there will be only collisional heating and no Hα emission. For shocks to occur the Mach number M of the fluid velocity v (M = v/v_sound) should be higher than unity. I.e. the relative velocity between the stellar wind and ISM needs to be higher than the speed of sound in the ISM. Note however that the speed of sound in the low density isothermal warm neutral medium (WNM) is very low (v_sound ~ 1 km s^-1) and thus even a slow AGB wind (v_w ~ 10 km s^-1) will move supersonically. Because the sound speed (adiabatic case; γ = 5/3) scales with (P/ρ)^1/2 (or T^1/2)) the sound speed will be lower in the diffuse cold neutral medium (CNM). In the non-adiabatic regime the ram pressure balance is given by ${\mathit{\rho }}_{\mathrm{w}}{\mathit{v}}_{\mathrm{w}}^{\mathrm{2}}\mathrm{=}{\mathit{\rho }}_{\mathrm{ISM}}{\mathit{v}}_{\mathrm{ISM}}^{\mathrm{2}}$. Assuming that the layers mix and that the post-shock cooling is efficient (i.e. instantaneous cooling), the thickness of the dense shell is negligible with respect to the distance from the star. This is, for example, valid for a slow stellar wind interacting with a hot low-density medium (Borkowski et al. 1992). In this approximation the stand-off distance, R₀, defined as the distance between the star and the apex of the contact discontinuity or bow shock, is given by (e.g. Baranov et al. 1971; Dyson 1975; Raga et al. 1995, Wilkin 1996): ${\mathit{R}}_{\mathrm{0}}\mathrm{=}\sqrt{\frac{\mathit{Ṁ}{\mathit{v}}_{\mathrm{w}}}{\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\mathrm{ISM}}{\mathit{v}}_{\mathit{\star }}^{\mathrm{2}}}}$(2)with Ṁ the rate of mass loss and v_w the velocity of the isotropic stellar wind (with respect to the rest-frame of the star), ρ_ISM the mass density of the ambient ISM, v_⋆ the relative space velocity of the star with respect to the ISM. The standard bow shock morphology consists of a forward shock separating the unshocked and shocked ISM, a wind termination shock that separates the free-streaming wind from the shocked wind, and between them a contact discontinuity separating the shocked wind from the shocked ISM (see e.g. Weaver et al. 1977; and Lamers & Cassinelli 1999). If the cooling of the shocked stellar wind is inefficient, a thick hot, low-density gas layer will exist between the free-flowing wind and the bow shock. The termination (reverse) shock and the bow shock, delineating the contact discontinuity, both travel away from this contact discontinuity; the bow shock in forward direction and the termination shock in backward direction with respect to the relative motion of the star in a stationary ISM.

As mentioned above, for efficient cooling, the physical thickness of the shocked region remains small and thus the two shock fronts delineating the contact discontinuity are not resolved. In Sect. 6 we present hydrodynamical simulations of wind-ISM interaction which explore the effect of varying the parameters used in Eq. (2) (Ṁ, v_⋆, n_H) as well as other parameters such as the dust-to-gas ratio and the temperature of the ISM. These simulations also confirm that, despite the formation of turbulent instabilities and inefficient cooling in some cases, Eq. (2) gives in general a rather accurate prediction for the stand-off distance.

Adopting appropriate values for the mass-loss rate, stellar wind velocity, the star’s velocity, and the ISM density (Sect. 5), we can directly predict the stand-off distance, R₀ from Eq. (2). This predicted value for R₀ can then be compared directly to the measured R₀ obtained by de-projecting the observed minimum distance between the star and the bow shock outline. Vice versa, we can use the measured R₀ and Eq. (2) to derive n_ISM. This is important as the local ISM density is difficult to determine observationally (see Sect. 5).

For example, let us apply this to α Ori. The star is at a distance, d = 197 pc (z =  −14 pc), moving at an intermediate velocity of 28.3 km s^-1 through the ISM (Ueta et al. 2008; van Marle et al. 2011b). Its mass-loss rate is 3 × 10^-6   M_⊙ yr^-1 and the stellar wind velocity, v_w = 14.5 km s^-1 (Table 1). From the observed de-projected stand-off distance, R₀ of 5.0′ or 0.3 pc, the local ISM density is then found to be 4.2 cm^-3. This estimate is a factor of two higher than the density derived from the models of the ISM (Eq. (5), n_H = 1.7 cm^-3). The predicted stand-off distance is close to that obtained from a detailed simulation of α Ori (see e.g. van Marle et al. 2011b, and simulation A in Sect. 6).

The 2D-shape of the contact discontinuity (and thus also the termination shock and bow shock) can be solved analytically in the optically thin approximation (Mac Low et al. 1991; Raga et al. 1995, Wilkin 1996): $\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathrm{0}}\frac{\sqrt{\mathrm{3}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\theta }\cot \mathit{\theta }\mathrm{)}}}{\sin \mathit{\theta }}$(3)with θ the latitudinal angle from the apex of the bow shock as seen from the position of the central star, and R₀, the stand-off distance defined in Eq. (2). Alternatively, it can be written as $\mathit{y}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{3}}{\mathit{R}}_{\mathrm{0}}\sqrt{\mathrm{1}\mathrm{-}\mathit{z}\mathit{/}{\mathit{R}}_{\mathrm{0}}}$ (e.g. Mac Low et al. 1991; Fig. 6). Equation (3) is valid only for the case of a relative star-ISM motion in the plane of the sky (i.e. a “side view” of the bow shock or a star-ISM motion inclination of i = 0° with respect to the plane of the sky). In this case, R₀ is directly the angular separation, A, of the star to the observed bow shock outline (Fig. 6). Similarly, when there is no inclination (i = 0°), the angular distance between the star and the surface of the bow shock parabola at an angle of θ = 90° (perpendicular to the apex direction from the star) can be defined as B, where $\mathit{B}\mathrm{=}\mathit{R}\mathrm{\left(}\mathrm{90}°\mathrm{\right)}\mathrm{=}\mathit{y}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathrm{0}}\sqrt{\mathrm{3}}$ (Wilkin 1996 and Fig. 6). In this geometry, R₀ can also be derived from measuring the angular distance between the star and the bow shock interface at an angle of θ = 90° with respect to the star-ISM motion direction (B; as illustrated in Fig. 6), since R(90°) = y(0) = ${\mathit{R}}_{\mathrm{0}}\sqrt{\mathrm{3}}$ (Wilkin 1996).

In general the relative star-ISM motion will be inclined with respect to the plane of the sky, thus i ≠ 0. In order to quantify the effects of an inclined bow shock with respect to the plane of the sky, we simulated the appearance of the three dimensional bow shock surface. This hollow paraboloid wilkinoid surface is defined by rotating the two dimensional shape, i.e. the R(θ) curve, around its axis of symmetry (which is defined by the line connecting the star and the bow shock apex). We numerically construct a 3D wilkinoid surface. These points are then projected onto the plane of the sky under a given inclination angle to compute the observable outline of the bow shock in the optically thin approximation. This outline thus traces the loci of tangential line of sight to the bow shock paraboloid at a given inclination angle. Four such constructed wilkinoids with an R₀ of 1 and star-ISM motion inclinations with respect to the plane of the sky of i = 0,25,50, and 70° are shown in Fig. 6 to illustrate the change in shape with increasing inclination angle. These simulations show indeed that for an increasingly higher inclination i (corresponding to a smaller viewing angle with respect to the line-of-sight), the projected size of the bow shock outline becomes increasingly larger, which is qualitatively consistent with earlier bow shock simulations for compact H ii regions (Mac Low et al. 1991). As noted above, in Fig. 6 we also indicate the two observable quantities, A and B. A is the projected – minimal – distance between the star and the bow shock outline in the direction of relative motion, which can be measured directly from the observed image. B is the distance perpendicular to A, i.e. at θ = 90°, from the star to the bow shock outline. It is important to note here that due to the inclination of the bow shock surface with respect to the line-of-sight, a “pseudo” apex of the observable outline appears where the line-of-sight becomes tangent to the rotated 3D-Wilkin paraboloid. The observable A is in fact the distance from the star to this “pseudo”-apex projected onto the plane of the sky (i.e. not to the real apex of the wilkinoid). At high inclination, but less than 90°, both the “pseudo”-apex and true apex of the observable outline would appear – in the model – as two brightness enhancements. We refer the reader to Figs. 5 and 11 in Mac Low et al. (1991) for an excellent illustration of this effect. However, due to superposition of the true apex and the very bright central star the latter will not be detectable for our stars. In addition, the faint “pseudo”-apex will also be difficult to observe. For i = 90° there is no tangent to the paraboloid, thus only the bow shock apex would theoretically appear as a spherical central peak brightness distribution, where it not for the fact that it will not be visible due to superposition with the very infrared bright central star. For the i = 0° (side view) the simulation gives indeed A = R₀ (i.e. the projected and de-projected stand-off distance are the same) and $\mathit{B}\mathrm{=}{\mathit{R}}_{\mathrm{0}}\sqrt{\mathrm{3}}$ (also here the true apex is observed).

In Fig. 7 the A and B values, in units of R₀, are given as a function of the inclination i. Measuring both A and B would then, in principle, directly constrain both R₀ and i. Unfortunately, measuring both quantities in the detected bow shocks is often not possible due to irregularities and incompleteness of the bow shock. Also the change of outline with varying i is a small effect, thus making the exact inclination angle and stand-off distance highly degenerate. For inclination angles i < 45 ° the change in A or B is less than 10% (Fig. 7 and Mac Low et al. 1991). Only at large inclinations (or small viewing angles) does the observed projected distance A to the bow shock outline differ significantly from R₀, the de-projected stand-off distance. Consequently, using these relations and the kinematical properties of the measured stars and its local ISM, inclination-corrected R₀ values can be derived by measuring A and/or B together with a calculated or assumed i. Thus, first a plausible inclination is derived from the space motion, secondly A is measured from the observed image and thirdly, R₀ is obtained from the R₀-A relation in Fig. 7. Ideally the star’s kinematical properties (proper motion and radial motion) would give the heliocentric 3D space motion vector of the star, wile the observed bow structure would – independently – yield the heliocentric 3D orientation of the bow, whose apex’s 3D orientation represents the heliocentric relative space motion of the star, with respect to the ISM. Thus, the difference of the two would yield the heliocentric 3D ISM flow vector. However, in order to circumvent the high degeneracy between i and R₀ we use the space motion vector inclination angle i to constrain R₀ from the observed bow shock shape, which thus neglects a potential ISM flow (see Sect. 5.1 for a discussion on the local LSR ISM velocity). Table 1 gives the de-projected R₀ obtained via wilkinoid fitting as discussed above for the objects which show emission in the shape of an arc or shell, and thus potentially trace the dust emission outline of a bow shock.

4.2. Special cases of “wind-ISM” and “wind-wind” interaction

In the case of a (nearly) stationary star (with respect to the local ISM and compared to the stellar wind velocity), spherical symmetry of the stellar wind-ISM interaction is preserved (e.g. Libert et al. 2007). Nevertheless, there is a relative velocity difference between the stellar wind and the ambient medium. The interaction of a spherical outflow with external matter leads to the formation of a region of compressed material within two spherical boundaries (Weaver et al. 1977). As before, this region consists of a termination shock (where the freely expanding (supersonic) stellar outflow is abruptly slowed down by compressed material), a contact discontinuity (separating the circumstellar and interstellar matter) and the bow shock (external boundary at which the external medium is compressed by the expanding shell). It is thus a special case of a “wind-ISM” interaction scenario described in the previous section. The total mass of the detached shell would consist of both circumstellar and interstellar matter. The latter can be estimated from the total hydrogen mass that would have been present in a sphere of radius equal to the radius of the detached ring (Sect. 3.3).

However, the fact that the “external” material of the stellar wind is interacting with is moving together with the source (thus keeping the spherical symmetry) may also suggest that in fact the “external” matter is not genuine ISM but rather material remaining from an older mass-loss event, perhaps during the RGB phase. The thermal pulse scenario leads to an interaction between a slow and a fast wind, i.e. “wind-wind” interaction. In this case, v_⋆ from Eq. (2) is the relative velocity between the slow and fast wind. In the rest frame of the star, the second wind is computationally equal to a flowing ISM. A detached dust shell is believed to result from an intense – relatively short – episode of enhanced mass-loss rate and wind velocity (e.g. initiated by a He-shell flash), whose wind has a higher outflow velocity and will thus interact with the previous slower wind. The subsequent sharp drop in mass-loss rate and outflow velocity for a few thousand years directly following the thermal pulse or He-flash will lead to a detached dust shell (e.g. Olofsson et al. 1990; Vassiliadis & Wood 1993; Steffen & Schönberner 2000; Mattsson et al. 2007). The (relative) jump in velocity and mass-loss rate from pre-flash to flash-peak is possibly the most critical parameter governing the formation of a detached shell (Mattsson et al. 2007). The latter scenario seems more appropriate for the formation of geometrically thin detached molecular shells with high expansion velocities showing little interaction with the surrounding medium. Indeed, the detection of bow shocks together with a small detached ring gives further support to this hypothesis. Another possibility could be that the ring represents the (inner) termination shock associated to the bow shock located further outwards. Other scenarios, such as constant mass-loss in combination with non-isotropic mass-loss events and clumpy dust formation have also been detected in the form of non-concentric spherical shells (Decin et al. 2011).

5. Stellar, circumstellar and interstellar properties

The relation (Eq. (2)) between the stand-off distance of the bow shock region and the star’s mass-loss rate, Ṁ, terminal wind velocity, v_w, peculiar velocity (v_⋆), and the ISM density (n_H) is powerful in its simplicity. For observed bow shocks around stars with known (observed) mass-loss properties and space motion, the local ISM density – which is often the most difficult to determine accurately – can be inferred directly from the measured stand-off distance. On the other hand, for stars with known stellar mass-loss, known space motions and (assumed) ISM densities, the stand-off distance can be predicted. In this section we review the relevant and available stellar properties of the observed AGB stars and supergiants as well as the properties of the local ISM. The interplay between these various physical properties determines the final shapes and sizes of the wind-ISM interaction zone, in particular the contact discontinuity as discussed in the Sect. 4.1. Adopting the stellar, circumstellar and interstellar properties presented in this section for each object in our survey, both the predicted ISM densities for observed bow shocks as well as predicted stand-off distances for all targets in the survey are given in Tables 1 to 3.

5.1. Stellar distance and (relative) space velocity

The parallax, proper motion, and the radial local standard of rest (LSR) velocities give a direct estimate of the distance, d, and the absolute peculiar space velocities, v_⋆ (as per Johnson & Soderblom 1987), for the majority of the AGB stars and red supergiants. The LSR velocities have been corrected for the solar motion v_⊙ [(U,V,W) = (11.1, 12.24, 7.25) km s^-1] (Schönrich et al. 2010)¹. Parallaxes from re-processed Hipparcos data (van Leeuwen 2007) were taken, when available, if the uncertainties are less than 20%. For several (distant) stars, the parallactic distances derived from Hipparcos data have large uncertainties or have not been obtained at all. For these targets distances can be derived, for example, from the (pulsation) period-luminosity(K) relation, and observed apparent K magnitudes. A generic Galactic P-K relation valid for both oxygen- and carbon-rich (Mira) variables was established by e.g. Whitelock et al. (2000, 2008): M_K =  −3.51 [log   P − 2.38]  − 7.25. Vlemmings et al. (2003) obtained VLBI parallaxes for four Miras that were more consistent with distances derived via the P-K relation than those obtained with Hipparcos. Tables 1 to 3 give the derived peculiar velocities and peculiar absolute proper motions with their associated position angles.

Radial (LSR) velocities are taken from CO line surveys where available. Primary sources are De Beck et al. (2010) and Menzies et al. (2006), which have good agreement for sources in common to both studies (typically the differences are less than 3 km s^-1).

At this point we assume that the relative peculiar velocity between the ISM and the star, v_⋆, is determined entirely by the star space velocity with respect to the local standard of rest (LSR) (i.e. the stationary ISM). In other words, we assume that there is no flow of the ISM itself. This could be a simplification for some cases where the ISM may have an appreciable flow velocity as matter is being blown away by super-bubbles (see Ueta et al. 2009a,b for the case of α Ori).

To test this, we estimate the local LSR ISM velocity (v_r,v_l,v_b) from the galactic rotation (Oort 1928, Feast & Whitelock 1997) at the position of the stars and correct v_⋆ to obtain the local v_{⋆ − ISM}. For the majority of the objects the v_{⋆  − ISM} is within  ± 20% of v_⋆, thus decreasing/increasing the predicted stand-off distance by the same factor. Only in a few cases does this correction give a significant change in peculiar velocity; for U Cam, V Hya, T Lyr, and μ Cep this gives relative space velocities that are approximately a factor of 1.5 higher, and for Y Pav a factor of 2 lower.

The average space velocity, v_⋆, for our sample of stars is  ~36 km s^-1 (Tables 1 to 3), consistent with the average velocity of  ~30 km s^-1 for Galactic AGB stars (Feast & Whitelock 2000).

5.2. Mass-loss rate

The mass-loss rate of evolved stars is important for the enrichment of the ISM. It also directly affects the shape and size of the stellar wind-ISM interaction (Sects. 4 and 6). Mass-loss rates can be estimated directly via modelling of observed CO line profiles (e.g. Knapp et al. 1998; Groenewegen & de Jong 1998; Groenewegen et al. 2002; De Beck et al. 2010). In general the derived (gas) mass-loss rates depend on the line strength, terminal gas expansion velocity (${\mathit{v}}_{\mathrm{w}}^{\mathrm{2}}$), and distance (d²). The average gas mass-loss rate from  ~300 Galactic carbon stars is 1.1 × 10^-5   M_⊙ yr^-1 (Groenewegen et al. 2002). The latter should probably be reduced by an order of magnitude since new Hipparcos results show that the distances adopted by Groenewegen et al. (2002) are over-estimated by a factor of two to four. Using CO-derived mass-loss rates, an empirical relation has been established between the mass-loss rate and luminosity variability period (e.g. De Beck et al. 2010): log(Ṁ) = –7.37 + 3.42  ×  10^-3 P (for P ≤ 800 days) and log(Ṁ) = –4.46 (for P ≥ 850 days), with Ṁ in units of M_⊙ yr^-1. However, the scatter on this relation is rather large with a typical uncertainty of a factor of 10 (see De Beck et al. 2010, for details). Similar results are obtained for carbon Mira variables by Groenewegen & de Jong (1998): log(Ṁ) = 4.08 log P – 16.54. Adopted mass-loss rates – based on CO observations or luminosity period – are included in Tables 1 to 3. If neither Ṁ nor period are available, we adopt a generic value of 5 × 10^-7   M_⊙ yr^-1. Uncertainties in Ṁ arise predominantly from inaccurate distance determinations for distant stars.

5.3. Stellar wind velocities

For the (terminal) wind velocity, v_w, we adopt the terminal velocity of the CO envelope (e.g. De Beck et al. 2010; Tables 1 to 3 and references therein). Groenewegen et al. (2002) find an average expansion velocity of 18.7 ± 6.1 km s^-1 for a sample of 330 carbon stars. For a set of 24 oxygen-rich and 13 carbon-rich AGB stars, De Beck et al. (2010) find mean v_w values of 14.5 and 15.4 km s^-1, respectively. These values are 4–10 km s^-1 higher than those of Ramstedt et al. (2006) for a set of 77 oxygen-rich and 61 carbon-rich stars. For carbon-rich stars, this average is  ~3 km s^-1 lower than the one found by Groenewegen et al. (2002). If no CO terminal velocities are available, we adopt a generic value of v_w = 15 km s^-1, for both oxygen- and carbon-rich stars.

5.4. Local interstellar medium density

The ambient (uniform) mass density of the ISM is defined as ρ_ISM = μ_H  m_H  n_H, with μ_H = 1.4, the mean nucleus number per hydrogen atom for the local medium, m_H is the mass of the hydrogen nucleus, and n_H the interstellar hydrogen nucleus density. Although the ISM is far from uniform, there is a general dependence between the interstellar matter (space-averaged) volume density, n_ISM, and z, the distance from the Galactic plane².

The different phases of the ISM have different average densities, scale-heights and filling factors (see e.g. Savage 1995; Boulanger 2001; Ferrière 2001; Wooden et al. 2004; Cox 2005). Dickey & Lockman (1990) used radio observations of H i to derive the structure of the CNM, WNM and WIM (warm ionised medium):. ${\mathit{n}}_{\mathrm{H}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{0.39}{\mathrm{e}}^{\mathrm{-}{\left(\frac{\mathit{z}}{\mathrm{127}\mathrm{pc}}\right)}^{\mathrm{2}}}\mathrm{+}\mathrm{0.11}{\mathrm{e}}^{\mathrm{-}{\left(\frac{\mathit{z}}{\mathrm{318}\mathrm{pc}}\right)}^{\mathrm{2}}}\mathrm{+}\mathrm{0.06}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{\left|}\mathit{z}\mathrm{\right|}}{\mathrm{403}\mathrm{pc}}}\mathit{.}$(4)

Based on UV observations of H i and H₂, Savage et al. (1977) and Diplas & Savage (1994) arrive at a similar result for the WNM as Dickey & Lockman (1990). An alternative relation is given by Loup et al. (1993) which is based on Spitzer (1978) (density n_H(z = 0) = 2 cm^-3) and Mihalas & Binney (1981) (scale height of 100 pc): ${\mathit{n}}_{\mathrm{H}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{2.0}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{\left|}\mathit{z}\mathrm{\right|}}{\mathrm{100}\mathrm{pc}}}\mathit{.}$(5)

These two relations are shown further below in Fig. 11 (Sect. 7). Thus, we obtain some first insights into the local ISM density for our objects (Tables 1 to 3). Evidently, one should be cautious using these relations for specific cases as there is structure in the ISM on all spatial scales, with different phases that have different filling factors. For example, diffuse to molecular clouds (n_H ~ 10 − 1000 cm^-3) have an order of magnitude lower filling factor than the low density WNM and WIM surrounding it. Thus, statistically about 10% of the objects in our survey could be moving through a denser medium than inferred from the above relation for the WNM. Generally, these relations do indicate that n_H drops as one moves away from the Galactic plane, leading to – on average – lower volume densities of the ISM. We use Eq. (5) to get a first estimate of the local ISM density for the objects in this survey (Tables 1 to 3). Accurate measurements of the local ISM density for all stars in our sample would be desirable, but are currently – to the best of our knowledge – unavailable. In fact, we will show that for certain cases the observations of bow shocks can potentially be used, via the measurements of the stand-off distances, to derive estimates of the local ISM density.

5.5. Circumstellar chemistry, spectral type and binarity

In addition to the above stellar properties that have a direct impact on the theoretical stand-off distance via Eq. (2), other possibly pertinent information such as spectral types of the central star and dominant circumstellar chemistry (oxygen versus carbon rich), and binarity is also included in Tables 1 to 3. If and how these stellar attributes could affect the occurrence and shaping of bow shocks is discussed further in Sect. 7.

6. Hydrodynamical models of interaction between the slow stellar wind of a moving star and the ISM

Hydrodynamical simulations offer the opportunity to explore the effect of varying physical properties of either the stars (e.g. mass-loss rate, wind and (relative) space velocity) and/or the ISM (density, temperature) in a coherent systematic way. For example, Steffen & Schönberner (2000) and Libert et al. (2007) used hydrodynamical simulations to show that a brief episode of increased mass-loss rate could give rise to an expanding, geometrically thin shell. Wareing et al. (2006b) applied numerical simulations of a two-wind model to explain the observed structure around R Hya. Simulations by Wareing et al. (2007a) and Wareing et al. (2007b) indicate that a higher mass-loss rate (for similar ISM density and space velocity) will result in more pronounced Kelvin-Helmholtz instabilities (see their Fig. 1). Wareing et al. (2006a) obtain bullet-shaped emission structures simulating a PN (their Figs. 3 and 5). Villaver et al. (2003) include time-dependent mass-loss in simulations of wind-ISM interaction, leading to time-dependent stand-off distances.

One particular strength of hydrodynamical simulations is that they allow studying the formation, growth and dissipation of fluid instabilities – like Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instabilities – important in wind-ISM shock interactions. RT instabilities occur when a dense, heavy fluid is accelerated by a light fluid. Normally two plane-parallel fluid layers are meta-stable, but a small perturbation can destroy this delicate equilibrium. This manifests itself as so-called inter-penetrating “RT-fingers”. Such instabilities could be quenched by a restoring force such as a magnetic field, thus preventing these instabilities to grow (Chandrasekhar, 1961). The KH instability results from a velocity shear between two fluid layers. This flow of one fluid over another will induce a centrifugal force which leads to changes in pressure which amplifies the ripple. Together with a RT instability, the KH instability will form structures in the shape of mushroom caps on the end of the RT fingers. The KH time-dependent turbulent eddies can form a complex structure arising in a steady flow, if they become large enough to influence that large-scale morphology of the shocked gas. Other instabilities that can occur are the non-linear thin shell instability (Vishniac 1994) and the transverse acceleration instability (Dgani et al. 1996) as shown numerically by e.g. Blondin & Koerwer (1998) and Comeron & Kaper (1998).

6.1. Hydrodynamical simulations

Here we present a series of seven simulations of the interaction between the ISM and the circumstellar medium of moving, evolved stars in order to find out how the morphology of the bow shock varies with the various stellar wind and ISM parameters. The different pertinent parameters for the simulations are summarised in Table 6. For our hydrodynamical simulations we use the MPI-AMRVAC code (Keppens et al., 2011). This code solves the conservation equations of hydrodynamics on an adaptive mesh (AMR) grid. We use a 2-D cylindrical grid in the r-z plane. The basic resolution is set at 80 grid points per parsec, but allows four additional levels of refinement, each doubling the effective resolution. This gives us a maximum effective resolution of 1280 grid points per parsec. Since some of the models require a larger physical space we use grids of different sizes for those simulations. In those cases we increase the number of grid points on the basic level to maintain the same resolution. The expanding stellar wind is inserted by filling a small circle with wind material. The motion of the star is handled by giving the ambient medium around the star a velocity in the z-direction. Therefore, we are simulating the wind interaction in the frame-of-reference of the star. Since all models are 2-D and the simulated space lies along the direction of motion of the star, the snapshots (Figs. 8 through 10) show projections for i = 0.

To the basic conservation equations we have added the effect of optically thin radiative cooling. This is necessary since some of the shocks are strongly radiative, which changes the morphology of the shocked gas. In the case of simulation G, where the ISM has a high temperature, we put a lower limit on the cooling in the ISM, so that in the ISM the temperature does not fall below 8000 K. More importantly, we have also added a (simplified) dust component (van Marle et al. 2011b). This is done by using a two-fluid approximation, with the dust represented as a gas without internal pressure. For simplicity we only include dust in the wind and neglect the dust component of the ISM. This gives us a set of five partial differential equations (not counting the vector components) for the pressure balance for the gas and dust. For the gas we have conservation of mass, momentum and energy, while for the dust we only have conservation of mass and momentum. The appropriate equations are given in van Marle et al. (2011b). Radiative losses depend on the hydrogen and electron particle densities (derived from ρ assuming full ionisation with hydrogen mass and involve a temperature dependent cooling curve Λ(T). This cooling curve has been calculated with the CLOUDY code (Ferland et al., 1998) and includes radiative losses through IR radiation from the dust. We also include a drag force linking the gas and dust, which is derived from a combination of Epstein’s drag law for the subsonic regime and Stokes’ drag law for the supersonic regime (Kwok 1975; see also van Marle et al. 2011b). The drag force depends on the dust particle density and radius (0.005 μm), the velocity difference between gas and dust, as well as the thermal speed of the gas (Kwok 1975). For further details we refer the reader to van Marle et al. (2011b).

Since our equations are purely hydrodynamical, we do not take into account magnetic fields or the effect of radiation on the dust particles. We also neglect destruction and creation of dust particles. Radiative processes are not included and thus the dust temperature cannot be treated appropriately (currently only collisional heating of dust could be accounted for). Also the gas is purely heated by collisions, and cooled radiatively; photo-ionisation is not included. Introducing radiative cooling by dust will cause the bow shock to become thinner and more unstable.

Adding destruction processes such as collisional heating or UV irradiation could, in principle, be included, but at a high cost in the required computing time. Both mechanisms are not expected to be very efficient as grains are not easily destroyed, because they are very effective radiators (emitting all heat immediately in the infrared) and the shocks discussed here are not very strong, and the radiation field would not be strong enough to destroy grains. Furthermore, dust can also be created in the high density shock regions, thus lowering the effective dust destruction rate. Further details are given in van Marle et al. (2011b).

The input and grid parameters are summarised in Table 6. We start with a basic model (simulation A), which has input parameters based on the observations of α Ori (Ueta et al., 2008). Using this model as a starting point, we vary individual parameters to investigate the effect on the morphology of the bow shock. Note that for a stationary star, v_⋆ = 0 km s^-1, the shock between the wind and ISM will drive a spherical shell of wind material sweeping up the ISM material, thus leading to detached shells/rings (not shown).

6.2. The morphology of the bow shock

As shown in Figs. 8 through 10, the collision between the expanding wind and the moving (in the rest-frame of the star) ISM creates a bow shock. The location and morphology of the shocked gas depend on the exact input parameters. These figures show the gas density and dust particle density in the circumstellar environment after the bow shock has reached its (semi-) permanent equilibrium distance from the star. Typically, this takes between 50 000 and 150 000 years, depending on the parameters of the simulation. Our basic model (Simulation A, Fig. 8) shows a standard bow shock morphology: a forward shock separating the unshocked and shocked ISM, a wind termination shock that separates the free-streaming wind from the shocked wind, and between them a contact discontinuity separating the shocked wind from the shocked ISM, with the temperatures in the shocked wind typically lower than in the shocked ISM due to the higher density of the former. This pattern is repeated in all subsequent models, with the exception of the fast moving star (Simulation E, Fig. 10) and the high ISM temperature cases (simulation G, Fig. 10), both of which show a chaotic bow shock with strong turbulent behaviour.

Behind the star a relatively structureless – empty – region remains. This is only filled by the free-streaming wind of the star. Note that in our simulation the wind velocity is always less than the stellar velocity. Exceptions occur if the shocked gas region has a high thermal pressure, which pushes it into the empty region (Simulation C, Fig. 9 and simulation G, Fig. 10), or if the shocked gas is very turbulent (Simulation E, Fig. 10).

The location of the shock can be approximated by the wind and ISM conditions as described in Eq. (2). In Table 6 we give the location of the contact discontinuity (R₀) as well as the location of the forward shock of the shocked-gas region (R_s), where the ISM transition from unshocked to shocked gas is located. In some cases R₀ has to be approximated since instabilities can make it a time-dependent value. For simulation A (the basic model), the analytical approximation puts the contact discontinuity at 0.42 pc. In our simulation the position of the contact discontinuity deviates from this prediction, lying at about 0.26 pc (while the forward shock lies further ahead at  ~0.35 pc). The difference can be explained as a result of radiative energy loss in the shock. The scaling, as described by Eq. (2), appears to be quite accurate. For example, for an order of magnitude decrease in ISM density (simulation F, Fig. 10) or increase/decrease in mass-loss rate (simulation B, Fig. 8), the location of the bow shock moves away/closer with a factor of about $\sqrt{\mathrm{10}}$.

The dust grains, which start with the same velocity as the wind, are initially carried along. When the wind reaches the termination shock, the gas slows down abruptly. Because they are not subject to the shock, the dust grains keep moving, but are now subject to drag force due to the difference in velocity with the gas. This causes them to slow down over time. As a result they tend to pile up at the contact discontinuity (van Marle et al. 2011b). Once the grains reach the contact discontinuity, they tend to follow the local instabilities and are eventually carried downstream. Those grains, which crossed the contact discontinuity, will usually go downstream faster, since in most of our simulations the shocked ISM has a higher velocity than the shocked wind. The influence of the grain size on the gas-dust interaction is investigated in van Marle et al. (2011b).

6.3. Instabilities in the bow shock

In most of our simulations, instabilities occur along the contact discontinuity. These instabilities consist of a combination of RT and KH effects (see above). In our case, the (relatively) low-density-shocked ISM exerts a force on the denser shocked wind, which leads to the formation of RT fingers. As a result, the wind material starts to flow into the shocked ISM. Here the wind material is subject to a sheer-force due to the relative motion of the star with respect to the ISM. These instabilities are enhanced by the dust (albeit a secondary effect), which tends to continue in a straight line in an attempt to cross from the shocked wind into the shocked ISM, dragging the gas along with it. This is best observed in the case of simulation B (Fig. 8). The velocity difference between the shocked gas layers causes KH instabilities to develop wherever the RT instabilities cause a local displacement of gas. These lead to the characteristic “cyclonic” features (see Figs. 8–10). The presence of dust slows down the formation of KH instabilities, since the circular motion of the KH instabilities has difficulty overcoming the inertia of the dust grains. Therefore, they are most visible in simulation D (Fig. 9). Also, the KH instabilities can form more easily if the low density region (the shocked ISM) is extended, so that the instabilities do not hit the forward shock (simulation G, Fig. 10).

In the case of the star with a low mass-loss rate (simulation C, Fig. 9), such instabilities are absent. The low mass-loss rate leads to a reduced density in the shocked wind. As a result, the density difference between the shocked wind and the shocked ISM is small, which reduces the RT effect. Also, due to the low density, the radiative cooling, which scales with the density squared, is less effective (van Marle & Keppens, 2011). This is important because radiative cooling tends to favour the formation of small, high density clumps, which in turn can serve as a start of other instabilities. Because the shocks are nearly adiabatic, the shocked-gas region is wide and its thermal pressure pushes material into the area behind the star, rather than leaving it empty.

The fast-moving star (simulation E, Fig. 10) has a far more irregularly shaped bow shock. In this case the shocks on both sides are highly radiative, leading to a very compressed shocked gas region. This, combined with the strong ram-pressure from the ISM, causes the entire shell to become unstable, leading to a ragged form which will change considerable over time as local instabilities grow to a size where they dominate the entire structure of the shell. This turbulent motion also causes the gas (and the dust grains) to end up behind the star.

7. Discussion

7.1. Comparison between observed and hydrodynamical bow shock morphology

The simulations presented in Sect. 6 only cover part of the parameter space spanned by all stellar, circumstellar and interstellar properties. In particular, the simulations focus on the effects of increasing the star’s peculiar velocity with respect to the ISM, changing mass-loss rate, increasing the ISM temperature and lowering the ISM density.

One evident issue revealed by the simulation is that the shocked gas region can have considerable spatial extent (this occurs when the shocks are adiabatic – as in most cases presented here – as opposed to radiative). In the approximation of efficient cooling, the shock region will be thin and unresolved; however, the simulations show that this is not necessarily valid for all cases and can lead – depending also on the exact cooling law – to extended interaction zones. Depending on the exact morphology and the temperatures in the shocked gas region, one could – theoretically – be observing either the forward shock, reverse shock or contact discontinuity. However, the current far-infrared observations are most sensitive to cold dust grains. Due to the limited spatial resolution, it is not clear which region is actually represented by the far-infrared emission. It may represent the contact discontinuity (i.e. the entire unresolved shocked region), the bow shock, or the termination shock (Libert et al. 2007). Our simulations in Sect. 6 show that most dust grains pile-up at the contact discontinuity where they tend to follow the local instabilities which carry the grains downstream. Thus, the simulations suggest that the observed far-infrared dust emission primarily traces the contact discontinuity. On the other hand, if grains are destroyed in the shock they will be found in the unshocked stellar wind region and thus far-infrared dust emission would delineate the termination shock (and the contact discontinuity if the shock region is physically thin enough). This might prove useful in explaining the spatial offset between observed UV emission and far-infrared emission in the bow shock of CW Leo (Ladjal et al. 2010). Further work on the destruction and formation of dust grains (c.q. alteration of the dust size distribution) in shocks is warranted. As stipulated above, the simulations do not include radiative transfer and can thus not provide appropriate dust temperatures needed to simulate the infrared emission.

RT and KH instabilities

The presented PACS infrared observations reveal RT and KH instabilities in astrophysical bow shocks. RT fingers can be seen in the shock region of R Scl. KH “wiggles” or density knots are more common, and can be seen in the bow shock regions of, for example, UU Aur, R Hya, X Pav, EP Aqr, μ Cep, R Leo, RT Vir, X Her, and V1943 Sgr.

Almost spherical bow shocks spanning a large azimuthal angular range correspond closest to simulations A, B, and F, the nominal, high Ṁ, and low n_H cases, respectively. For low n_H (simulation F) the KH instabilities are more pronounced, while the RT instabilities are more evident in the high Ṁ case. The smooth, spherical bow shocks of α Ori, UU Aur, and X Pav are similar to the morphology of the simulation A (nominal case; Fig. 8). UU Aur and X Pav show back flow emission also seen in the simulation. The “wiggles” in the bow shocks of R Hya and R Leo correspond to those seen in simulation F (Fig. 10), which contradicts the high ISM density derived for R Leo (see Sect. 7.5). Both X Her and V1943 Sgr have high space velocities leading to morphologies much like that seen in simulation E (Fig. 10).

Further comparisons are hampered on the one hand by the (unknown) projection effects in the observations and on the other hand by the, currently, limited parameter space covered with the few simulations at hand. A more extensive, multi-parameter grid of simulations and/or a detailed modelling of individual cases will be required to match observations with simulations to confirm whether or not the latter accurately predict the former.

7.2. Proper motion and inclination

In most cases the proper motion vector (if known) points roughly (within  ~20°–30°) in the direction of the observed (Class I) bow shock apex. Exceptions are μ Cep and o Cet. One possible explanation is that, like o Cet, μ Cep is also a long-period (tens of years) binary, for which the present-day proper motion is not pointing in the same direction as it was when the bow-shock matter was expelled. This conjecture needs to be further tested by further observation and modelling of these systems and their observed bow shocks.

For the majority of the Class II objects, the proper motion direction is aligned roughly along the symmetry axis (perpendicular to the arc-star-arc line). Only for V Pav does the proper motion point to the north-west arc. For Class IV the proper motion direction coincides with the location of extended irregular emission. The proper motion for Class III objects is not consistently aligned with, for example, shell features. In a few cases (TT  Cyg, U Ant) there is some stronger shell emission in the direction of the proper motion, but this is not observed for the other detected rings.

There is no clear indication of any dependency on the inclination angle, i. All five classes include sources with a range of inclinations between almost 0° up to 60–70°. Except for VY UMa (i =  −1°), there are no objects with small inclination angles among the “eyes” (all have i ≥ 28°), whereas the “fermata”, “rings”, and “irregular” all have a more uniform distribution over inclination angles. However, we stress that at this point the sample is too small to derive any further conclusions.

For about one third of the “non-detection” class, the proper motions are not known and thus their space velocities are, at best, lower limits based on radial velocities. This is in part a distance bias. The objects with non-detections are on average further away and consequently their proper motions have been too small to be measured. This means also that the predicted stand-off distances are upper limits.

7.3. Binary interaction, circumstellar chemistry, and magnetic fields?

Binary interaction

There is little evidence of any binary interaction in our sample. Even though for several stars in our sample, there is solid evidence of close companions (see Tables 1 to 3), none show suspicious features in their bow shock interfaces. However, many of these appear in the category “irregular”, which could indicate that such wind-wind-ISM interaction is only visible on small spatial scales. It is difficult to arrive at any firm conclusion as the number of binaries in the survey is low. Excluding the optical “visual” binaries, as these are not yet confirmed, out of the 25 objects in the “fermata” class only one is a known binary (Mira) which indeed reveals a peculiar shape of the interaction zone (Mayer et al. 2011). Three others, W Aql, EP Aqr, and θ Aps are solid candidates for binarity. None of the “rings” are binaries. Noteworthy, three out of six of the “irregulars” are confirmed binaries, and two of the six “eyes” are potentially binaries. Finally, also six out of 30 non-detections are binary systems. Including visual binaries adds three binaries to each of the “fermata”, “eyes”, and “non-detection” classes, and four binaries to the “rings”. In this case, five out of seven “eyes” show evidence of binarity, whether this means there is a connection between binarity and the “eyes” morphology can not be confirmed nor excluded at this point.

Oxygen-rich versus carbon-rich chemistry

One key property of the AGB stars and red supergiants in our sample is their circumstellar envelope chemistry. This can be either carbon-rich (the dust is predominantly made up of amorphous carbon, graphite and silicon carbide) or oxygen-rich (the dust is composed of silicates and oxides). Although we do not find a strong dependence on the presence versus absence of bow shocks and rings with respect to chemistry, there does appear to be a distinction between the shape of the extended emission found around carbon and oxygen-rich objects. The “fermata” and “irregular” classes include a high fraction, 22 out of 31, of O-rich stars. The “eyes” and “rings”, on the other hand, include only few, three out of 22, O-rich stars. It seems that at least the “ring” structure is typical of C-rich chemistry, as already pointed out by Olofsson et al. (2000). A further argument favouring such a link between “ring” structure and C-richness is that all “fermata” stars also displaying a “rings” structure are C-rich stars (noted with an asterisk in Table 1, and are also listed in Table 3).

Could the circumstellar envelope structure be related to chemistry through differences in mass loss rates? It has been shown (Knapp 1985) that oxygen-rich stars have a wider range of mass-loss rates than C-rich stars for which the Ṁ is globally more homogeneous. On average, mass loss rates of C stars are higher than those of MS and S stars (Guandalini 2010). Apart from the fact that the mass-loss rate standard deviation of O-rich stars is indeed found to be somewhat larger than that of C-rich stars, no clear dichotomy emerges between the mass loss rates of O-rich stars and C-rich stars in the present sample.

Furthermore, the higher emissivity (Draine & Lee 1984; Li & Draine 2001) and lower gas-to-dust ratio of 160 for oxygen-rich dust compared to the higher ratio of 400 for carbon-rich dust (see e.g. Knapp 1985; Heras & Hony 2005) will, for O-rich dust, give rise to brighter infrared emission for a given total dust and gas mass. One could then argue that the O-rich stars mass loss is more easily detected even with a patchy structure. This could perhaps explain why O-rich stars are detected preferentially among both “irregulars” and “fermata” classes.

An alternative explanation is that the “Ring” (and possible the “Eyes”) morphology, are explained through (i) the interaction between a fast wind sweeping out matter from a previous slower wind, or (ii) a phase of drastically enhanced mass loss, e.g. caused by a thermal pulse. If this explanation holds, and if one assumes that most thermal pulses events are followed by a third dredge-up (TDUP) episode, then the discriminant character would not be O-rich or C-rich, but instead pre-TDUP or post-TDUP (i.e., thermally-pulsing (TP) AGB) stars). All intrinsic S-type stars, C-type stars and technetium-rich M stars in Tables 1 to 3 are TP-AGB stars. In fact, all “Ring” stars are TP-AGB stars.

Stellar and interstellar magnetic fields

Some wind-ISM models include the effects of an interstellar magnetic field on the shaping of the bow shock region. Heiligman (1980) predicted that the ISM magnetic pressure could deform the outer shells of very large spherical halos of PNe into “lemon” shapes, with the axis of symmetry inclined to the ordered ISM magnetic field direction. Soker & Dgani (1997) and Dgani (1998) elaborated on this scenario predicting different bow shock shapes depending on the relative magnitudes of the star’s space velocity, v_⋆, the wind velocity, v_w, and the Alvén speed, v_A. The double arc objects (Class II) such as VY UMa and AQ Sgr could perhaps provide evidence for this scenario, although both are at relatively high latitudes and their symmetry axes are not aligned with the Galactic plane (assuming that at these locations the magnetic field is also parallel to the Galactic plane). It is also possible that the star has its own magnetic field and thus may exhibit axis-symmetric mass loss. Preliminary simulations (Van Marle, priv. comm.) indicate that double arcs as well as jets (tentatively identified for VY UMa and W Hya) could be formed by a (rotating) star with a magnetic field.

7.4. Detached spherical shells

The observed detached shell objects (Class III) show quite a range of expansion ages calculated out of their measured or assumed expansion velocities, distances and angular sizes. Ranging from the youngest shells like U Cam or R Scl with ages of the order of 1000 years and much older ones like AQ And, Y CVn or UX Dra with ages of a few 10 000 years we see also a significant trend in the detectability in fossil mm-CO. Whereas all objects with shell ages below some 1000 years are prominent in CO, no detached shell older than 10 000 years was successfully detected yet. This could be understood in terms of photo-dissociation by the interstellar radiation field in the older extended shells with sizes of 0.2 pc or more (Kerschbaum et al. 2010). Arimatsu et al. (2011) find that the detached shell of U Ant detected with PACS is not the outermost one found with AKARI. Thus, the shell detected with PACS is either a density enhancement due to a two-wind interaction or it could represent the termination shock, similar to the scenario proposed for Y CVn by Libert et al. (2007).

Another interesting finding with respect to Class III objects are the objects that have both a detached shell and a bow shock interaction region, with R Scl being the showcase object. Seeing the bow shock far out off the detached shell more or less excludes one of the explanations of detached shells, namely their interpretation as being the wind-ISM interface (Libert et al. 2007).

7.5. Comparison between predicted and observed R₀: Implications for the ISM density

For a number of objects in our sample, previous studies have reported on the detection of wind-ISM bow shocks with IRAS, Spitzer or AKARI. The measured deprojected angular stand-off distance obtained with PACS agree very well with earlier results (Table 7). The predicted values for R₀ (as function of $\sqrt{{\mathit{n}}_{\mathrm{H}}}$) give ISM densities that are somewhat higher than those computed from Eq. 5. For α Ori, Ueta et al. (2008) predict n_H = 1.5 − 1.9 cm^-3 which is indeed in line with n_H = 1.9 cm^-3 derived from Eq. (5), but higher than that inferred from the stand-off distance, n_H = 4.6 cm^-3 (Table 7). Ladjal et al. (2010) derive n_H ≥ 2 cm^-3 for CW Leo, which is a factor of 2 higher than derived from both our predicted stand-off distance as well as from Eq. 5 (n_H = 0.5–1.0 cm^-3). For o Cet the observed features delineate most likely the termination shock and not the contact discontinuity, this measure is thus a lower limit to R₀, which will be larger than 1.2′(e.g. more like 2′ as observed in the UV) leading to a lower density of n_H = 0.2 cm^-3.

Adopting the values for the different parameters given in Tables 1 to 3 yields R₀ less than 8′ for all objects with detected bow shock interaction. For all detected wind-ISM interaction objects, the predicted stand-off distances range from about 0.02 to 0.9 pc, in line with the observed de-projected values. This gives some credibility in the predicted stand-off distances (from Eq. (2)) for the objects in Class X (“non-detections”). These results suggest that additional bow shocks should have been detectable for other stars within 300 pc, in particular for R Dor, Y Lyn, RS Cnc, and possibly HD 100764 and RX Boo. We note that higher sensitivity observations could possibly still reveal the expected extended emission.

The uncertainty of the predicted stand-off distance depends on the different stellar and ISM parameters. The space velocity, v_⋆ (for nearby stars), as well as the wind velocity, v_w are well determined observationally with small (~10%) uncertainties. The uncertainties in the mass-loss rates can be an order of magnitude, in particular for the targets with inaccurate distance estimates (Sect. 5) which introduces errors of a factor of three in R₀. Finally, the local ISM density is difficult to derive observationally. Without knowledge of the exact phase of the ISM the star is traversing, estimates of the density could be off by several orders of magnitude (see the discussion on densities, scale heights, and filling factors in Sect. 5). Table 1 lists n_H derived from equating the observed de-projected R₀ with the theoretical R₀, taking n_H as the unknown. The comparison of n_H derived from Eq. (5) with that derived from the observed and predicted R₀ shows that these generally agree within a factor of three (Fig. 11). In a few cases, the derived ISM density is an order of magnitude higher (n_H ≈ 5 to 35 cm^-3) than given by Eqs. (4) and (5). Possibly, for these cases, the star is moving through a diffuse cold medium (CNM) with typical densities n_H ~ 10−100 cm^-3, an order of magnitude higher than those of the WNM (Sect. 5). Alternatively, this apparent discrepancy could also be resolved if the ISM itself has a peculiar flow velocity of v_⋆-v_ISM with respect to the local standard of rest. If the peculiar velocity between the two media is higher, the density required to arrive at the same stand-off distance is lower, since v_⋆ ∝ 1/$\sqrt{{\mathit{n}}_{\mathrm{H}}}$. On the other hand, the derived ISM density also scales linearly with the adopted Ṁ, thus an order of magnitude over-estimate of the mass-loss rate leads to the same order of magnitude over-estimates of n_H. This may explain the high n_H for some objects with high adopted Ṁ, such as W Aql, but not easily for others, such as R Leo, which have already low values for Ṁ. For R Leo, the turbulent features in the bow shock indicate a higher space velocity which is at odds with the observed v_⋆ = 15 km s^-1. If indeed the relative star-ISM space velocity for R Leo is higher, the inferred ISM density will be correspondingly lower (now ${\mathit{n}}_{\mathrm{H}}\mathrm{\propto }{\mathit{v}}_{\mathit{\star }}^{-2}$). For example, for a much higher peculiar velocity of 45 km s^-1 (more in line with the bow shock morphology), n_H = 1.6 cm^-3 which is close to n_H = 1.1 cm^-3 given by Eq. (5). In Sect. 5 we roughly estimated the LSR velocity of the local ISM and found that the corrections to v_⋆ are small ( ≤ 20%) for most objects, with a few exceptions. Only in two cases do these higher relative velocities lead to much smaller values for the ISM density; n_H = 10.1 and 13.5 cm^-3 for μ Cep and U Cam, respectively. The ISM velocity correction does not alter the high densities obtained for e.g. R Leo, W Aql and EP Aqr. Another explaination could be that this ‘denser’ medium is a remnant of an earlier, slower stellar wind. How and if the (tentative) binarity of both W Aql and EP Aqr plays a role is as of yet not clear.

7.6. Presence versus absence of wind-ISM bow shocks

In the previous sections, we have established different morphological classes of bow shocks and examined the basic physics giving rise to a bow shock and discussed the various parameters affecting their size and shape. In this section we explore different properties of the stellar objects as well as the ISM in order to understand which conditions are required to be able to detect bow shocks and/or detached shells around AGB stars and red supergiants.

For 50 out of 78 (63%) AGB stars and red supergiant, we detect detached extended far-infrared emission suggestive of a bow shock, detached ring or irregular extended emission. Limiting the sample to nearby objects (d ≤ 500 pc) eliminates many uncertainties and the “detection-rate” improves to 43 out of 56 objects, or 78%. Restricting the distance further to, for example, d ≤ 300 pc, only marginally improves the detection rate to 28 out of 34 (80%).

To test whether the “detection” and “non-detection” samples have similar or statistically different distributions as a function of d, z or n_H, v_⋆, i, Ṁ, and R₀, we use the two-sample Kolmogorov-Smirnov (K-S) test. The results are given in Table 8. A value close to unity indicates that the two samples (detection and non-detection) have a high probability of being from the same parent distribution. A low probability suggests their distributions differ significantly, thus indicating the parameter is decisive in observing bow shock interaction in our survey. In order to eliminate the effect of large uncertainties in the distances, we perform the K − S test on a distance limited (d ≤ 500 pc) sample. The K − S probabilities of “detections” versus “non-detections” are particularly low (and thus indicative that the two samples have different distributions) for n_H/z, v_⋆, and R₀. The corresponding distribution histograms for the most distinctive parameters, n_H, v_⋆, and R₀, are shown for the different classes (detection/non-detection) in Figs. 13 to 15. This points towards the scenario that the presence of bow shocks is strongly dependent on the stellar velocity (relative to the local medium), the local ISM density, and the resulting standoff distance. The last sets the size of the bow shock region and is apparently determined predominantly by v_⋆ and n_H and not so much by the star’s mass-loss properties (cq. evolutionary phase).

All objects with observed wind-ISM interaction zones (Class I to IV) have R₀  ≤  1 pc, while many class X sources have R₀  ≥  1 pc. The distribution of the observed de-projected R₀ values is very similar to the distribution of predicted R₀ for sources in Class I to IV. However, a quantitative comparison between the predicted and observed stand-off distances does not reveal a strong correlation.

Regardless of the absolute distance (and thus apparent angular stand-off distance) it appears that a necessary, but not sufficient, condition for the detection of bow shocks and rings is that their physical size has to be smaller than about 1 pc. This points towards a physical effect (e.g. reduced surface brightness) instead of an observational limit, although it does not explain the absence of bow shock emission around nine of the observed nearby AGB stars (with known v_⋆ and R₀ < 1pc) included in Class X.

Concluding, we find that bow shocks are detected and predicted to occur for most nearby objects (<500 pc) as well as for objects whose stellar and local ISM properties yield relatively small stand-off distances (<1 pc). For more distant objects the detection is (likely) hampered by lower sensitivity and lower spatial resolution, though we could not find any reasons to suggest that these objects would not have bow shock interaction. Indeed, several distant objects such as AQ And (~800 pc), S Cas (~940 pc), RZ Sgr (~730 pc) reveal evidence for wind-ISM interaction, but these detections might have been fortuitous.

8. Conclusions

This paper presents a morphological inventory of “wind-ISM” bow shocks and “wind-wind” interactions detected in the far-infrared with Herschel/PACS. Five main classes were identified.

• I

  “Fermata”. These objects are characterised by a large arc or shell-like structure spanning an angle of at least 90°. Several objects in this class show the presence of turbulence, i.e. Rayleigh-Taylor or Kelvin-Helmholtz instabilities. Some objects, such as X Pav, EP Aqr, and X Her, reveal the presence of back flowing material in the wake (behind) the star, leading to a distinct bullet shape.

• II

  “Eyes”. These are characterised by two arcs on opposing sides of the central object. These arcs are elliptical and non-concentric. In two cases the arcs connect and form a “lemon” shape morphology.

• III

  “Rings”. This class includes the well known detached shell objects, which are typical of wind-wind interaction and dust pile-up. Younger, smaller rings have a co-spatial counterpart molecular gas ring, while the larger dust shells do not.

• IV

  “Irregular”. This class includes all sources that show extended (non-detached), irregular emission at 70 or 160 μm.

• X

  “Non-detection”. This final class includes all objects for which we do not detect any extended (c.q. resolved) far-infrared emission in the form of a bow shock or ring in the PACS images.

Oxygen rich stars give rise predominantly to “fermata” or “irregular” morphologies. There is tentative evidence that all “ring” stars are thermally pulsing AGB stars. We identify a few cases for which both a detached (inner) shell is found within the (outer) bow shock, suggesting that the detached shell is not due to the wind-ISM interaction, but rather suggest a wind-wind scenario.

The presence or absence of bow shocks is determined by the stand-off distance that depends on relative values of the stellar parameters (v_⋆, v_w and Ṁ) as well as the properties of the ISM (n_H). This is indicated by both the theoretical approximation and the hydrodynamical simulations. The distribution of stand-off distances for the entire survey sample shows a clear separation between objects with and without detected extended emission. Indeed most of the objects (44) assigned to Classes I to IV have R₀ <  1.0 pc. No extended bow shock or detached shell emission has been detected for the 15 objects with predicted R₀ > 1 pc and/or d > 1000 pc. Limiting the sample to d < 500 pc, we derive a R₀ < 1.0 pc, and thus predict the presence of bow shocks, for 8 stars in the “non-detection” class: TW Hor, V Eri, R Dor, R Lep, RS Cnc, HD 100764, RY Dra, and RX Boo. Many of these show extended shells at about 1 MJy sr^-1 sensitivity in Spitzer and AKARI maps (private communication: Ueta for Spitzer and Izumiura for AKARI).

The angular size of the predicted stand-off distances for 6 other objects in Class X are larger than the obtained PACS image maps. In particular, all-sky survey missions, such as AKARI (Murakami et al. 2007) and WISE (Wright et al. 2010) could help in observing extended bow shocks and shells around these objects. Observations of bow shocks and detached shells for distant (>1 kpc) AGB stars and red supergiants require thus higher sensitivities to detect extended emission, but not necessarily higher spatial resolution as most have predicted stand-off distances larger than 1′.

The observed infrared emission indicates the presence of moderate amounts of dust (5 × 10^-8 to 3 × 10^-6   M_⊙). This implies that there is less than 10^-3   M_⊙ of dust and gas present in the bow shock interaction region. The observed dust and gas masses are, in most cases, lower than the that of potentially swept-up ISM material.

This survey represents only a first step in fully characterising and understanding the formation and shaping of bow shocks and detached shells around AGB stars and red supergiants. It is clear that additional simulations trying to represent as closely as possible the different observed morphologies are required in particular to understand the formation of instabilities. Detailed observations of the bow shock spectral energy distribution (dust) and line emission (gas) is needed to quantify the physical conditions and composition of the material in these shocked regions. Potentially, detailed observations of bow shocks around larger samples of AGB stars can be used to independently probe the local ISM density (and possibly the magnetic field) distribution in the Galaxy.

---

Acknowledgments

This research has made extensive use of the Simbad and Vizier, operated at CDS, Strasbourg, France, as well as of NASA’s Astrophysics Data System. We thank Wang Ye at the Department of Physics & Astronomy, University of Kentucky, for providing us with the radiative cooling curve. This work was supported in part by the Belgian Federal Science Policy Office via the PRODEX Programme of ESA (No. C90371). F.K. acknowledges funding by the Austrian Science Fund FWF under project number P23586-N16, R.O. under project number I163-N16. A.J.v.M. acknowledges support from FWO, grant G.0277.08 and K.U. Leuven GOA/09/009. PACS has been developed by a consortium of institutes led by 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). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI (Italy), and CICT/MCT (Spain).

References

All Tables

All Figures