Proper motions of five OB stars with candidate dusty bow shocks in the Carina Nebula

TL;DR

This study measures the proper motions of five OB stars with candidate bow shocks in the Carina Nebula using HST data, revealing that these shocks are often shaped by ambient gas flows rather than stellar motion, and suggesting in-situ star formation.

Contribution

First direct proper motion measurements of OB stars with candidate bow shocks in Carina, clarifying the origins of observed bow shocks and the formation history of the region.

Findings

Stars are not runaway stars; their motions do not point away from the nebula's center.

Large-scale gas flows influence the shape of bow shocks more than stellar motion.

Stars' motions suggest in-situ formation rather than migration from central clusters.

Abstract

We constrain the proper motions of five OB stars associated with candidate stellar wind bow shocks in the Carina Nebula using HST ACS imaging over 9--10 year baselines. These proper motions allow us to directly compare each star's motion to the orientation of its candidate bow shock. Although these stars are saturated in our imaging, we assess their motion by the shifts required to minimize residuals in their Airy rings. The results limit the direction of each star's motion to sectors less than 90 degrees wide. None of the five stars are moving away from the Carina Nebula's central clusters as runaway stars would be, confirming that a candidate bow shock is not necessarily indicative of a runaway star. Two of the five stars are moving tangentially relative to the orientation of their candidate bow shocks, both of which point at the OB cluster Trumpler 14. In these cases, the large-scale flow of the interstellar medium, powered by feedback from the cluster, appears to dominate over the motion of the star in producing the observed candidate bow shock. The remaining three stars all have some component of motion toward the central clusters, meaning that we cannot distinguish whether their candidate bow shocks are indicators of stellar motion, of the flow of ambient gas, or of density gradients in their surroundings. In addition, these stars' lack of outward motion hints that the distributed massive-star population in Carina's South Pillars region formed in place, rather than migrating out from the association's central clusters.

Tables

Table 1. Table 1: HST data log.

ERO	Star ID	R.A.	Dec.	Spectral	V	ACS	Date 1	Exp. time 1	Date 2	Exp. time 2
No.a		(J2000)	(J2000)	type	(mag)	field		(s)		(s)
2	ALS 15206	10:44:00.9	-59:35:46	O9.2 Vb	10.7d	TR14	2005 Jul 17	2 $\times$ 500	2015 Jun 28	2 $\times$ 520
23	TYC 8626-2506-1	10:44:30.2	-59:26:13	O9 Vb	10.9e	TR14	2005 Jul 17	2 $\times$ 500	2015 Jun 28	2 $\times$ 520
24	CPD-59 2605	10:44:50.4	-59:55:45	B1 Vc	11.1f	POS27	2006 Mar 18	2 $\times$ 500	2015 Mar 12	2 $\times$ 560
25	HDE 305533	10:45:13.4	-59:57:54	B1 Va	10.6f	POS26	2006 Mar 16	2 $\times$ 500	2015 Mar 12	2 $\times$ 560
31	HD 93576	10:46:53.8	-60:04:42	O9.5 IVb	9.6d	POS20	2006 Mar 15	2 $\times$ 500	2015 Mar 11	2 $\times$ 455

Table 2. Table 2: Local proper motions of stars associated with candidate bow shocks.

ERO No.	Star ID	$\delta x$	$\delta y$	${\mu }_{\alpha }\cos \delta$	${\mu }_{\delta }$	v${}_{\text{T}}$a	Position angle
		(pixels)	(pixels)	(mas yr-1)	(mas yr-1)	(km s-1)	(deg E of N)
2	ALS 15206	-0.39 (0.16)	0.35 (0.28)	2.0 (0.8)	1.8 (1.4)	29 (17)	48 (25)
23	TYC 8626-2506-1	-0.27 (0.29)	-0.45 (0.39)	1.4 (1.4)	-2.3 (2.0)	29 (27)	149 (35)
24	CPD-59 2605	0.36 (0.37)	0.46 (0.20)	-2.0 (2.1)	2.5 (1.1)	35 (25)	322 (31)
25	HDE 305533	-0.39 (0.18)	0.40 (0.19)	2.2 (1.0)	2.2 (1.1)	34 (16)	45 (19)
31	HD 93576	0.26 (0.04)	0.04 (0.07)	-1.4 (0.2)	0.2 (0.4)	16 ( 5)	279 (15)

Full text

††pubyear: 2017††pagerange: Proper motions of five OB stars with candidate dusty bow shocks in the Carina Nebula–Proper motions of five OB stars with candidate dusty bow shocks in the Carina Nebula

1 Introduction

