EDP Sciences
Free Access
Issue
A&A
Volume 582, October 2015
Article Number A24
Number of page(s) 20
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201526159
Published online 30 September 2015

© ESO, 2015

1. Introduction

Most stars form in clusters (Lada & Lada 2003) and, in particular, massive stars are predominantly found in clusters and associations (Gies 1987). Furthermore, many O stars not currently in a cluster have been identified as runaway stars (e.g. Gvaramadze et al. 2012). Feedback from massive stars in clusters is important for the enrichment and evolution of the interstellar medium (ISM) of gas-rich galaxies (Krumholz et al. 2014), and it is a crucial ingredient in models of multiple stellar populations in globular clusters (e.g. Gratton et al. 2012). Globular clusters formed early in cosmic time, during the epoch of first galaxy formation, and their multiple populations (e.g. Piotto et al. 2007) demonstrate that somehow gas was able to accumulate, cool, and become gravitationally unstable within the star cluster (D’Ercole et al. 2008; Palouš et al. 2014). It is not at all clear how this happens; a number of theories and their shortcomings are discussed in depth in Bastian et al. (2015). It is fair to say that there is no consensus on how globular clusters formed. Because of this, researchers have begun to look at present-day massive star clusters to try to gain some insight (e.g. Bastian et al. 2013; Bastian & Strader 2014), with the result that little evidence of cold gas and dust has been found in the intracluster medium of nearby massive star clusters.

Models for multiple populations that involve mass loss from massive stars (as opposed to asymptotic giant branch stars, e.g. D’Ercole et al. 2008) all concentrate on winds from hot stars (e.g. Cantó et al. 2000; Wünsch et al. 2008; de Mink et al. 2009; Krause et al. 2013). The difficulty these models face is how to cool down the hot shocked wind to temperatures low enough to allow star formation (Palouš et al. 2014). Surprisingly, these models do not explicitly consider red supergiants, even though Galactic stars with initial mass M ≲ 30 M lose more mass as red supergiants than they do as main sequence stars (Chiosi & Maeder 1986). Their winds are already cold (temperature 100500 K; e.g. Schuster et al. 2009; Le Bertre et al. 2012), neutral, dusty, and dense, and so if this wind material can remain in a star cluster then it is the obvious candidate to host second generation star formation. Furthermore, evolutionary models of very metal-poor stars with initial masses above 100 M show that they may evolve to red supergiants already on the main sequence (Marigo et al. 2003; Szecsi et al. 2015; Köhler et al. 2015; Yoon et al. 2012). For this reason we consider it important to study mass loss from red supergiants in star clusters, to constrain whether their wind material can remain cold and dense in such a harsh environment.

A growing number of red supergiants are found to be surrounded by nebulae with emission lines from ionized gas, showing that their winds are photoionized by external radiation fields. This was first proposed for NML Cyg (Morris & Jura 1983), then for IRS 7 in the Galactic centre (Yusef-Zadeh & Morris 1991), and more recently W26 (Wright et al. 2014) and IRC 10414 (Gvaramadze et al. 2014; Meyer et al. 2014a). WOH G64 in the Large Magellanic Cloud also shows similar nebular emission lines (Levesque et al. 2009). Mackey et al. (2014) proposed that Betelgeuse’s wind is also externally photoionized, to explain the presence of a static neutral shell around the star (Le Bertre et al. 2012). We expect on this basis that many red supergiants have photoionized winds. Mackey et al. (2014) showed that the winds of red supergiants could be decelerated and accumulated into a dense shell, if the wind is exposed to an external ionizing radiation field from all sides. This could have important consequences for the fate of the mass lost during the red supergiant phase, and also for any subsequent evolutionary phases and the observational characteristics of the star’s eventual supernova explosion (Mackey et al. 2014). If this mechanism is effective, then the mass lost by red supergiants could remain within the parent star cluster for much longer than if the wind is freely expanding, increasing the likelihood that secondary star formation could happen.

Westerlund 1 (Westerlund 1961, 1987) is probably the most massive young star cluster in our Galaxy (Brandner et al. 2008), and its age (4 ± 1 Myr) is such that it contains a population of extreme supergiant stars (Clark et al. 2005). Only at this age do star clusters contain the hottest (Wolf-Rayet stars) and coldest (red supergiants) stages of massive star evolution simultaneously. The red supergiant W26 (RA: 16 47 05.403; Dec: −45 50 36.76; SIMBAD identifier: Cl* Westerlund 1 W 26) is surrounded by a bright radio nebula (Clark et al. 1998; Dougherty et al. 2010), also seen in mid infrared (Clark et al. 1998) and in narrow-band optical emission (Wright et al. 2014). The wind of W26 is so bright in typical tracers of photoionized gas (Hα and [N II] forbidden-line emission) that a number of authors have suggested that it could be externally photoionized by radiation from the O-type and Wolf-Rayet (WR) stars in the cluster (Dougherty et al. 2010; Wright et al. 2014; Mackey et al. 2014). The main alternative explanation for the ionized nebula around W26 is that it is hydrodynamically confined by the pressure of the intracluster gas. This could be ram pressure from stellar winds or a global cluster wind (Dougherty et al. 2010), or thermal pressure from thermalised shocked wind bubble and recent supernovae that must have occurred in Westerlund 1 (Muno et al. 2006; Kavanagh et al. 2011).

External photoionization and hydrodynamic confinement must both be present, because the combined winds of the evolved stars will create a high pressure environment (Muno et al. 2006) and there are many hot stars in the cluster that emit ionizing photons, e.g., 24 WR stars (Crowther et al. 2006). The challenge then is rather to understand which process is producing the ionized nebula in W26, and by consequence in the winds of red supergiants in other star clusters. The radius at which a stellar wind becomes ionized is well defined for a specified external radiation field (Morris & Jura 1983), and the hydrodynamic confinement radius is set by the balance between the ram pressure of the stellar wind and the thermal/ram pressure of the environment. Unfortunately there are no good estimates of the mass-loss rate or wind velocity of W26, and we do not know the physical distance separating it from the other stars (such as the supergiant sgB[e] star, W9) that are nearby in projection.

To make progress, we take the existing spectroscopic data on W26 and its circumstellar environment, and compare this to the predictions of numerical simulations of an externally photoionized wind, following the method of Mackey et al. (2014). The main observables are the radial velocity of the emitting gas, the intensity of this emission, and the line ratio, as a function of position along the slit in the Hα and [N II] lines (in this paper we focus entirely on the brighter [N II] line at 6584 Å). As shown below, these data already provide contraints that are difficult to satisfy with a spherically symmetric photoionized wind.

The observational data, reduction methods and results are described in Sect. 2. The results of numerical simulations are described and compared with observations in Sect. 3. We discuss our results in Sect. 4 and conclude in Sect. 5.

2. Dataset and analysis processing

W26 was observed with the FOcal Reducer and low dispersion Spectrograph (FORS 2; Appenzeller & Rupprecht 1992; Appenzeller et al. 1998) attached to the Cassegrain focus of the ESO Very Large Telescope of the Paranal Observatory, Chile. The observations were carried out on the 16th of April 2011 under the ESO program ID 087.D-0673(A). The data were obtained in long-slit mode using the 1200R and 1028Z grisms and, a slit size of 0.3 arcsec to obtain the highest resolution available. This instrumental configuration offers a wavelength range between ~5700−9300 Å and a resolving power of R ~ 7000 (Negueruela et al. 2010; Clark et al. 2013). The slit location and orientation is shown in Fig. 2 of Wright et al. (2014), where the data were first presented.

The data reduction was performed employing IRAF1 (Tody 1993) and IDL routines following the standard steps and tasks for long-slit spectral data reduction (consisting of bias subtraction, flat-fielding, spectral extraction with and without background subtraction, and wavelength calibration) and analysis.

2.1. W26 stellar parameters

We used the background-subtracted, wavelength-calibrated FORS2 spectra of W26 to roughly estimate the effective temperature (Teff). We did this by comparing the observed spectrum with synthetic fluxes calculated with marcs stellar atmosphere models2 (Gustafsson et al. 2008) with spherical symmetry and solar metallicity (Grevesse & Sauval 1998) spanning between 3300 K and 4500 K in Teff, in steps of 100 K, and 0.5 and +1.0 in log g, in steps of 0.5 dex.