Feedback from massive stars impacts their surroundings on scales ranging from the shaping of their immediate circumstellar environment to the reionization of the universe. Stellar wind bow shocks, falling on the former end of that scale, provide important information about a star’s history and environment. Bow shocks are produced when the relative velocity between a star and the surrounding interstellar medium (ISM) is supersonic (Baranov et al., 1971; van Buren & McCray, 1988). They typically appear as arc-shaped features in optical line emission (e.g., Kaper et al., 1997; Bally et al., 2000; Brown & Bomans, 2005; Brownsberger & Romani, 2014) and/or thermal infrared continuum emission from dust (e.g., van Buren et al., 1995; Noriega-Crespo et al., 1997; Comerón & Pasquali, 2007; France et al., 2007; Gáspár et al., 2008; Peri et al., 2012; Winston et al., 2012; Kobulnicky et al., 2016). These features mark the sweeping-up of ambient material between the stellar wind termination shock and a second shock from the supersonic motion. The orientation of the bow shock arc depends on the direction of the relative motion, although it can be skewed by density gradients in the environment (Wilkin, 2000). The arc’s standoff distance from the star depends on the pressure balance between the stellar wind and the ISM and hence on the magnitude of their relative motion and the density of the ISM. Bow-shock-like structures may also be produced where dust in a photoevaporative flow is stalled by radiation pressure rather than the stellar wind (Ochsendorf et al., 2014a; Ochsendorf et al., 2014b; Ochsendorf & Tielens, 2015; Ochsendorf et al., 2015). The asymmetric stellar wind bubbles of slower-moving stars may also have a similar appearance in the mid-infrared (Mackey et al., 2015, 2016).

Arc-shaped structures around massive stars have commonly been considered a marker of high stellar velocities, under the assumption that the relative motion between star and ISM is dominated by the absolute motion of the star (van Buren et al., 1995; Kaper et al., 1997; Gvaramadze & Bomans, 2008; Gvaramadze et al., 2010; Kobulnicky et al., 2010; Gvaramadze et al., 2011a, b). The typical velocity of an O-type star relative to its surroundings is $\sim 10$ km s*-1* (Blaauw, 1961; Cruz-González et al., 1974; Gies & Bolton, 1986; Tetzlaff et al., 2011), comparable to the speed of sound in an H ii region, but 20–30% of O-type stars are “runaways” with velocities $\gtrsim 40$ km s*-1* (Blaauw, 1961; Cruz-González et al., 1974; Stone, 1991; Tetzlaff et al., 2011). The high speeds of runaway stars are imparted through dynamical interactions in a cluster (Poveda et al., 1967; Gies & Bolton, 1986; Fujii & Portegies Zwart, 2011), through the disruption of a binary system when the companion star explodes as a supernova (Blaauw, 1961), or through a two-step scenario involving both processes (Pflamm-Altenburg & Kroupa, 2010). Runaways make up 50–100% of field O-type stars, the O-type population found outside clusters and associations (de Wit et al., 2005; Schilbach & Röser, 2008; Gvaramadze et al., 2012).

The question of whether all field O-type stars are runaways, or whether a small fraction formed in isolation, is of key importance to our understanding of massive star formation. The monolithic collapse model (McKee & Tan, 2003; Krumholz et al., 2005; Krumholz et al., 2009) permits truly isolated massive star formation, albeit rarely, while the competitive accretion model (Zinnecker, 1982; Bonnell et al., 2001a, b; Bonnell et al., 2004) requires that massive stars form exclusively in clusters. In observational studies, the presence of a bow shock or candidate bow shock is sometimes taken as a clue that a given massive field star did not form in situ. For example, HD 48229 and HD 165319 were part of the $4\pm 2$% of all O-type stars identified by de Wit et al. (2004, 2005) as likely candidates for isolated massive star formation. Bow shocks were later discovered around both sources (Gvaramadze & Bomans, 2008; Gvaramadze et al., 2012), calling into question their origins in the field.

While 70% of bow shocks and bow-shock-like structures are located in relatively isolated environments consistent with runaway stars (Kobulnicky et al., 2016), the rest are found around OB stars in clusters and associations. These stars have sometimes been interpreted as runaway interlopers from other regions (Gvaramadze et al., 2011a). However, the assumption that the motion of the ambient ISM is negligible relative to that of the star may not always be valid, particularly in and around giant H ii regions. In many cases, bow shock orientations suggest that feedback-driven ISM flows are relevant. Povich et al. (2008) observed that bow shocks in the massive star-forming regions M17 and RCW 49 are oriented toward those regions’ central clusters, suggesting that global expansion of the H ii regions is the dominant component of the relative star–ISM velocity. Similarly, several bow shocks in Cygnus OB2 point toward the association’s interior (Kobulnicky et al., 2010), as do more than half of the candidate bow shocks in the Carina Nebula (Smith et al., 2010b; Sexton et al., 2015). The Galactic Plane survey of Kobulnicky et al. (2016) found that roughly 15% of infrared bow shocks are pointed at H ii regions, while another 8% face bright-rimmed clouds; they also noted that bow shock orientations are correlated on small scales, indicative of the influence of external forces. Povich et al. (2008) refer to such feedback-facing bow shocks as “interstellar weather vanes,” tracing photoevaporative flows off local dense gas and/or large-scale gas motions driven by cluster feedback. Kobulnicky et al. (2016) call them “in-situ bow shocks,” reflecting their origin around presumably non-runaway OB stars.

When the motion of the star dominates over the motion of the surrounding ISM, as it does for runaway stars, the bow shock is expected to point in the direction of the star’s motion. van Buren et al. (1995) surveyed bow shocks around known runaway stars and found that the bow shocks were preferentially aligned with their host stars’ proper-motion vectors. However, they used proper motions measured in an absolute reference frame, not corrected for Galactic rotation and solar peculiar motion and thus not necessarily representative of a star’s motion relative to the surrounding ISM. More recent surveys by Peri et al. (2012); Peri et al. (2015), again of bow shocks around known runaways stars, did correct proper motions for Galactic rotation and noted a similar, albeit qualitative, tendency for alignment. Individual runaway stars are also often observed to be moving in the direction of their bow shocks (e.g., Moffat et al., 1998, 1999; Comerón & Pasquali, 2007). But what about bow shock around stars that have not already been identified as runaways? Kobulnicky et al. (2016) compiled a sample of bow shocks without any selection on their host stars’ kinematics. They found that more than 50% of the host stars with significant measured proper motions had velocity–bow shock misalignments of more than 45°, although again, they were working with absolute proper motions rather than local. The relationship between stellar motion and bow shock orientation for stars in clusters and associations remains largely unexplored.

To further investigate this relationship, we measure local proper motions for five massive stars in the Carina Nebula (listed in Table 1), each of which is associated with a candidate bow shock from Sexton et al. (2015, which includes objects first identified by ). Smith et al. (2010b) and Sexton et al. (2015) identified a total of 39 “extended red objects” (EROs) in the Carina Nebula. These EROs exhibit extended, often arc-shaped, morphology in Spitzer Infrared Array Camera (IRAC) 8.0 $\mu$m images. Nine of the Sexton et al. (2015) EROs are clearly resolved arcs and are classified as morphological bow shock candidates; one of our stars (ALS 15206) is associated with one of these sources (ERO 2). Another eight of the Sexton et al. (2015) EROs lack resolved morphologies at 8.0 $\mu$m but have infrared colors that rule out emission from young stellar objects (YSOs) and polycyclic aromatic hydrocarbons (PAHs). The remaining four of our stars are associated with sources in this category, known as color bow shock candidates.

Our five target stars reside in the Carina Nebula: their visual magnitudes (see Table 1), spectral types, and extinctions (Povich et al., 2011) confirm that they are unlikely to be foreground or background objects. The Carina Nebula is home to nearly 70 O-type and evolved massive stars (Smith, 2006a), including some of the earliest known O-type stars (Walborn et al., 2002a) and the luminous blue variable $\eta$ Carinae (Davidson & Humphreys, 1997). At 2.3 kpc (Smith, 2006b), it is one of the closest and least-extincted massive star-forming regions. Its two central clusters, Trumpler (Tr) 14 and Tr 16, contain about half of its massive-star population. The rest is spread across $\sim$30 pc, mostly in a region of ongoing star formation known as the South Pillars (Smith et al., 2000). Emission-line profiles show that feedback from the central clusters is driving the expansion of multiple shells of ionized gas (Damiani et al., 2016), resulting in a global expansion of the H ii region at $\pm$15–20 km s*-1* (Walborn & Hesser, 1975; Walborn et al., 2002b, 2007). It is easy to envision that the inward-facing orientations of many of Carina’s candidate bow shocks are the result of this supersonic, feedback-driven ISM expansion (Sexton et al., 2015) or that they are shaped by interaction with dense photoevaporative flows. Here, we explore whether those interpretations are valid and to what degree these bow shocks are shaped by the motion and structure of the ISM versus the motion of their driving stars.

The organization of this paper is as follows: In Section 2, we describe our multiepoch Hubble Space Telescope (HST) observations, our image alignment procedure, and our method for measuring proper motions. We present our results and compare the stellar motions to the orientations of their associated bow shock candidates in Section 3. Section 4 discusses the implications and limitations of our results, and Section 5 summarizes our conclusions.

2 Observations and Analysis

2.1 HST ACS Imaging

We have conducted a large-scale multiepoch survey of the Carina Nebula using the Wide Field Camera (WFC) of HST’s Advanced Camera for Surveys (ACS). All observations were made with the F658N filter, which captures emission from H$\alpha$ and [N ii] $\lambda$6584. Our imaging coverage is shown in Figure 1, where each small rectangle is one orbit made up of three overlapping pairs of CR-SPLIT exposures. Orbital pointings were designed to target features of particular interest in star formation (pillars, Herbig-Haro objects, etc.) as well as the central Tr 14 and Tr 16 clusters. The pointings in Figure 1 are labelled according to their designations in the HST data archive.

The first epoch of our ACS observations was taken in 2005–2006 (GO-10241 and GO-10475, PI: N. Smith; see Smith et al. 2010a). The same set of observations was repeated in 2014–2015 (GO-13390 and GO-13791, PI: N. Smith). This second epoch was designed to replicate the first as closely as possible in pointing and position angle in order to minimize position-dependent systemic errors when measuring proper motions. Owing to changes in the HST Guide Star Catalog between epochs, we were unable to duplicate the orientation angle of the central segment of the Tr 14 mosaic (marked in orange in Figure 1) and of Positions 25 and 30. Those observations were rotated by $\sim$180°.

Also marked in Figure 1 are the locations of EROs from Sexton et al. (2015). Our survey serendipitously imaged the stars associated with seven EROs. However, as discussed below, we were unable to constrain the proper motions of the two ERO-associated stars in the central part of the Tr 14 mosaic (the part for which the orientation angle changed beween epochs), leaving us with a sample of five. Full details of the observations of each of these five stars are given in Table 1.

2.2 Image Alignment and Stacking

Our image alignment procedure, which adapts the methods of Anderson et al. (2008a, b), Anderson & van der Marel (2010), and Sohn et al. (2012), is described in detail in Reiter et al. (2015a, b) and Kiminki et al. (2016). In summary, we find the positions of uncrowded, unsaturated stars in individual exposures and use those positions to relate each image to a master, distortion-free reference frame. We use the program img2xym_WFC.09x10 (Anderson & King, 2006), which uses an array of effective point spread functions (PSFs) and has the option to fit a spatially constant perturbation PSF to account for telescope breathing and other focus changes. The measured stellar positions were then corrected for geometric distortion (Anderson, 2006).