Given the very low resolution of the FORS2 spectrum, hence the heavy line blending, there was no continuum available that could be used to properly normalise the observed spectrum. For this reason, we normalised the observed spectrum using each model flux at a time as a reference before determining the total deviation between observed and model spectrum. We obtained a best fitting effective temperature of Teff = 3600 K, which corresponds to a spectral type of about M2.5 (van Belle et al. 1999), in agreement with previous spectral type determinations (Clark et al. 2010); a comparison of the model with the FORS spectrum is shown in Fig. 1. The formal uncertainty in Teff is 200 K, but systematic uncertainties in the modelling are likely to be larger. It is important to recall here that W26 is a known variable star with a spectral type varying between M2 and M5 (Clark et al. 2010).

2.2. Radial velocity and intensity distribution maps

The kinematic study presented here is based on the radial velocity distributions of Hα6562.79 Å and [N II] 6583.45 Å along the spatial direction of the FORS 2 long-slit spectra. This analysis was done with the 1200R grism data since these two transitions can only be reached in its wavelength coverage. To explore the entire spatial information available in the spectra, each row of the science images was sliced and individually calibrated by the equivalent row in the calibration arc-lamp spectrum. The targeted lines were fitted by a Gaussian function per row, identifying the position of the emission peaks and their amplitudes with respect to the spectral continuum. The velocities are presented in the heliocentric reference system. Previous studies reported the possible presence of small instrumental flexures that might lead to shifts in the wavelength calibration as large as 5 km s-1 (Bagnulo et al. 2013). We corrected for such a possible instrumental bias by further calibrating the spectrum using two telluric sky emission lines at 5889.951 and 6300.304 Å.

thumbnail Fig. 1

Comparison between observed FORS spectra for W26 (black solid line) and MARCS model fluxes (dashed red line) for a Teff = 3600 K supergiant. The wavelength gap where no spectrum was obtained arises from a gap between the two grisms.

Open with DEXTER

thumbnail Fig. 2

Spectral observations of W26 and neighbouring stars. a) Radial velocity of the Hα (black dots) and [N II] (blue triangles) lines along the slit in the W26 co-moving reference system (~−32 km s-1), with position measured relative to W26. The directions along the slit are also indicated: WNW for Westnorthwest and ESE for Eastsoutheast. b) Line brightness (intensity) in the two spectral lines along the slit (arbitrary units). c) Fitted Gaussian linewidth, σ in the two spectral lines along the slit. The three stars included in the slit, and no further away than 0.4 arcmin from W26, are labelled according to Negueruela et al. (2010).

Open with DEXTER

Figure 2 displays the radial velocity (vlos, panel a) and line brightness (b) spatial distributions of Hα and [N II] along the FORS 2 slit up to a maximum distance of 0.4 arcmin from W26. The radial velocity was shifted to the W26 Heliocentric radial velocity co-moving frame, ≈−32 km s-1 (Negueruela, priv. comm.), obtained using data from the AAOmega spectrograph on the Anglo-Australian Telescope, and reduced using the techniques described in Marco et al. (2014), with a precision of 24 km s-1. We derived uncertainties of 4 and 7 km s-1 for Hα and [N II] radial velocities, respectively.

We also measured the linewidths of the Hα and [N II] lines along the slit, shown in panel (c) of Fig. 2, and from the Gaussian fits we obtained σ ≈ 30 km s-1 for both lines. For Hα most of the measurements are in the range 3040 km s-1, whereas for [N II] they are mostly in the range 2040 km s-1. This is consistent with the instrumental resolution of FORS 2. The only exception is near W9, where the very bright Hα line has a width σ ≈ 55 km s-1. There is an apparent trend that σ decreases from left to right across the nebula around W26, or with increasing distance from W9. It is unclear if this is significant because we are so close to the instrumental resolution; higher resolution spectra could confirm this. Such a trend might be expected if the wind from W9 were directly colliding with the wind of W26 (see Sect. 4).

For photoionized gas at T ≈ 104 K, Hα has a thermal linewidth σ ≈ 10 km s-1, and the [N II] line is significantly narrower. We therefore cannot hope to resolve any substructure in the nebular emission, for example blue- and redshifted components of the wind of W26. This means that the only information we can reliably obtain from the spectral lines are the radial velocity of the peak and the line brightness. It should be noted, however, that the observed line could be a blend of multiple unresolved lines with differing radial velocities. Furthermore, the line brightness is uncalibrated which makes absolute flux/luminosity measurements impossible, and so we are limited to comparing the relative line brightness at different locations along the slit. Spectroscopy with significantly higher spectral resolution, able to resolve multiple line components in photoionized gas, would provide much stronger constraints on models.

W9 has an enormous Hα luminosity (Westerlund 1987) with a measured equivalent width W = −520 ± 30 Å (Clark et al. 2005), which they point out is unprecedented for an early type, emission line star (or W = −640 ± 40 Å; Clark et al. 2013). The Hα line is much brighter than the [N II] emission lines. The diffuse emission in the ISM far away from W9 and W26 also has stronger Hα than [N II] emission. In contrast, the nebula around W26 has [N II] emission significantly brighter than Hα, as noted by Wright et al. (2014). This is not generally obtained in shock models with solar abundances, which predict intensity ratios [N II]/Hα< 1 (e.g. Raymond 1979), and so can be a signature of enhanced nitrogen abundance. Stellar evolution models predict that the nitrogen abundance in red supergiant winds should be enhanced by a factor of 35 (Brott et al. 2011), and this has been measured in circumstellar nebulae around a number of evolved massive stars including the runaway red supergiant IRC 10414 (Gvaramadze et al. 2014; Meyer et al. 2014a, and references therein). Alternatively, ionization fronts can show enhanced [N II] emission under certain conditions ([N II]/Hα ~ 1−3), even for solar abundances (Henney et al. 2005).

The three data points closest to W26 have significantly larger uncertainty than the surrounding points, because the continuum flux from W26 itself fills in the line at radial velocities near zero. This biases the radial velocity measurement of the line to the blue, and so these three points blueshifted by 3050 km s-1 are probably spurious (the line intensities in these points are similarly more uncertain than neigbouring points). Aside from these points, the radial velocity of the nebula around W26 is almost uniform at ≈−20 km s-1.

In the low-density limit, assuming that all N is in the form of N+ and that H is fully ionized, the emissivity ratio of [N II] to Hα can be expressed as (Dopita 1973) (1)where we assume T ≫ 350 K, and where f(N) is the number fraction of N (relative to H) for which 6.8 × 10-5 is the solar photosphere abundance (Asplund et al. 2009). In our models we assume f(N) = 2 × 10-4, appropriate for nitrogen-enriched winds of rotating red supergiants (Brott et al. 2011), leading to a pre-factor of 0.341 instead of 0.116. For typical H ii region temperatures (700010 000 K) the ratio is between 0.5 and 1.5 for solar abundance, but always >1 for the enhanced N abundance.

The low density limit applies (neglecting collisional de-excitation) because in our simulations the electron number density, ne, is always less than the critical density for the 6584 Å line of [N II], which is ne,crit = 8.6 × 104 cm-3 (Osterbrock 1989, Table 3.11). For reference Dougherty et al. (2010) estimate ne ≈ 5 × 103 cm-3 from radio observations, if one assumes a constant density in the emitting region. Comparing this estimate to the critical density, we expect a reduction in [N II] emission of <10% because of collisional de-excitation. Taking into account the possibility of doubly ionized N and of collisional de-excitation, any constraints we obtain on the nitrogen abundance are lower limits.

3. Simulations of an externally photoionized wind

Mackey et al. (2014) showed that external photoionization of the winds of red supergiants can have strong hydrodynamical effects, in particular decelerating the wind into a dense and massive shell that could contain up to 35% of the mass lost. This completely changes the density and velocity structure of the wind, and hence its emission properties. The wind structure was calculated by Mackey et al. (2014) using analytic and numerical methods, in both cases assuming an isothermal equation of state (with temperatures of T = 102 and 104 K in the neutral and ionized parts of the wind, respectively). This allowed them to predict 21 cm emission from the neutral shell around Betelgeuse, but this is a special case because it is so close to Earth that the wind can be spatially resolved at 21 cm.

Winds around other red supergiants have been detected (and spatially resolved) in optical nebular lines (Wright et al. 2014; Gvaramadze et al. 2014). We wish to model these lines but they are very senitive to temperature (see Eq. (1)), so it is important to calculate the gas temperature self-consistently with heating and cooling functions, particularly near the ionization front where most of the emission is concentrated (e.g. Henney et al. 2005). The method we use to do this is described below.