A master reference frame with a pixel scale of 50 mas was constructed for each orbital pointing, aligned with north in the $+y$ direction. The six overlapping images from each epoch of a given pointing were stacked into two reference-frame master images (one per epoch) using the stacking algorithm of Anderson et al. (2008a). Object positions in the master images are directly comparable between epochs to an alignment accuracy of approximately 1 mas ($\sim$ 1 km s*-1* over a 9–10 year baseline at the distance of the Carina Nebula). We found that including a perturbation PSF in fitting stellar positions did not improve the alignment precision, but we address other possible effects of HST focus changes in Section 2.3.

In all cases, the master reference frames are not tied to an absolute proper-motion zero point. Instead, the zero point is based on the average motion of several hundred well-measured stars in the image. In other words, the bulk motion of the Carina Nebula is removed, as are smaller differences in the large-scale motion of Carina’s clusters and subclusters. Features that are locally stationary, like bow shocks, are expected to be stationary in our reference frames, allowing direct measurement of the motion of stars relative to their surroundings.

2.3 Measuring Local Proper Motions of Saturated Stars

With the images from two epochs on the same reference frame, measuring local proper motions for unsaturated stars in our stacked images is as simple as comparing their PSF-derived positions between epochs. However, most of the OB stars observed, including those associated with candidate bow shocks, are saturated in our ACS images (which were all $\sim$500 s long). We were unable to reconstruct the PSF core to perform traditional astrometry. Instead, we used the positions of the extended Airy rings, which are clearly visible for these stars in these deep, high-resolution images. The left column of Figure 2 shows the first-epoch image of all five ERO-associated stars for which we measured proper motions.

As marked in Figure 2, we identified four regions in the outer PSF of each star, avoiding diffraction spikes, saturation bleeding, and close companions. We then found the pixel offset in $x$ (west-east) and $y$ (north-south) by which the first epoch needed to be shifted in order to minimize the sum over those four regions of the absolute value of the flux difference between epochs. The best-fit offset was computed using the AMOEBA algorithm, the IDL implementation of the downhill simplex function minimization method (Nelder & Mead, 1965; Press et al., 1992). AMOEBA requires an input estimate. As recommended by Press et al. (1992), we run the algorithm twice, giving it a random starting estimate on the first run and then starting the second run at the best-fit parameters of the first. The resulting best-fit difference images (unshifted second epoch minus best-fit shifted first epoch) are shown in the middle column of Figure 2. Although the Airy rings do not disappear completely in the difference images, there are no systematic differences in residual flux between quadrants. The best-fit offset for each star, in pixels over the 9–10 year baseline, is plotted in the right column of Figure 2.

This method does not produce accurate results for saturated stars in fields that experienced significant rotation between epochs. The asymmetry of the ACS WFC PSF (Anderson & King, 2006; Mahmud & Anderson, 2008) causes the flux distribution in the outer PSF to be orientation-dependent. When the fields are rotated into alignment, the asymmetric flux distribution introduces an apparent shift of up to several pixels. Consequently, as mentioned in Section 2.1, we were unable to measure the true shifts of the two ERO-associated stars in the central part of our Tr 14 mosaic (see Figure 1). These stars were removed from further analysis and are not shown in Figure 2; our final sample consists of the five ERO-associated stars listed in Table 1.

To characterize the uncertainties in our fits, we ran several different tests. First, we adjusted the size and placement of the boxes used to calculate the residuals, then refit. These adjustments proved to have a negligible effect on the resulting best-fit offset. Second, we repeated the full fitting process 100 times and measured the standard deviation among the results. These ranged from 0.03 to 0.4 pixels depending on the star. Finally, we applied 100 random [$x$,$y$] offsets to the first-epoch images and repeated the full fitting process again for each, to determine how well we could recover the expected (artificial $+$ true) offsets. The standard deviation in the offsets recovered via this approach ranged from 0.04 to 0.23 pixels. For each star, we adopt the greater of the two uncertainty values as the formal uncertainty.

Changes in HST focus from both short-term thermal breathing and long-term non-thermal effects (e.g., Cox & Lallo, 2012) could induce an apparent offset between epochs by redistributing flux in the PSF. To evaluate the magnitude of this effect, we downloaded a Tiny Tim model PSF (Krist et al., 2011) for each star, at its observed chip position, for the appropriate focus value taken from the HST focus model (di Nino et al., 2008; Niemi & Lallo, 2010; Cox & Niemi, 2011). The true shift between “epochs” of Tiny Tim models is zero, so any measured shift would be a false positive. We ran the pair of Tiny Tim models for each star through our fitting procedure and measured apparent offsets of 0.035 to 0.22 pixels, which are illustrated by the shaded gray regions in the third column of Figure 2. In most cases, the focus-induced shifts are small and/or distributed roughly evenly about the origin. For HDE 305533 (ERO 24), however, the focus changes induced a systematic $-x,+y$ offset. Removing this shift would reduce the magnitude of the observed proper motion of HDE 305533 by roughly half, but would have only a small effect on its direction of motion.

3 Results

As the plots in the right column of Figure 2 demonstrate, each of the five stars associated with candidate bow shocks traveled no more than $\sim$0.5 pixels (25 mas) in any direction over their 9–10 year baselines. The measured pixel offsets are given in Table 2 along with the corresponding proper motion components, the total transverse velocity, and the position angle of the proper motion vector. The best-fit local transverse velocities range from 16 to 35 km s*-1*; the red arrows in Figure 1 show the expected travel distances over 105 yr. However, the uncertainties on most of the measured velocities are relatively large: most of the stars have motion consistent with zero within 1–2 $\sigma$. Only HD 93576 has motion significant at the 3$\sigma$ level, in the $x$ direction, although it has negligible $y$ (north–south) motion. We argue in Section 4.1 below that the true proper moions are likely on the smaller side of the allowed ranges. Even so, the results for all five stars constrain their directions of motion to sectors less than 90° wide.

In Figure 4, we compare the local proper motions of the stars to the orientations of their associated candidate bow shocks. The latter were determined by Sexton et al. (2015) based on the peaks of the 8.0 $\micron$ flux. (ERO 25, associated with HDE 305533, does not have a measured orientation.) In these three-color Spitzer IRAC images, the candidate stellar wind bow shocks appear as extended red (8.0 $\micron$) features, while nearby stars are prominent in blue (3.6 $\micron$) and green (4.5 $\micron$). We indicate the stars’ motions with white arrows (lengths arbitrarily scaled for visiblity) and show the range of possible directions with dotted white lines. The orientations of the candidate bow shocks, where known, are denoted by cyan arrows, and the outer yellow arrows show the directions to the various OB clusters.

The uppermost panels in Figure 4 show ALS 15206 (with ERO 2 from Sexton et al. 2015) and TYC 8626-2506-1 (ERO 23). Both of these stars are closer to Tr 14 than to Tr 16 (see Figure 1), both are associated with candidate bow shocks pointing at Tr 14, and both have proper motions directed tangentially to the orientation of their candidate bow shocks. The radial velocity of ALS 15206 is poorly constrained, as it is a probable spectroscopic binary (Kiminki et al. in preparation), but is consistent with being drawn from the radial velocity distribution of Tr 14 (Penny et al. 1993; García et al. 1998; Kiminki et al. in preparation). No radial velocity data exist for TYC 8626-2506-1. Thus based on the proper motions of their associated stars, the relative motion shaping EROs 2 and 23 appears to be dominated by the motion of the surrounding ISM, expanding outward from Tr 14. Unseen density gradients may also play a role, but the motion of the stars themselves do not look to be influencing the directions of these candidate bow shocks. They may truly be acting as “weather vanes,” tracing the large-scale flows of the ISM.

The middle row of Figure 4 shows CPD-59 2605 (ERO 24) and HDE 305533 (ERO 25). These stars’ candidate bow shocks are not arc-shaped at IRAC resolutions; Sexton et al. (2015) were able to measure an orientation for ERO 24 but not for ERO 25. Both stars are in Carina’s South Pillars region, and are roughly 7.5 pc northeast of the nominal center (Wu et al., 2009) of the sparse open cluster Cr 228. ERO 24 points north toward Tr 16, suggesting that it is influenced by feedback-driven outflows. However, its associated star (CPD-59 2605) has a local proper motion to the northwest, consistent with the orientation of the candidate bow shock within the uncertainties. It is thus not possible to distinguish between the effects of ISM flows and stellar motion in the case of ERO 24, as both may be relevant to shaping that feature. No radial velocity data exist for CPD-59 2605.

Relative to its surroundings, HDE 305533 (ERO 25) is moving to the northeast, away from the WNH star HD 93131 and the small group of late O-type stars that make up the center of Cr 228. Its path hints at an ejection from Cr 228, although at its observed speed it would have covered the 7.5 pc from Cr 228 in just 220,000 yr (but see discussion in Section 4.1 below on the likelihood that our measured proper motions are upper limits). The age and extent of Cr 228 are also poorly constrained, as it has often been considered an extension of Tr 16 (Walborn, 1995; Smith & Brooks, 2008), while X-ray data show it to be a discrete collection of groups and subclusters without a clear center (Feigelson et al., 2011). The origin of HDE 305533 is therefore not clearly evident. Its radial velocity (-18 km s*-1*; Levato et al., 1990) is typical for the massive stars in Cr 228 and the South Pillars region (Levato et al. 1990; Kiminki et al. in preparation) and comparable to the radial velocity of the surrounding gas pillars (Rebolledo et al., 2016).