Our models provide the first predictions for nebular emission from ionization fronts in red supergiant winds, including both radiation-hydrodynamics and postprocessing with radiative transfer. The wind properties of red supergiants are often poorly constrained from both theory and observations (see e.g. the review by Mauron & Josselin 2011). Optical spectroscopy of nebular lines offers a new avenue to measure the wind velocity and density (and hence mass-loss rate, ), through line intensities and radial velocities. This is of crucial importance for constraining the evolution of evolved stars.

3.1. Input physics for simulations

We use the pion code (Mackey 2012) in spherical symmetry with a wind expanding from the origin (r = 0) and an ionizing radiation (extreme ultraviolet, EUV, with > 13.6 eV) source at radius r = ∞. The method largely follows previous calculations of photoionization-confined shells around red supergiants (Mackey et al. 2014), except that we use non-equilibrium heating and cooling of the gas (Henney et al. 2009), non-equilibrium ionization of hydrogen, and multi-frequency EUV radiation. This calculates gas temperature more realistically than in Mackey et al. (2014). The implementation is almost exactly as described in Mackey et al. (2015), where the scheme was used to describe expanding H ii regions around O stars; we refer the reader to Sect. 2 of this work for details. The rate equation of hydrogen ionization is solved including photoionization, radiative recombination, and collisional (electron impact) ionization. We make the approximation that helium is singly ionized whenever hydrogen is, so the electron fraction is 1.1 times the H+ fraction (for helium number fraction nHe/nH = 0.1). We consider a blackbody radiation spectrum with radiation temperature Tr = 4 × 104 K, in the energy range 13.654.4 eV, and also far-ultraviolet (FUV) heating radiation following Henney et al. (2009) with equal photon flux in FUV and EUV photons. This algorithm for thermal and ionization evolution captures the strong cooling in the neutral gas to tens or hundreds of Kelvin and the photoheating to T ≈ 7000−10 000 K in photoionized gas, including spectral hardening (and extra heating) near the ionization front (e.g. Henney et al. 2005), which is important because [N II] and Hα have different temperature scaling (Eq. (1)).

3.2. Wind properties for W26

There is considerable uncertainty in the wind properties of W26 (e.g. Wright et al. 2014), so we have run simulations for two different mass-loss rates, = 2 × 10-5 and 10-4 M yr-1, and four different wind speeds, v = 15, 20, 25, and 30 km s-1. According to the Mauron & Josselin (2011) reformulation of the mass-loss rate of Nieuwenhuijzen & de Jager (1990), we can estimate (2)Using Teff = 3600 K (estimated above) and L = 3.5 × 105L (Wright et al. 2014) we obtain = 4.2 × 10-5 M yr-1. This is only a very approximate estimate, however, as individual red supergiants often have mass-loss rates a factor of 3 larger or smaller than empirical relations (Mauron & Josselin 2011). Furthermore, Mauron & Josselin (2011) show that measured mass-loss rates for the same object often differ by a factor of a few depending on the technique used.

Fok et al. (2012) compare near- and mid-infrared photometry for W26 with models and estimate = 2.2 × 10-5 M yr-1, although their model requires W26 to have a total luminosity (106.03L) that is much larger than the above-quoted estimates. Dougherty et al. (2010) study the nebula around W26 with radio observations, but did not obtain a measurement of the mass-loss rate. They estimate that the nebular mass of W26 is comparable to that surrounding VY CMa (Smith et al. 2001). VY CMa has a time-averaged ≈ 4 × 10-5 M yr-1 (Decin et al. 2006; Mauron & Josselin 2011) from gas emission measurements, but dust emission gives much larger estimates: = (1.6−4) × 10-4 M yr-1 (Danchi et al. 1994; Smith et al. 2009; Mauron & Josselin 2011). The discrepency may arise partly because the mass loss is variable. On the other hand, Dougherty et al. (2010) estimate = 2 × 10-5 M yr-1 for W237, another (presumably coeval) red supergiant in Westerlund 1, and then note that it has a similar nebular mass to W26. These results show that the two mass-loss rates that we use are reasonable estimates for W26, although they are not hard lower or upper limits.

The distance from W26 of the edge of the Hα emission is about 0.02–0.04 arcmin on the Westnorthwest side (hereafter WNW) and 0.04–0.06 arcmin on the Eastsoutheast side (hereafter ESE); see Fig. 2. We adjust the ionizing flux to set up simulations that have a photoionized nebula at a separation from W26 that is appropriate for each side. Westerlund 1 contains a number of Wolf-Rayet and O stars with ionizing photon luminosities Q0 ~ 1049 s-1, providing ionizing fluxes at W26 on the order of Fγ ~ (1010−1013) cm-2 s-1 (see Sect. 4.6), so the fluxes that we use are about what we expect for the environment of W26.

3.3. Simulation parameters

We have 16 simulations for two mass-loss rates, two positions for the ionization front, and four wind velocities. The identifier (ID) of each simulation and the parameters used are listed in Table 1. In the main text we present only the simulations with = 2 × 10-5 M yr-1 and ionizing fluxes chosen for the ESE side of the nebula (i.e. only v is varied). We focus on the ESE side of the nebula because it extends further from the star and hence the data are less contaminated by the central star. The results from the simulations with = 10-4 M yr-1 are very similar to those with = 2 × 10-5 M yr-1 except that the wind in the former case is five times denser than the latter, and so the line emission is approximately 25 times brighter. It is not possible to calibrate the data in absolute flux, and so we do not know the absolute line brightness and cannot constrain the wind density (but see Sect. 4). For completeness, the other simulations are presented in Appendices A and B.

All simulations use 10 240 uniformly spaced grid zones from r = 9.5 × 1015 cm (0.0031 pc) to 6.172 × 1017 cm (0.2 pc). This is more than sufficient to resolve the structure of the flow and the internal structure of the shocked shell. The simulations were run for 0.1 Myr, after which the star has lost 10 M for = 10-4 M yr-1 or 2 M for = 2 × 10-5 M yr-1. For v ≥ 20 km s-1 the simulations had reached a steady state at this time. The shocked shell has not reached a steady state for simulations with v = 15 km s-1 (it is still accumulating mass), but this does not affect the ionized gas line emission. We checked this by running some simulations for 0.2 Myr and comparing the resulting emission properties.

Table 1

Parameters used for the 16 spherically symmetric radiation-hydrodynamics simulations.

3.4. Structure of the photoionized wind

The structure of the photoionized wind as a function of radius is shown in Fig. 3 for simulation M5V15E. The ionization front is D-type (dense; Kahn 1954) because the wind velocity satisfies v ≤ 2ai, where ai is the isothermal sound speed in photoionized gas. It has the classical structures of (from right to left) the ionization front, a dense shocked shell, and a shock in the neutral gas. Mackey et al. (2014) assumed isothermal neutral gas, and so the shock was very compressive, whereas here the shock is adiabatic with a factor of 4 increase in density. Then there is a post-shock cooling region where the temperature decreases (from r ≈ 0.043−0.0455 pc), and density increases to compensate (velocity also decreases to conserve mass). Finally there is a dense cooled shell bounded by the cooling layer and the ionization front.

thumbnail Fig. 3

Radial structure of the wind in simulation M5V15E: we show the H number density nH, H+ fraction, gas radial velocity vr, and temperature T on a logarithmic scale as a function of distance from the red supergiant.

Open with DEXTER

Simulations with v = 20 and 25 km s-1 are shown in Figs. 4 and 5, respectively. For v = 25 km s-1 the wind is too fast to sustain a shocked region and the ionization front is R-type (rarefied), characterised by weak density and velocity changes across the front (Kahn 1954). The temperature increases smoothly in the wind from a few hundred to about 10 000 K. The ionization front is also much broader because there is no dense shell to absorb all of the photons. Simulation M5V20E (Fig. 4) is a transition case on the boundary between D-type and R-type ionization fronts. There is a shock in the wind, but the shock cannot propagate far upstream from the ionization front and so there is too little time for the post-shock gas to cool. The separation between the shock and the ionization front is calculated in Mackey et al. (2014) and goes to zero as v → 2ai. In this case there is a shocked shell, but it cannot accumulate significant mass.

thumbnail Fig. 4

As Fig. 3 but for simulation M5V20E.

Open with DEXTER

thumbnail Fig. 5