Finally, the bottom panel of Figure 4 shows HD 93576, the binary system (Levato et al., 1990) associated with ERO 31. HD 93576 lies on the outskirts of the small open cluster Bochum 11 (Bo 11), located in the southeastern part of the South Pillars. Bo 11 is home to an estimated 1000 stars (Dias et al., 2002), including the O5 supergiant HD 93632 (Sota et al., 2014) and three additional O-type stars besides HD 93576. Photometric analysis indicates that the cluster is 3–5 Myr old (Fitzgerald & Mehta, 1987; Patat & Carraro, 2001; Preibisch et al., 2011); the presence of an O5I star suggests that 3 Myr is more likely. As Figure 4 shows, the proper motion vector of HD 93576 is closely aligned with the orientation of its candidate bow shock, which in turn points nearly directly away from the center of Bo 11. This configuration suggests that HD 93576 was ejected from Bo 11 and that its subsequent supersonic motion produced the observed candidate bow shock. Its systemic radial velocity (-8 km s*-1*; Kiminki et al. in preparation) is commensurate with the radial velocities of the other massive members of Bo 11 (Levato et al. 1990; Kiminki et al. in preparation) and the nearby dense gas (Rebolledo et al., 2016). But its observed proper motion (15 km s*-1*) and current position (1.9 pc from the center of Bo 11) indicate an ejection date just 130,000 yr ago. Perhaps HD 93576 was ejected 2–3 Myr after the formation of Bo 11 (possible; see Oh & Kroupa, 2016). Alternately, it could have originated outside the Carina Nebula and have a coincidental agreement in radial velocities, or the magnitude of its proper motion could be smaller than measured (see discussion in Section 4.1). In addition, its candidate bow shock is also generally directed toward the interior of the Carina Nebula, so a contribution from ISM flows driven by cluster feedback cannot be ruled out regardless of the origin of the star.

4 Discussion

4.1 Proper Motions as Upper Limits

For four of the five candidate bow shock host stars in our sample, we measure local proper motions of $\sim$30 km s*-1*, with associated uncertainties of $>15$ km s*-1*. (This total includes HDE 305533, whose observed motion may include a contribution from focus changes as described in Section 2.3.) Several lines of reasoning support the interpretation of these measured velocities as upper limits, with the true proper motions lying on the small side of the allowed range. First, the typical velocity of an O-type star relative to its surrounding is $\sim$ 10 km s*-1* (Blaauw, 1961; Cruz-González et al., 1974; Gies & Bolton, 1986; Tetzlaff et al., 2011). Of course, the stars in our sample are arguably not typical, given their association with candidate dusty bow shocks. Space velocities of 30 km s*-1* may qualify them as runaway stars, depending on the choice of runaway classification criteria. None of the five stars measured here are moving with trajectories that could have originated in Tr 14 or Tr 16, although HDE 305533 (ERO 25) and HD 93576 (ERO 31) may have come from the smaller open clusters Cr 228 and Bo 11, respectively. An object moving at 30 km s*-1* would cover 60 pc in 2 Myr (the average estimated age of Tr 14/16; Walborn, 1995; Smith, 2006a; Preibisch et al., 2011), and these stars are all significantly closer than that to any possible clusters of origin in the Carina Nebula. It is possible that all four of the stars with measured proper motions of $\sim$30 km s*-1* were ejected more recently, but that scenario would still not explain their directions of motion. Similarly, it is possible that all four are interlopers in the Carina Nebula, originating from another cluster, but the chance of encountering four such stars in our small sample is low. And as described in Section 3, the radial velocities of our sample stars, where available, agree with the radial motions of the surrounding stars and gas, consistent with more local origins.

In addition, speeds of 30 km s*-1* are inconsistent with the relative star–ISM velocities computed for Carina’s candidate bow shocks by Sexton et al. (2015). The pressure balance governing a standard bow shock makes it possible to estimate the relative star–ISM velocity as a function of measured standoff distance by making reasonable assumptions about stellar wind velocity, mass-loss rate, and ISM density. Sexton et al. (2015) measured the standoff distances of nine EROs in the Carina Nebula and found an average star–ISM velocity of 17 km s*-1*. For ERO 2 (associated with ALS 15206), the relative star–ISM velocity was a barely-supersonic 7 km s*-1*. Similar relative velocities for bow shocks in the massive star-forming region RCW 38 were reported by Winston et al. (2012). These numbers have substantial uncertainties due to the assumptions that go into their calculation, but they still suggest somewhat lower stellar velocities. Consider ERO 2 (ALS 15206): The orientation of the candidate bow shock indicates that the direction of the highest relative star–ISM velocity is to the northwest. We have measured that ALS 15206 is moving to the northeast, tangential to its candidate bow shock. If the relative star–ISM velocity in the direction of the candidate bow shock is on the order of 7 km s*-1*, the relative star–ISM velocity in a different direction cannot be substantially higher than that, although the picture may be complicated if there are density gradients in the ISM.

For these reasons, it is unlikely that the measured stars associated with candidate bow shocks are moving as fast as 30 km s*-1* relative to their surroundings. The local proper motions given here should thus be treated as upper limits. HD 93576 may be an exception, as its westward motion is measured at 3$\sigma$ significance (but this raises questions about its possible origin in Bo 11, as discussed in Section 3).

4.2 Comparison to Absolute Proper Motions

All five of the stars in our sample have proper motions listed in the USNO CCD Astrograph Catalog (UCAC4; Zacharias et al., 2013), and ALS 15206, TYC 8626-2506-1, and HD 93576 also have proper motions in the Tycho-2 Catalogue (Høg et al., 2000) and Gaia Data Release 1 (DR1; Gaia Collaboration et al., 2016a, b; Lindegren et al., 2016). (Note that Tycho-2 and Gaia DR 1 are not wholly independent measurements, as the latter incorporates positional information from the former.) The UCAC4, Tycho-2, and Gaia DR1 proper motions are measured in an absolute reference frame and are therefore not directly comparable to the local proper motions measured here. We would expect to see a roughly constant offset between these absolute proper motions and our local ones, with that offset representing the bulk motion of the Carina Nebula relative to the Sun. We plot the available absolute proper motions for each star, along with our measured local proper motions, in Figure 5. Contrary to expectations, Figure 5 does not show a consistent offset between local and absolute proper motions. The UCAC4 proper motions in particular do not follow any apparent trend relative to the local proper motions or the Tycho-2 and Gaia DR1 data. The differences between catalogues suggest that there may be systematic effects in the literature measurements that are not taken into account in the published uncertainties.

For further comparison, we correct the Gaia proper motions, where available, to the rest frame of the Carina Nebula in two ways. First, we formally correct for Galactic rotation and solar peculiar motion, as in Moffat et al. (1998, 1999) and Comerón & Pasquali (2007). We adopt Oort’s constants $A=15\pm 1$ km s*-1* kpc*-1* and $B=-12\pm 1$ km s*-1* kpc*-1* (Feast & Whitelock, 1997; Elias et al., 2006; Bovy, 2016) and components of the solar peculiar velocity $(U_{\sun},V_{\sun},W_{\sun})=(10\pm 1,12\pm 1,7\pm 1)$ km s*-1* (Feast & Whitelock, 1997; Elias et al., 2006; Schönrich et al., 2010; Tetzlaff et al., 2011). The corrected proper motions are plotted in Figure 6. For all three stars, the corrected Gaia proper motions are $\leq 1.3$ mas yr*-1* ($\leq 14$ km s*-1*), supporting our interpretation that these stars are not runaways. The corrected Gaia proper motion of ALS 15206 (ERO 2) is, like our measured motion, directed to the northeast, tangential to the orientation of its candidate bow shock. The corrected Gaia motion of TYC 8626-2506-1 (ERO 23) is to the southwest, into its candidate bow shock, although its 1$\sigma$ uncertainties overlap with those of our measured motion to the southeast. The corrected Gaia motion of HD 93576 (ERO 31) is also consistent with our data, although the Gaia results indicate a smaller velocity to the west, suggesting a longer time since ejection from Bo 11.

We also perform an empirical correction to the local reference frame: we compute the weighted mean proper motion of the 38 O-type stars in the Carina Nebula in Gaia DR1 (which is roughly half the total O-type population of the region; e.g., Smith, 2006a; Gagné et al., 2011; Alexander et al., 2016) and subtract that from the absolute Gaia proper motions of the three sample stars. The results are consistently $\sim 1.1$ mas yr*-1* ($\sim 12$ km s*-1*) offset from the results of formally correcting for Galactic rotation and solar peculiar motion. For ALS 15206 (ERO 2) and TYC 8626-2560-1 (ERO 23), the empirical correction brings the corrected Gaia proper motions into better agreement with our results. For HD 93576 (ERO 31), the empirical correction produces worse agreement with our results and suggests that the star is moving to the east, away from its candidate bow shock and toward Bo 11. The different parts of the Carina Nebula may have different large-scale motions not properly accounted for in these corrections. Future Gaia data releases, extending into Carina’s intermediate-mass population, will allow more precise and locally-specific corrections to the local reference frame.

4.3 Interpreting Bow Shocks in Giant H ii Regions

In our subsample of bow shock candidates in the Carina Nebula, EROs 2 and 23 face the OB cluster Tr 14, while EROs 24 and 31 point more generally toward Tr 14 and 16. The majority of bow shock candidates in the full Sexton et al. (2015) sample also point in toward the clusters rather than out as would be expected for runaway stars. Sexton et al. (2015) hypothesized that these candidate bow shocks are markers of large-scale ISM flows driven by cluster feedback. The ionized gas in the Carina Nebula is known to be globally expanding at 15–20 km s*-1* (Walborn & Hesser, 1975; Walborn et al., 2002b, 2007), with multiple local centers of expansion including Tr 14 (Damiani et al., 2016). Feedback-driven outflows are also thought to explain the orientations of bow shocks in other massive star-forming regions (Povich et al., 2008; Winston et al., 2012) and the correlation of bow shock orientations on small angular scales (Kobulnicky et al., 2016). Our results are broadly compatible with this interpretation, but they indicate that the factors influencing a bow shock or bow-shock-like structure cannot be deduced solely from its orientation.

The associated stars of EROs 2 and 23 are not moving in the direction of their infrared arcs, which are thus likely shaped by feedback from Tr 14. These two objects confirm that, at least in this environment, apparent bow shock orientation does not always follow stellar motion. In contrast, the associated stars of EROs 24, 25, and 31 are moving roughly toward Tr 14 and 16. Both stellar motion and ISM flows could be relevant in setting the orientation of these three candidate bow shocks, demonstrating that cluster-facing bow-shock-like structures are not necessarily clear markers of the motion of the ISM.

Feedback may also be affecting bow shocks in giant H ii regions in other ways besides large-scale outflows. For instance, the pile-up of dust in bow shock arcs depends on the presence of dust in the surrounding H ii region, as the hot winds of OB stars do not make dust effectively. This dust may originate in photoevaporative flows off nearby molecular pillars, driven globally by ionizing radiation from the central clusters and locally by individual OB stars. Kobulnicky et al. (2016) found that eight percent of bow shocks across the Galactic Plane face bright-rimmed clouds, suggesting they are shaped by local photoevaporative flows. The arc-shaped dust waves around $\sigma$ and $\lambda$ Ori are also thought to be driven by photoevaporative flows off the edge of ionized bubbles (Ochsendorf et al., 2014a; Ochsendorf et al., 2014b; Ochsendorf & Tielens, 2015; Ochsendorf et al., 2015). Density gradients in the ISM can also affect bow shock symmetries (Wilkin, 2000) or create infrared arcs via uneven heating.