As Fig. 3 but for simulation M5V25E.

Open with DEXTER

Figures 35 also plot , which determines (along with temperature) the emissivity of the nebular emission lines. The peak value of is always right at the ionization front, and drops by an order of magnitude between r = 0.044 pc and r = 0.065 pc. We therefore expect most of the line emission to arise from a relatively thin shell of ionized gas around the ionization front.

3.5. Properties of the Hα and [N II] line emission

We project the spherically symmetric wind onto the plane of the sky and trace rays through the wind for various impact parameters, b, the radius of closest approach for each ray to the red supergiant. The emission and absorption coefficient for the Hα and [N II] lines are calculated for every grid zone and tabulated. Then each ray is split up into line segments corresponding to the grid zones, and the radiative transfer equation is solved analytically along each segment assuming a constant source function (see e.g. Rybicki & Lightman 1979), first inwards to r = b and then back outwards. To calculate velocity information we obtain the line-of-sight radial velocity of the gas from the gas dynamics, and the emission is then thermally broadened about this value using the gas temperature. We use 501 velocity bins between −50.1 and 50.1 km s-1 to ensure that the lines are always well resolved in velocity.

The radial velocity of peak emission is then extracted for each pixel, and the line brightness is obtained by summing the emission at all velocities for each pixel. The line ratio is the ratio of these line intensities. We convert from impact parameter, b, to an angular distance from W26, θ, by assuming a distance to Westerlund 1 of 3.55 kpc (Brandner et al. 2008), consistent with the distance derived more recently from eclipsing binary stars (3.7 ± 0.6 kpc) by Koumpia & Bonanos (2012). We also spatially smooth the intensity maps to a resolution of 1 arcsec at the distance of Westerlund 1, for easier comparison with W26.

The results are plotted in Fig. 6 for simulation M5V15E, with the observational data overplotted, as a function of distance from W26 measured in arcminutes. Line intensities for [N II] and Hα are plotted in panel a, the line ratio in panel b, and radial velocity of the line peaks in panel c. The relative normalisation of the simulation to the data is chosen by setting the mean [N II] line intensity to unity between the two radii shown by vertical dotted lines in panel a. Note again that the simulation results assume nitrogen is enhanced in the stellar wind by a factor of 2.8 (see Sect. 2). The spectral data close to W26, within about 0.015 arcmin, are less reliable than at larger separations because of stellar continuum subtraction (see Sect. 2).

This simulation has a D-type ionization front with a dense static shell, so the ionization front is much denser than in the higher velocity simulations. This leads to a thinner region of peak emission because the newly ionized gas accelerates, and so its density (and hence emissivity) decreases more rapidly with radius than for a constant velocity wind. This simulation therefore has the most limb-brightened emission, compared with simulations with faster winds.

thumbnail Fig. 6

Hα and [N II] spectral lines from the simulation M5V15E, with = 2 × 10-5 M yr-1, v = 15 km s-1, and Fγ = 1.70 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The dashed lines show the simulation, solid lines smoothed to 1 arcsec resolution at the distance of Westerlund 1, and the dotted lines with points show the observations to the ESE of W26. The panels show a) line brightness; b) line ratio [N II]/Hα; and c) radial velocity of the peak emission, all as a function of distance from the star. Vertical dotted lines in panel a) show the range of radii used to set the normalisation of the simulated and observed [N II] line brightness.

Open with DEXTER

The peak emission in panel a is at the same position as the observations by construction, because we constrained the ionizing flux to locate the ionization front at this radius. The observed emission is much less peaked than the synthetic data for both lines, although this can partly be explained by the finite spatial resolution of the observations. The degree of limb brightening constrains the thickness of the shell of brightly emitting gas, with a thinner shell being more limb-brightened. The data thus appear to imply a thicker shell of brightly-emitting gas compared with what we find from this simulation.

The [N II] brightness decreases more rapidly with distance from W26 than the synthetic data, but this is less clear with the Hα brightness because the ratio between maximum and minimum observed brightness is not as large as for [N II]. The observed diffuse Hα emission at large radii is a constant component from Westerlund 1 and not related to the W26 nebula (see Fig. 2), even though its brightness is similar to our predictions at separation 0.070.1 arcmin. Subtraction of this constant component complicates any estimate of the rate of decrease at these larger radii, but it may be that the observed Hα also decreases more rapidly with radius than we predict. The rate that the brightness decreases with distance from the peak constrains the density of the flow as a function of radius if the gas is isothermal (assuming spherical expansion). The data then indicate that either nitrogen begins to be doubly ionized further from the star, or the emitting gas density decreases much more rapidly than our models predict. This could arise if the wind of W26 is hydrodynamically confined by much hotter intracluster gas (see Sect. 4).

The brightness ratio, [N II]/Hα, is plotted in panel b and shows that within the bright part of the nebula the synthetic and observed data have a broadly constant line ratio that is 2 for the observations and 2.5 for the simulations. If the wind were not enriched in nitrogen, the simulations would predict a line ratio of 0.9 in this region, inconsistent with the results. The observations therefore indicate that nitrogen is enhanced by a factor of ~23 in the nebula. If we consider the nitrogen abundance to be a free parameter, we obtain good agreement between simulations and observations for f(N) = 1.60 × 10-4, or 2.35 times the solar abundance (Asplund et al. 2009). This enhancement is strong evidence that the nebula contains gas processed by the CNO cycle, and hence that it consists of mass lost by W26, at least out to a distance of 0.06 arcmin.

The line ratio decreases in the simulations from the ionization front (θ ≈ 0.045 arcmin) outwards, whereas the observed line ratio stays constant to the edge of the bright nebula and then decreases sharply from 1.5 to 0.5 between θ = 0.06 and θ = 0.08 arcmin. The synthetic data never decrease to such a small line ratio, so again we suggest for this low-emission gas that either nitrogen is more highly ionized or the gas is not enriched in nitrogen (i.e. it is not part of the wind of W26). Whatever the explanation, clearly there is something happening in the region between θ = 0.06 and θ = 0.08 arcmin that is not seen in the synthetic data.

This is further shown by looking at the radial velocity of the line peak as a function of position, shown in panel c. The inner part of the nebula has radial velocity vrad ≈ −15 km s-1, similar to the simulated lines. The observations then show that the brightest emission at 0.0350.055 arcmin is also the most blueshifted, whereas the simulations show that the radial velocity of the brightest emission is zero. This is the biggest discrepency between our model and the data: both the line ratio and the radial velocity stay approximately constant throughout the bright nebula, while the simulations show a change at the brightest part of the nebula. The decrease in observed line ratio is spatially correlated with the decrease in radial velocity.

Figure 7 shows the same plots for simulation M5V25E, with a faster wind such that the ionization front is R-type, with only weak density and velocity changes in the wind (see Fig. 5). The agreement of the line brightness with observations is better than for M5V15E, in that the line emission is not so strongly peaked at the ionization front, but the decrease in intensity at θ> 0.05 arcmin still has the wrong shape for [N II]. The line ratio also agrees well for θ ≤ 0.06 arcmin. We see the same discrepencies for this simulation as for M5V15E, however. The simulation predicts that the line goes to zero radial velocity at the emission peak, whereas the data show the maximum radial velocity at the peak. The same holds true for the line ratio, which in the observational data remains large as long as the radial velocity is significantly different from zero.

thumbnail Fig. 7

Hα and [N II] spectral lines from the simulation M5V25E, again after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. 6.

Open with DEXTER

The situation is somehat different on the opposite side of W26, Westnorthwest (WNW) in the direction of the extreme sgB[e] star W9. Here the nebula is confined closer to W26, possibly because of the proximity of W9. The observed data shows no limb-brightened peak in the nebular emission, rather just a shoulder of emission between 0.02 and 0.03 arcmin. The results of a comparison between simulation and observation are very similar to what we have already discussed for the ESE side of the nebula, and are shown in Appendix B.

3.5.1. Wind velocity as traced by Hα and [N II] lines

Figure 8 plots the blueshift of the two lines as a function of wind velocity for the four sets of simulations (two values of and two sets of ionizing fluxes). For slow winds with D-type ionization fronts (v ≲ 22 km s-1) the ionized gas is accelerating outwards from the outer edge of a static shell, and so the line velocity is independent of wind velocity. In this velocity range the Hα line is blueshifted by about 2 km s-1 more than the [N II] line for all simulations. This is a temperature effect: Eq. (1) shows that the two lines will be produced from different depths in the ionization front (see e.g. Hester et al. 1996; Henney et al. 2005). We also find that the simulations with denser winds ( = 10-4 M yr-1) are blueshifted by about 2 km s-1 compared with simulations that have less dense winds. This is a density effect because within the ionization front, the heating rate scales linearly with density, whereas the cooling rate is quadratic for collisional coolants, implying that the temperature structure within the front should be density-dependent. The ionization front thickness is inversely proportional to density, and so its internal structure is also expected to depend on . We find that for a given temperature in the photoevaporation flow, the denser winds have an expansion velocity about 2 km s-1 larger than the less dense winds.

This result for D-type ionization fronts likely depends weakly on the temperature of the ionizing radiation field: harder radiation fields create hotter ionization fronts, faster photoevaporation flows, and hence the lines will probably be more blueshifted. The metallicity of the wind is also important, because lower metallicity gas has a higher equilibrium temperature, which will also drive faster photoevaporation flows.

For faster winds (v ≳ 22 km s-1) the ionization front becomes R-type, and this has very little effect on the wind dynamics except for a slight deceleration. The ionized gas exiting the ionization front has a velocity given by (Mihalas & Mihalas 1984, chap. 106) (3)where the term (ai/v)2< 0.25 for R-type ionization fronts. This explains why the line velocity is always less blueshifted than v, and also why both lines approach v for larger v. Again, [N II] is less blueshifted than Hα because it forms deeper in the ionization front before thermal re-acceleration of the wind can take place. We conclude that these lines are useful but imperfect tracers of the wind velocity for red supergiants, with their accuracy increasing as wind velocity increases.

thumbnail Fig. 8

Comparison between wind velocity, v, and the blueshift of the Hα and [N II] spectral lines for lines of sight close to the star. The boundary between D-type and R-type ionization fronts is at v ≈ 22 km s-1. The red lines are for Hα and the blue for [N II].

Open with DEXTER

3.5.2. Limb brightening of the Hα and [N II] lines

thumbnail Fig. 9

Comparison between wind velocity, v, and the degree of limb brightening of the Hα and [N II] spectral lines, defined as the ratio of the peak brightness to that at lines of sight near the star (r = 0). The boundary between D-type and R-type ionization fronts is at v ≈ 22 km s-1. The red lines are for Hα and the blue for [N II]. Panel a) shows the raw simulation data, and panel b) when smoothed to a spatial resolution of 1 arcsec. The lines in panel b) follow the legend in panel a).

Open with DEXTER

Figures 6 and 7 show that the slow wind simulations have more limb-brightened emission than faster winds. In Fig. 9a this is quantified with plots of the ratio, rLB, of the peak to the central intensity for all 16 simulations and for both spectral lines. This confirms the trend that the two slow-wind simulations (v = 15 and 20 km s-1) with D-type ionization fronts have stronger limb brightening than the two fast-wind simulations (v = 25 and 30 km s-1) with R-type ionization fronts. The slow-wind simulations have rLB ≈ 3.5−5.5 for Hα and rLB ≈ 4.5−6.5 for [N II]. The fast-wind simulations have rLB ≈ 3.3−4.3 for Hα and rLB ≈ 4−5.1 for [N II]. This arises, as noted above, because the emitting region is thinner for the slow-wind simulations, and limb brightening is stronger the thinner the emitting shell. [N II] emission arises from a thinner shell of emitting gas within the ionization front than Hα emission, and so it has larger rLB than Hα for the same reason.

The simulations with denser winds have smaller rLB, which arises because of self-absorption. Lines of sight towards the star have very little absorption (from the front side; the back side is invisible) because the wind density decreases outwards. Lines of sight at the limb pass through much more material, however, and the opacity is significant for the simulations with the densest winds.

Figure 9b shows, however, that the limb brightening is much decreased when the data are smoothed to a resolution of 1 arcsec, with most values of rLB lying between 2 and 3.5. We still see the trend that denser winds have less limb-brightening, but the difference between Hα and [N II] is almost completely gone. Also, the trend of decreasing rLB with increasing v is no longer present, and is even reversed for simulations with dense winds. The smoothed values are still somewhat larger than the observed ratio for W26, rLB ≈ 1.5.

4. Discussion

4.1. Constraining wind properties using photoionized winds

Photoionized red supergiant winds can be extremely bright, being very strongly irradiated by nearby stars. The peak brightness of the simulated nebulae in Hα (smoothed to a spatial resolution of 1 arcsec for the distance of Westerlund 1) ranges from, Imax ≈ 10-13 erg cm-2 s-1 arcsec-2 for M5V30E to Imax ≈ 10-11 erg cm-2 s-1 arcsec-2 for M4V20W. This is, for example, 200−20 000 times brighter than the photoionized bow shock around the red supergiant IRC 10414 (Meyer et al. 2014a), and 1–100 times brighter than the Hα emission from the brightest of the Eagle Nebula pillars (Hester et al. 1996). The overall brightness scales approximately with the square of the electron density (hence the wind density). This can potentially provide a good constraint on the mass-loss rate of red supergiants, in the same way that radio data gives the emission measure and hence electron density (with assumptions about geometry).

Using the brightness estimate from Wright et al. (2014) for the circumstellar nebula around W26, and assuming the same extinction as they find for the nearby triangular nebula (8.4 mag at Hα), the mean intrinsic brightness of the circumstellar nebula is I(Hα) ~ 6 × 10-12 erg cm-2 s-1 arcsec-2. This puts the nebula at the upper end of the emission range of our models, and may be evidence that the mass-loss rate is closer to our upper value than our lower value. Our simulations have similar electron densities (few × 103 cm-3, see Figs. 35) to the radio estimates from Dougherty et al. (2010), so this agreement is not surprising. This holds even though our model does not match the data very well, because the mass-loss rate is the primary parameter determining the gas density at a given distance from the star.

Winds that are slow enough to allow D-type ionization fronts (≲ 22 km s-1 in these models) have blueshifted lines for which the velocity is determined by the ionization front and not by the wind. We find line blueshifts of 16−19 km s-1, although this depends somewhat on wind density (and probably also on the radiation temperature). For faster winds an R-type ionization forms, and this has only weak dynamical effects on the wind. In this regime, nebular lines become good tracers of the wind velocity, although the line velocity is always an underestimate. This is potentially very useful for estimating wind terminal velocities because the best estimates (usually from masers) come from the wind acceleration region itself and require models for the turbulent velocities in the wind (e.g. for NML Cyg, Zubko et al. 2004; Nagayama et al. 2008). Fok et al. (2012) detected maser emission from the wind of W26, but were not able to measure the wind velocity with their data. Advances in instrumental capabilities will soon enable measurements of the wind velocity using similar methods as for NML Cyg, so an independent constraint from nebular emission at larger radii is a useful check on these models.

4.2. Comparison with W26

We have shown that there is some agreement between observations and the predictions of the spherically symmetric photoionized wind described in Mackey et al. (2014), but that there are also important differences that seem difficult to reconcile with any spherically symmetric wind model. The whole photoionized nebula seems to be significantly blueshifted wherever there is significant emission, favouring an interpretation where W26 is on the near side of Westerlund 1 and its wind is being swept towards us by some process (either radiative or hydrodynamical) associated with winds, radiation, or supernovae from the star cluster.

This seems more similar to the one-sided photoionization of the wind of NML Cyg (by the Cyg OB2 association; Morris & Jura 1983) than to the diffuse ionizing radiation field that is invoked to explain the photoionization of Betelgeuse’s wind (Mackey et al. 2014). If the Mackey et al. (2014) model is correct, then Betelgeuse is irradiated by an isotropic ionizing radiation field, creating a spherically symmetric ionization front in the star’s wind and producing a static shell. Betelgeuse is in this relatively quiet environment because it is a runaway star that is now far from its place of birth. NML Cyg is very different, being near a massive and young OB association, and so experiencing a very anisotropic radiation field. W26 experiences an even more extreme environment because it is still within (at least in projection) the most massive star cluster in our Galaxy, that contains many WR and O-type stars. It follows then, that we should not be surprised if W26 were to experience an asymmetric radiation field, and if a multi-dimensional extension of the Mackey et al. (2014) model were required to explain the data. We discuss this and other alternative explanations in the following subsections.

4.3. Bow shock