We inspected Spitzer images of the Carina Nebula to assess the relationship between Carina’s EROs, its molecular gas, and the distribution of warm dust. There is a possible tendency for EROs to be closer to dense pillars than expected from a random distribution, but there is no correlation between ERO orientation and the direction to the nearest pillar. ERO 31, for example, lies just 45″ ($\sim$0.5 pc) from the edge of a prominent pillar, but points almost directly away from it. Multiband Imaging Photometer (MIPS) images at 24 $\micron$ reveal complex warm dust structures throughout the nebula, including around EROs 24 and 31. However, the origin and impact of these structures with respect to the EROs is unclear. Higher-resolution mid-infrared imaging is required to tease out the effects of density gradients and photoevaporative flows in shaping Carina’s EROs.

In any case, our main result is unaffected: in a giant H ii region, the orientation of bow-shock-like structures may be determined by the ISM, by stellar motion, or by some combination of factors. It is worth reiterating that none of the five stars in our study are runaways from Tr 14 or 16. While this result is unsurprising given the orientation of their candidate bow shocks, it confirms that stars with bow-shock-like structures are not automatically runaways. The statistical preference for alignment between stellar motion and bow shock orientation, particularly among known runaway stars (van Buren et al., 1995; Kobulnicky et al., 2016) suggests that stellar motion does dominate over ISM flows for bow shocks far from feedback-generating clusters. But within associations, assumptions about the implications of bow shocks and bow-shock-like structures (e.g., Kobulnicky et al., 2010; Gvaramadze et al., 2011a) should be made with caution.

4.4 Implications for the Origins of OB Associations

The local proper motions of CPD-59 2605 (ERO 24), HD 305533 (ERO 25), and HD 93576 (ERO 31) can also shed light on the origins of the distributed massive-star population in the Carina Nebula. Nearly half of Carina’s massive stars, including the WNH star HD 93131, are spread across roughly 20 pc in the South Pillars (Smith, 2006a). Some of these massive stars are associated with small open clusters (Bo 11, Cr 228) and other groups and subclusters of young stars (Smith et al., 2010b; Feigelson et al., 2011). However, Herschel imaging detected no massive protostars in the region, suggesting that the ongoing star formation in the South Pillars is limited to low- and intermediate-mass stars (Gaczkowski et al., 2013).

In the classic picture of clustered star formation (e.g., Lada & Lada, 2003), massive stars rarely form in a distributed mode as seen in the South Pillars. Instead, massive stars are born in clusters that may subsequently become unbound after gas dispersal and expand into OB associations (Tutukov, 1978; Hills, 1980; Lada & Lada, 1991, 2003). In this picture, one would expect the Carina Nebula’s distributed massive stars to have formed in the central Trumpler clusters and drifted out to their current locations over several Myr.

Our proper motion results are inconsistent with this expectation, as all three of the massive South Pillars stars measured here are moving toward the Trumpler clusters, not away. These stars’ kinematics suggest that they were born in the South Pillars, possibly in one of the smaller open clusters, and support a model of star formation in which OB associations form directly as loose aggregates (e.g., Efremov & Elmegreen, 1998; Clark et al., 2005). A similar result has been observed for the Cyg OB2 association based on its substructure and lack of global expansion (Wright et al., 2014, 2016). Further investigation of stellar kinematics in the South Pillars is needed to confirm this interpretation of the Carina Nebula’s distributed population.

5 Conclusions

Using HST ACS imaging with 9–10 year baselines, we have measured the local proper motions (i.e., relative to the surrounding stars) of five OB stars associated with candidate bow shocks in the Carina Nebula. Because these stars are highly saturated in our data, we use precisely-aligned images to measure the shift in each star’s Airy rings between epochs. The results are largely upper limits, but we are able to constrain the direction of each star’s motion for comparison to the orientation of its candidate bow shock.

Stellar wind bow shocks are formed when the relative velocity between star and ISM is supersonic, but the bow shock alone does not indicate which component of the relative velocity dominates. Are bow shocks indicators of fast-moving runaway stars or do they mark the large-scale flow of the ISM? In our sample of five, we find two cases where the latter is likely the case, as the stars are moving at a tangent to the arc of their candidate bow shocks. In the other three cases, we conclude that the possible influences of ISM flows, ISM structure, and stellar motion cannot be separated, and that multiple factors could be relevant for each object. We consequently caution against overinterpreting the orientation of bow shocks and bow-shock-like structures in giant H ii regions like the Carina Nebula.

In addition, none of the five stars measured here are runaways from the central OB clusters of the Carina Nebula, although two may have been ejected from smaller open clusters in the South Pillars. This finding emphasizes that bow shocks and bow-shock-like structures in giant H ii regions are not definite markers of runaway stars. It also suggests that the distributed massive-star population in the Carina’s South Pillars formed along with the distributed low- and intermediate-mass population; the resulting OB association is not the expanding remnant of an embedded cluster but a loose collection of many small groups and clusters.

Acknowledgements

The authors would like to thank Jay Anderson for providing us with his suite of PSF-fitting and image alignment software, and for his valuable instruction, guidance, and technical support. We also thank the anonymous referee for a constructive review. Support for this work was provided by NASA grants GO-13390 and GO-13791 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.