If the nebula is a bow shock then it must be oriented away from us, i.e., W26 is embedded in an ISM with a bulk flow directed approximately towards us. The nebula is blueshifted by about 25 km s-1 relative to W26, which represents the velocity of the wind material that has been entrained in the bulk flow by the bow shock. This radial velocity then measures the relative velocity between the star and the ISM. The standoff distance, RSO, is then less than or equal to the smallest radial extent of the nebula, which we estimate to be 0.037 pc, and is given by (4)where ρ0 is the ISM density and v0 = 25 km s-1 the flow velocity of the ISM past W26. The wind is photoionized, so we assume this accelerates the ionized part of the wind to at least v = 35 km s-1 (Meyer et al. 2014a). We then use our upper and lower estimates of for W26 to estimate the value of ρ0 that is required to produce the observed nebula. We find ρ0 = 2.15 × 10-20 g cm-3 (or n0 = 104 cm-3, where n0 is the number density of H and He nuclei) for = 10-4 M yr-1, or 5 times smaller for = 2 × 10-5 M yr-1. This seems an unreasonably large ISM density for the intracluster medium of Westerlund 1 (X-ray emission suggests n0 ≈ 0.5 cm-3; see below), arguing against a straightforward bow shock interpretation for the nebula. For reference to the following discussion, the ram/thermal pressure required to confine the wind to RSO is p = 1.3 × 10-7 dyne cm-2 if = 10-4 M yr-1, or p = 2.7 × 10-8 dyne cm-2 if = 2 × 10-5 M yr-1, both calculated for v = 35 km s-1.

To further test this model, we relax the assumption that the shocked wind of W26 is accelerated to the full velocity of the bulk flow, which might be possible if there is little mixing at the contact discontinuity. This then allows a much faster bulk flow with a lower density, but to decrease the flow density to n0 ~ 1 cm-3 we require a gas velocity, v0 ≳ 103 km s-1. This sort of gas flow in the intracluster medium is only likely in the immediate vicinity of an extreme driver of momentum, e.g. a massive star wind or a supernova. In comparison to other massive star clusters, Povich et al. (2008) estimate ram pressures of dyne cm-2 for bow shocks around O stars in the M17 and RCW 49 massive-star-forming regions. This shows that the nebula around W26 is probably even more extreme than any of the bow shocks studied by Povich et al. (2008), if the bow shock interpretation is correct, suggesting again that perhaps we need an extreme source of momentum very nearby, rather than a cluster-scale gas outflow.

4.4. Thermal pressure confinement

Thermal pressure confinement is a promising mechanism, because X-ray observations have revealed a hot, diffuse, intracluster medium (Kavanagh et al. 2011). This is produced by the shocked wind of the WR and O stars in Westerlund 1, by the numerous supernovae that have already occurred, or a combination of both. These authors estimate the H number density and temperature of the X-ray emitting plasma to be nH = 0.5 cm-3 and kT = 3.7 keV (T = 4.3 × 107 K), respectively. For a fully ionized plasma containing 1/11 helium by number, this gives a thermal pressure of p0 = 2.2nHkT = 6.5 × 10-9 dyne cm-2. Balancing this against the ram pressure of the wind from W26 gives the radius of the wind termination shock, Rts, of (5)For our two mass-loss rates we obtain Rts = 0.17 pc ( = 10-4 M yr-1) and Rts = 0.075 pc ( = 2 × 10-5 M yr-1), both of which are larger than the observed extent of the nebula (0.037 pc). For the lower mass-loss rate, however, Rts is not much larger than the observed nebula, and a factor of 4 increase in the external pressure would bring the numbers into agreement (or a 4 × decrease in to 5 × 10-6 M yr-1, although this is unlikely for such a massive red supergiant). Some asymmetric pressure is still required, for example a subsonic bulk flow (e.g. Mackey et al. 2015), to create the asymmetric nebula, so purely thermal-pressure confinement is not a viable explanation.

4.5. Wind-wind interaction

The star W9 is 0.23 arcmin from W26 in projection, or 0.24 pc at the distance of Westerlund 1. It is the strongest radio emitter in Westerlund 1 and is surrounded by a dense wind with ≈ 3.3 × 10-4 M yr-1 and v ~ 200 km s-1 (Dougherty et al. 2010, model 3). If the winds of the two stars collide, then the wind of W26 will be shocked and decelerated at a distance, Rwc, from W26 that is given by (6)where R0 is the separation of the two stars and c is the ratio of the wind momenta: (7)For the W9 wind properties above, and using = 10-4 M yr-1 and v = 35 km s-1 for W26, we find Rwc = 0.19R0. For W26 with = 2 × 10-5 M yr-1, this reduces to Rwc = 0.093R0. The observed radius of the nebula is 0.037 pc, implying R0 = 0.20 pc for the larger estimate or R0 = 0.40 pc for the smaller . This is again comparable to the observed separation between the two stars, and so should be considered as a plausible alternative scenario to a simple photoionized wind.

Considering the pressure of the X-ray plasma quoted above (6.5 × 10-9 dyne cm-2), the wind of W9 would have a termination shock at Rts = 0.73 pc. The freely expanding wind region would have a diameter of 1.4 arcmin at the distance of Westerlund 1, covering a large fraction of the core of Westerlund 1. This means that W26 and many other stars in the cluster may well be embedded in the freely-expanding wind of W9, and could be experiencing a wind-wind collision.

4.6. Photoionized wind

We have shown that a spherically symmetric model for a photoionized wind has severe problems explaining the observations, but perhaps an asymmetric model could work. This was proposed by Morris & Jura (1983) for NML Cyg, to explain the arc-shaped H ii region on one side of the star. In this model there is one dominant radiation source so that the wind is photoionized mainly from one direction and the ionization front has a shape similar to a bow shock (for mathematical details, see Morris & Jura 1983). If the ionization front is D-type then the wind can be decelerated into a one-sided shell, and the asymmetry of the shell leads to non-radial flows away from the radiation source because of the rocket effect (Oort & Spitzer 1955). This requires multi-dimensional simulations to model quantitatively, beyond the scope of this paper, but we intend to investigate this in more detail in future work.

The EUV fluxes in Table 1 show what is required to ionize the wind of W26, and simple estimates show that nearby stars can provide these fluxes. The Wolf-Rayet star WR77o (WR B; Crowther et al. 2006) is about 0.5 pc from W26 (projected distance). According to Crowther (2007), it should have an ionizing photon luminosity, Q0 = 1049.4 s-1. For a distance of 0.5 pc, the EUV flux at W26 is Fγ = 8.4 × 1011 cm-2 s-1, comparable to the largest fluxes used in the spherically symmetric simulations. Its wind, by contrast, is probably not strong enough to shape the nebula, because Crowther (2007) estimate = 1.6 × 10-5 M yr-1 and v = 1300 km s-1, too weak to confine W26’s wind according to Eq. (6).

Much closer to W26, the O9 Iab star W25 has an ionizing photon luminosity of Q0 = (1048.8−1049.0) s-1 (Martins et al. 2005), and is only 0.09 pc from W26 in projection. This corresponds to an ionizing flux of Fγ = (6.4−10.0) × 1012 cm-2 s-1, and so could ionize the wind of W26 even deeper than the observed bright nebula, if the physical separation of the two stars is comparable to the projected separation. Another candidate is W9, discussed above in the context of wind-wind collision, but we do not have any reliable estimates of its ionizing photon luminosity.

4.7. Hot companion star?

Westerlund 1 was observed with Chandra in X-rays and analysis of point sources was presented by Skinner et al. (2006). W9 was clearly detected, which they argue is because W9 is a colliding-wind binary, but W26 (referred to as Ara A) was not detected although it was observed. If W26 was a binary system with a neutron star (or black hole) companion that provided the ionizing photons to produce the nebula, then it would also be bright in X-rays, and this is not observed. W26 is also too young to host a white dwarf companion, so we can exclude that a compact companion is ionizing the wind. It is also not a photometric variable star (Bonanos 2007) and, although it has spectral variations, these are consistent with other massive red supergiants and not indicative of duplicity (Clark et al. 2010).

4.8. Outlook

We see that there are many radiation sources that should photoionize the wind of W26 to at least the radius of the observed nebula, and potentially much closer to the star, although it is not yet clear if one-sided photoionization could produce the velocity signatures of the observed nebula. We have also shown that the measured X-ray plasma pressure could contribute to confining the wind. The large temperature difference between the photoionized red supergiant wind and the X-ray plasma of the intracluster medium may also set up a strong thermal conduction front at the contact discontinuity between wind and ISM, which has been shown to strongly modify the properties of O star bow shocks (Comerón & Kaper 1998; Meyer et al. 2014b). We suggest, however, that a wind-wind collision between the winds of W9 and W26 is the most obvious scenario that could produce the observed nebula. Naively, we expect that a bow shock model (possibly from a wind-wind collision with the wind of W9) is the only one that can explain why all of the bright emission is blueshifted.

It would be very useful to have higher resolution spectra, to reduce the uncertainties in the radial velocity of the spectral lines, and to search for substructure in the lines (e.g. redshifted absorption). Better spatial coverage of the nebula is also required. In particular spectra of the triangular nebula discovered by Wright et al. (2014) and the emission connecting it with the W26 nebula could reveal velocity gradients with distance from W26. This would show whether the gas is being accelerated away from W26, as expected for gas that is entrained by a wind. This would also provide more conclusive evidence for whether W26 is on the near or far side of Westerlund 1: if on the near side then the gas should be more blueshifted the further it is from W26.

5. Conclusions

We have presented the first predictions for nebular emission from photoionized winds of red supergiants, using spherically symmetric, radiation-hydrodymamical simulations and postprocessing radiative transfer. Although the results have been compared only to W26, they also have application to photoionized winds around other red supergiants. The simulations improve on the previous work of Mackey et al. (2014) in that the gas temperature is calculated self-consistently from mechanical and radiative heating/cooling, allowing us to make realistic predictions of temperature-sensitive nebular emission from ionized gas.

We produce position-velocity datasets for Hα and [N II] spectral lines from the simulations, and study the properties of this emission. The peak brightness of the simulated nebulae is 10-13−10-11 erg cm-2 s-1 arcsec-2, for mass-loss rates, = (0.2−1) × 10-4 M yr-1, and wind velocities, v = (15−30) km s-1. Simulations of slow winds that are decelerated into a dense shell show strongly limb-brightened line emission, and radial velocity of spectral lines is independent of the wind speed. Faster winds (≳ 22 km s-1) that cannot form a dense shell have less limb-brightening and the line radial velocity is a good tracer of the wind speed. All of the simulations have the most blue-shifted emission at the position of the star, and the radial velocity of the lines goes to zero at the region of brightest emission (the limb), as is typical for emission from an expanding shell.

If nitrogen is not enriched in the wind, then the predicted line ratio of [N II] to Hα is about 0.9. Nitrogen is enhanced in red supergiant winds of a factor of about 3 (cf. Brott et al. 2011), and for this we obtain line ratios of about 2.5−2.8.

We have compared spectra of the Hα and [N II] emission lines from the nebula around the red supergiant W26 in Westerlund 1 (Wright et al. 2014) with our simulations to investigate the origin and physical properties of the nebula. We also obtain a temperature for the star Teff = 3600 K from the same spectra. The brightness distribution of the emission and the line ratio of [N II] to Hα show similarities to the photoionization model, but there are discrepencies for the radial velocity of the emission. All of the observed bright nebular emission is blueshifted with respect to W26 by 1525 km s-1, and the brightest emission is the most blueshifted, neither of which can be reproduced by a spherically symmetric model of a photoionized wind. From the observed line ratio (2.25), we can conclude that the nebula is photoionized, however, and that it is enriched in nitrogen by a factor of about 2.35. This shows that the nebula is produced by the wind of W26, although the discrepencies with our model suggest that the outer edge of the nebula may be hydrodynamically confined rather than freely expanding.

We propose that a wind-wind collision scenario, where the wind of W26 is swept towards the observer by the stronger wind of the nearby star W9, may match the data better. Our results suggest that the circumstellar medium may have more structure even closer to the star than the observed nebula, because it does not appear as though we have detected the ionization front and its associated shell in the current observations. Multidimensional models and/or simulations, and spectra with higher resolution and better spatial coverage are necessary to progress further in understanding the circumstellar nebula around W26. Further study of this nebula will contribute to our understanding of what happens to the mass lost by red supergiants in star clusters. It is important to know whether it is rapidly heated and expelled from the cluster, or whether it can remain within the cluster for an extended period of time, remain cool and dense, and potentially be available for secondary star formation in the cluster.


1

Image Reduction and Analysis Facility (IRAF – http://iraf.noao.edu/) is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

2

http://marcs.astro.uu.se/

Acknowledgments

This project was supported by the Deutsche Forschungsgemeinschaft priority program 1573, Physics of the Interstellar Medium. L.F. acknowledges financial support from the Alexander von Humboldt foundation. The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUROPA at Jülich Supercomputing Centre (JSC). We thank I. Negueruela for radial velocity data for W26 and A. Whitworth for useful discussions on star formation within massive star clusters. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 087.D-0673(A). We are grateful to the referee for comments that significantly improved the presentation of our work.

References

Online material

Appendix A: Results for the ESE side of the nebula

Here we plot the line intensity, radial velocity, and line ratio for the simulations that are not discussed in the text. Figures A.1 and A.2 show M5V20E and M5V30E, respectively, and Figs. A.3A.6 show results for M4V15E, M4V20E, M4V25E, and M4V30E.

thumbnail Fig. A.1

Hα and [N II] spectral lines from the simulation M5V20E, with = 2 × 10-5 M yr-1, v = 20 km s-1, and Fγ = 2.78 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The solid lines show the simulation, and the dotted lines with points show the observations to the ESE of W26. The panels show a) line brightness; b) line ratio [N II]/Hα; and c) radial velocity of the peak emission, all as a function of distance from the star.

Open with DEXTER

thumbnail Fig. A.2

Hα and [N II] spectral lines from the simulation M5V30E, with = 2 × 10-5 M yr-1, v = 30 km s-1, and Fγ = 1.26 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER

thumbnail Fig. A.3

Hα and [N II] spectral lines from the simulation M4V15E, with = 10-4 M yr-1, v = 15 km s-1, and Fγ = 3.67 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER

thumbnail Fig. A.4

Hα and [N II] spectral lines from the simulation M4V20E, with = 10-4 M yr-1, v = 20 km s-1, and Fγ = 5.05 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER

thumbnail Fig. A.5

Hα and [N II] spectral lines from the simulation M4V25E, with = 10-4 M yr-1, v = 25 km s-1, and Fγ = 3.10 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER

thumbnail Fig. A.6

Hα and [N II] spectral lines from the simulation M4V30E, with = 10-4 M yr-1, v = 30 km s-1, and Fγ = 2.15 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER

Appendix B: Results for the WNW side of the nebula

Here we plot the line intensity, radial velocity, and line ratio for the simulations labelled “W” in Table 1. These model the nebula Westnorthwest (WNW) from W26, in the direction of the extreme sgB[e] star W9. As noted in the text, the nebula is confined closer to W26 on this side, possibly because of the proximity of W9. The “W” simulations have larger photon fluxes than the corresponding “E” simulations, so that the ionized region penetrates deeper into the wind of the red supergiant. The results of simulations are overplotted on the observed data in Fig. B.1 for simulation M5V15W, and similarly in Figs. B.2B.4 for M5V20W, M5V25W, and M5V30W. Figures B.5B.8 show results for the simulations M4V15W, M4V20W, M4V25W, and M4V30W, which have a denser wind because of a larger mass-loss rate.

The observed data shows no limb-brightened peak in the nebular emission, rather just a shoulder of emission between 0.02 and 0.03 arcmin. If we discard the uncertain two points closest to zero (at W26) then the emission decreases monotonically with distance from W26, in strong constrast to the predictions of simulations. This makes the discrepency with the limb-brightened emission from M5V15W (Fig. B.1a) even larger for this side of W26. The observed line ratio is slightly lower on this side of W26, but again is relatively constant across the nebula, decreasing only outside the bright part of the nebula (θ> 0.04 arcmin). Additionally, the radial velocity of the lines (Fig. B.1c) appears too low at all radii, but particularly for θ> 0.03 arcmin.

Simulation M5V25W (Fig. B.3) has better agreement with observations, for the same reasons as M5V25E: the ionization front emission is less peaked, the radial velocity is more blueshifted in agreement with observations, and the line ratio has reasonable agreement with observations. Again, however, the crucial problem is that the observed nebular emission is all blueshifted by about 25 km s-1, whereas in the simulation the brightest emission is neither blue- nor redshifted.

thumbnail Fig. B.1

Hα and [N II] spectral lines from the simulation M5V15W, after 0.1 Myr of evolution. The solid lines show the simulation, and the dotted lines with points show the observations to the WNW of W26. The panels show a) line brightness; b) line ratio [N II]/Hα; and c) radial velocity of the peak emission, all as a function of distance from the star. Vertical dotted lines in panel a) show the range of radii used to set the normalisation of the simulated and observed [N II] line brightness.

Open with DEXTER

thumbnail Fig. B.2

Hα and [N II] spectral lines from the simulation M5V20W, with = 2 × 10-5 M yr-1, v = 20 km s-1, and Fγ = 1.29 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

thumbnail Fig. B.3

Hα and [N II] spectral lines from the simulation M5V25W, with = 2 × 10-5 M yr-1, v = 25 km s-1, and Fγ = 8.29 × 1010 cm-2 s-1, again after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

thumbnail Fig. B.4

Hα and [N II] spectral lines from the simulation M5V30W, with = 2 × 10-5 M yr-1, v = 30 km s-1, and Fγ = 5.83 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

thumbnail Fig. B.5

Hα and [N II] spectral lines from the simulation M4V15W, with = 10-4 M yr-1, v = 15 km s-1, and Fγ = 1.70 × 1012 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

thumbnail Fig. B.6

Hα and [N II] spectral lines from the simulation M4V20W, with = 10-4 M yr-1, v = 20 km s-1, and Fγ = 2.34 × 1012 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

thumbnail Fig. B.7

Hα and [N II] spectral lines from the simulation M4V25W, with = 10-4 M yr-1, v = 25 km s-1, and Fγ = 1.44 × 1012 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

thumbnail Fig. B.8

Hα and [N II] spectral lines from the simulation M4V30W, with = 10-4 M yr-1, v = 30 km s-1, and Fγ = 9.95 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER

All Tables

Table 1

Parameters used for the 16 spherically symmetric radiation-hydrodynamics simulations.

All Figures

thumbnail Fig. 1

Comparison between observed FORS spectra for W26 (black solid line) and MARCS model fluxes (dashed red line) for a Teff = 3600 K supergiant. The wavelength gap where no spectrum was obtained arises from a gap between the two grisms.

Open with DEXTER
In the text
thumbnail Fig. 2

Spectral observations of W26 and neighbouring stars. a) Radial velocity of the Hα (black dots) and [N II] (blue triangles) lines along the slit in the W26 co-moving reference system (~−32 km s-1), with position measured relative to W26. The directions along the slit are also indicated: WNW for Westnorthwest and ESE for Eastsoutheast. b) Line brightness (intensity) in the two spectral lines along the slit (arbitrary units). c) Fitted Gaussian linewidth, σ in the two spectral lines along the slit. The three stars included in the slit, and no further away than 0.4 arcmin from W26, are labelled according to Negueruela et al. (2010).

Open with DEXTER
In the text
thumbnail Fig. 3

Radial structure of the wind in simulation M5V15E: we show the H number density nH, H+ fraction, gas radial velocity vr, and temperature T on a logarithmic scale as a function of distance from the red supergiant.

Open with DEXTER
In the text
thumbnail Fig. 4

As Fig. 3 but for simulation M5V20E.

Open with DEXTER
In the text
thumbnail Fig. 5

As Fig. 3 but for simulation M5V25E.

Open with DEXTER
In the text
thumbnail Fig. 6

Hα and [N II] spectral lines from the simulation M5V15E, with = 2 × 10-5 M yr-1, v = 15 km s-1, and Fγ = 1.70 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The dashed lines show the simulation, solid lines smoothed to 1 arcsec resolution at the distance of Westerlund 1, and the dotted lines with points show the observations to the ESE of W26. The panels show a) line brightness; b) line ratio [N II]/Hα; and c) radial velocity of the peak emission, all as a function of distance from the star. Vertical dotted lines in panel a) show the range of radii used to set the normalisation of the simulated and observed [N II] line brightness.

Open with DEXTER
In the text
thumbnail Fig. 7

Hα and [N II] spectral lines from the simulation M5V25E, again after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. 6.

Open with DEXTER
In the text
thumbnail Fig. 8

Comparison between wind velocity, v, and the blueshift of the Hα and [N II] spectral lines for lines of sight close to the star. The boundary between D-type and R-type ionization fronts is at v ≈ 22 km s-1. The red lines are for Hα and the blue for [N II].

Open with DEXTER
In the text
thumbnail Fig. 9

Comparison between wind velocity, v, and the degree of limb brightening of the Hα and [N II] spectral lines, defined as the ratio of the peak brightness to that at lines of sight near the star (r = 0). The boundary between D-type and R-type ionization fronts is at v ≈ 22 km s-1. The red lines are for Hα and the blue for [N II]. Panel a) shows the raw simulation data, and panel b) when smoothed to a spatial resolution of 1 arcsec. The lines in panel b) follow the legend in panel a).

Open with DEXTER
In the text
thumbnail Fig. A.1

Hα and [N II] spectral lines from the simulation M5V20E, with = 2 × 10-5 M yr-1, v = 20 km s-1, and Fγ = 2.78 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The solid lines show the simulation, and the dotted lines with points show the observations to the ESE of W26. The panels show a) line brightness; b) line ratio [N II]/Hα; and c) radial velocity of the peak emission, all as a function of distance from the star.

Open with DEXTER
In the text
thumbnail Fig. A.2

Hα and [N II] spectral lines from the simulation M5V30E, with = 2 × 10-5 M yr-1, v = 30 km s-1, and Fγ = 1.26 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER
In the text
thumbnail Fig. A.3

Hα and [N II] spectral lines from the simulation M4V15E, with = 10-4 M yr-1, v = 15 km s-1, and Fγ = 3.67 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER
In the text
thumbnail Fig. A.4

Hα and [N II] spectral lines from the simulation M4V20E, with = 10-4 M yr-1, v = 20 km s-1, and Fγ = 5.05 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER
In the text
thumbnail Fig. A.5

Hα and [N II] spectral lines from the simulation M4V25E, with = 10-4 M yr-1, v = 25 km s-1, and Fγ = 3.10 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER
In the text
thumbnail Fig. A.6

Hα and [N II] spectral lines from the simulation M4V30E, with = 10-4 M yr-1, v = 30 km s-1, and Fγ = 2.15 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. A.1.

Open with DEXTER
In the text
thumbnail Fig. B.1

Hα and [N II] spectral lines from the simulation M5V15W, after 0.1 Myr of evolution. The solid lines show the simulation, and the dotted lines with points show the observations to the WNW of W26. The panels show a) line brightness; b) line ratio [N II]/Hα; and c) radial velocity of the peak emission, all as a function of distance from the star. Vertical dotted lines in panel a) show the range of radii used to set the normalisation of the simulated and observed [N II] line brightness.

Open with DEXTER
In the text
thumbnail Fig. B.2

Hα and [N II] spectral lines from the simulation M5V20W, with = 2 × 10-5 M yr-1, v = 20 km s-1, and Fγ = 1.29 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text
thumbnail Fig. B.3

Hα and [N II] spectral lines from the simulation M5V25W, with = 2 × 10-5 M yr-1, v = 25 km s-1, and Fγ = 8.29 × 1010 cm-2 s-1, again after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text
thumbnail Fig. B.4

Hα and [N II] spectral lines from the simulation M5V30W, with = 2 × 10-5 M yr-1, v = 30 km s-1, and Fγ = 5.83 × 1010 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text
thumbnail Fig. B.5

Hα and [N II] spectral lines from the simulation M4V15W, with = 10-4 M yr-1, v = 15 km s-1, and Fγ = 1.70 × 1012 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text
thumbnail Fig. B.6

Hα and [N II] spectral lines from the simulation M4V20W, with = 10-4 M yr-1, v = 20 km s-1, and Fγ = 2.34 × 1012 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text
thumbnail Fig. B.7

Hα and [N II] spectral lines from the simulation M4V25W, with = 10-4 M yr-1, v = 25 km s-1, and Fγ = 1.44 × 1012 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text
thumbnail Fig. B.8

Hα and [N II] spectral lines from the simulation M4V30W, with = 10-4 M yr-1, v = 30 km s-1, and Fγ = 9.95 × 1011 cm-2 s-1, after 0.1 Myr of evolution. The lines and symbols are the same as in Fig. B.1.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.