A Massive Young Runaway Star in W49 North
Luis F. Rodriguez, Roberto Galvan-Madrid, Joel Sanchez-Bermudez, Christopher G. De Pree
DDraft version February 10, 2020
Typeset using L A TEX preprint style in AASTeX63
A Massive Young Runaway Star in W49 North
Luis F. Rodr´ıguez,
1, 2
Roberto Galv´an-Madrid, Joel Sanchez-Bermudez,
3, 4 andChristopher G. De Pree Instituto de Radioastronom´ıa y Astrof´ısica, Universidad Nacional Aut´onoma de M´exico,, Apdo. Postal 3-72(Xangari), 58089 Morelia, Michoac´an, M´exico. Mesoamerican Center for Theoretical Physics, Universidad Aut´onoma de Chiapas, Carretera Emiliano Zapata Km.4,Real del Bosque (Ter´an). 29050 Tuxtla Guti´errez, Chiapas, M´exico Instituto de Astronom´ıa, Universidad Nacional Aut´onoma de M´exico, Apdo. Postal 70264, Ciudad de M´exico 04510,M´exico. Max-Planck-Institut fr Astronomie, Knigstuhl 17, 69117 Heidelberg, Germany. Department of Physics & Astronomy, Agnes Scott College, 141 East College Avenue, Decatur, GA 30030, USA.
ABSTRACTWe analyzed high angular resolution 45.5 GHz images of the W49 North massivestar forming region obtained in 1998 and 2016 with the Very Large Array. Most ofthe ultracompact HII regions show no detectable changes over the time interval of theobservations. However, subcomponents B1, B2, G2a and G2c have increased its peakflux densities by values in the range of 3.8 to 21.4 %. Most interestingly, the cometaryregion C clearly shows proper motions that at the distance of the region are equivalentto a velocity of 76 ± − in the plane of the sky. We interpret this region as theionized bowshock produced by a runaway O6 ZAMS star that was ejected from theeastern edge of Welch’s ring about 6,400 years ago. Keywords: stars: formation – stars: massive – astrometry – ISM: individual (W49North) INTRODUCTIONDynamical interactions in young stellar systems have important consequences to cluster- and stellarevolution (Zinnecker & Yorke 2007; Reipurth 2000; Tan et al. 2014). An important signature of suchinteractions is the existence of a population of runaway stars (Poveda et al. 1967). A significantfraction of the field O-type stars within a few kpc of the Sun could be runaways (Ma´ız-Apell´aniz et al.2004). A few massive stars have been reported as runaways in optically-visible clusters such as R136 inthe Large Magellanic Cloud (Lennon et al. 2018) and Cygnus OB2 in our Galaxy (Comer´on & Pasquali2007), inferring travel times (cid:38) ∼
500 yr. Refined measurements [email protected] a r X i v : . [ a s t r o - ph . S R ] F e b confirmed the original findings (Rodr´ıguez et al. 2017) and unveiled additional runaway candidateswith similar ejection times (Dzib et al. 2017; Luhman et al. 2017). Other young runaway stars inthe vicinity of the Orion Nebula Cluster with longer ejection times have been reported by McBride& Kounkel (2019) and Dzib et al. (2017).In this Letter, we report on the existence of a massive runaway moving away from the Welch ringof Ultracompact (UC) HII regions (Welch et al. 1987; De Pree et al. 1997) in W49N, one of the mostmassive, concentrated, and luminous cluster forming regions in the Galaxy (Galv´an-Madrid et al.2013; Lin et al. 2016). We also report on the flux variability of the 7 mm continuum in several of theUC HIIs in the field, and discuss them in the context of the variability study at 3.5 cm presented inDe Pree et al. (2018). MULTIEPOCH DATAWe analyzed high-angular resolution Q-band VLA data. The older image comes from concatenatingdata of three projects (AD356, 1995 May 30, D-configuration; AD414, 1998 Apr 11, A-configuration;and TEST, 2001 Apr 09, B-configuration). Since most of the data comes from observations in the Aconfiguration we attribute to these concatenated data an epoch of 1998.28. The more recent Q bandimage was taken as part of Project 16B-022 on 2016 Sep 28 (epoch = 2016.74) with the VLA in theA configuration. The gain calibrator for 1998 Apr 11 and 2016 Sep 28 was the same, J1922+1530.However, the positions for the gain calibrator between these two epochs differed by ∼ (cid:48)(cid:48)
05. As wewill see below, the method used to align the two epochs will correct for these small differences, aswell as others coming from tropospheric phase noise and self-calibration.Both images were taken at a central frequency of 45.5 GHz, with a total bandwidth of 100 MHzfor the 1998.28 data and of 2048 MHz for the 2016.74 data. The 1998.28 data were calibrated usingthe standard AIPS (Astronomical Image Processing System) techniques, while the 2016.74 data werecalibrated using the CASA (Common Astronomy Software Applications) package of NRAO and thepipeline provided for VLA observations. The data were self-calibrated and the resulting imageshave angular resolutions of 0. (cid:48)(cid:48) × (cid:48)(cid:48) ◦ (cid:48)(cid:48) × (cid:48)(cid:48) ◦
5, for 1998.28 and2016.74, respectively. Both images were smoothed to an angular resolution of 0. (cid:48)(cid:48) × (cid:48)(cid:48) ∼ (cid:48)(cid:48) using the AIPS task MWFLT and following the prescription of Rudnick(2002). IMAGE ANALYSISComparison of the two images showed an offset of the order of 0. (cid:48)(cid:48) ∼ − (Zhang et al. 2013) and over the 18 year interval between the images it will accumulate to adisplacement of order 0. (cid:48)(cid:48)
01. Most likely the displacement between images is due to a combination ofthe tropospheric conditions during the observations, the small position displacements that can resultfrom self-calibration and the difference in position for the gain calibrator in the two epochs. To alignthe two images as best as posible we followed two procedures. The positions of the most compact andbright sources (components A, B2 and G2a in the nomenclature of Dreher et al. (1984) and De Preeet al. (1997)) were compared, finding that the 1998 image was displaced by ∆ α = − (cid:48)(cid:48) ± (cid:48)(cid:48) https://science.nrao.edu/facilities/vla/data-processing/pipeline ∆ δ = − (cid:48)(cid:48) ± (cid:48)(cid:48) (cid:48)(cid:48)
01 both in α and δ and searching for theresidual image with the smallest rms. A displacement of ∆ α = − (cid:48)(cid:48)
02; ∆ δ = − (cid:48)(cid:48)
11 was found to beoptimal and we adopted it to compare the two images.Finally, the absolute flux density calibration at Q band is known to present random variations ofthe order of 10%. Again we searched for the image subtraction that produced the smallest rms andthis was obtained multiplying the 1998.28 image by a correction factor of 1.09. The final differenceimage shows significant residuals only in association with five components: B1, B2, C, G2a and G2c.Bright sources such as components A1 and A2 cancel out in the subtraction of the two images (seeFigure 1), indicating that they have not experienced significant changes in position, morphology oremission. FLUX VARIABILITYThe images associated with components C, B, and G are shown in Figure 2. In the case of theimages of components B and G we find that the subcomponents B1, B2, G2a and G2c appear tohave increased their peak intensity by 3.8 ± ± ± ± (cid:48)(cid:48) (cid:48)(cid:48) n e , the small HII region could appear fainter at 3.6cm and brighter at 7 mm. Altogether, the results of De Pree et al. (2018) and ours are consistent withthe predictions of Galv´an-Madrid et al. (2011), who used the radiation-hydrodynamical simulationsof Peters et al. (2010) to estimate a ∼
15% probability for an increase- and a ∼
6% probability for adecrease in the 2-cm flux of UC/HC HII regions over a time span of 20 years. The number of sourceswith detected time variability is also roughly consistent with the estimations of Galv´an-Madrid et al.(2011), who concluded that ∼
10% of the sources should have detectable variations within a periodof ∼
10 years. RADIO SOURCE C AS A RUNAWAY FROM W49NIn contrast to the components discussed above, component C clearly shows proper motions, signaledby the shifted negative and positive contours in Fig. 3, present at a 10- σ significance level. Toestimate the angular size of the displacement, again we took the 1998.28 image and displaced it byincrements of 0. (cid:48)(cid:48)
001 both in α and δ . The initial image was subtracted to the displaced image until theresiduals were minimized. The best agreement, shown in Fig. 3, was obtained with a displacementof ∆ α = − (cid:48)(cid:48) ± (cid:48)(cid:48) δ = +0. (cid:48)(cid:48) ± (cid:48)(cid:48) s = 0. (cid:48)(cid:48) ± (cid:48)(cid:48)
002 ata position angle of 309 ◦ ± ◦ . This displacement is equivalent to a proper motion of 1.41 ± − .The direction of the proper motion of source C matches the axis of its arc morphology (see Figs.2 and 4), suggesting that the emission arises from a bowshock (e.g., Arthur & Hoare 2006). At adistance of 11.1 kpc (Zhang et al. 2013), the proper motions of component C over the time intervalof 18.46 yr are equivalent to a velocity in the plane of the sky of 76 ± − . The position angleof the proper motions suggests that the star was ejected from the eastern edge of Welch’s ring (seeFigure 4). The present angular distance between the component C and the eastern edge of Welch’sring is ∼ (cid:48)(cid:48) , implying that this star was ejected about 6,400 years ago. The other sources in theregion do not show significant evidence for proper motions, with an upper limit to their plane-of-thesky velocities of ∼
20 km s − .To further test the hypothesis of a bow shock from a runaway star, we modelled the emission profileusing the algebraic solution of Canto et al. (1996). First, a radial profile of the emission was extractedat 120 different position angles. The position of the emission centroid and its uncertainty at each P.A.was obtained through Gaussian fitting. The derived bowshock profile consists of the locus containingall the fitted emission centroids. This profile is then modeled according to the following equation: R ( θ ) = R csc θ (cid:112) − θ cot θ ) , (1)where R is the stand-off distance between the apex of the bowshock and the position of the starand θ the different opening angles of the bow-shock profile. R depends on the pressure-momentabalance between the relative motion of the stellar wind and the circumstellar medium as follows: R = (cid:115) ˙ mv w πρv ∗ , (2)where ˙ m is the mass-loss rate of the star, v w is the terminal velocity of the wind, ρ is the mean densityof the circumstellar medium, and v ∗ is the relative velocity of the star through the environment.To determine R we used a Monte-Carlo Markov-Chain model fitting using the emcee code(Foreman-Mackey et al. 2013), leaving ˙ m and the projected position angle of the bow shock inthe plane of the sky as free parameters. v w was set to 2000 km s − , a typical parameter for massivestars (e.g., Kudritzki & Puls 2000). ρ = µm H n H was obtained from previous mass estimations ofthe W49N molecular clump which give an average density n H ∼ cm − (Galv´an-Madrid et al.2013), and using µ = 2 . v ∗ is set to 76 km s − asdetermined from the observed proper motions here reported. It should be noted that most massivestars are part of binary systems and that the velocity adopted here is actually the velocity of thebowshock. The true velocity of the star with respect to the medium could be modulated by motionsaround a companion.Figure 5 displays the histogram of the best-fit values for the mass-loss rate. We obtain a mass-lossrate ˙ m =2.62 × − ± × − M (cid:12) /yr, which corresponds to a stand-off distance of R =2216 ± ◦ (E → N). Figure 5 also displays the bowshock with the best-fit profile over-plottedon the radio continuum map, as well as the derived position of the central star.We use the model-derived ˙ m to estimate that the runaway has a stellar mass M (cid:63) ≈ ± M (cid:12) (Dzibet al. 2013; Vacca et al. 1996), corresponding to a main-sequence spectral type of O6 ± K, we derive that an ionizing photon rate of 5 . × s − is requiredto keep the photoionization (Snell et al. 2019). An O8 ZAMS star can provide this photoionization(Dzib et al. 2013). An explanation for the apparent discrepancy between the spectral types derivedfrom the mass loss rate and the ionizing photon rate is that most of the ionizing photons of thecentral star escape to its SE, where there is no gas to stop the photons. If we assume that only 1/3 ofthe photons from the central star are absorbed by component C, an O6 ZAMS star is then requiredand excellent agreement is obtained between the two independent estimates. The assumption thatonly 1/3 of the ionizing photons are intercepted by component C is consistent with the geometryshown in the right panel of Fig. 5. In this Figure the opening angle of component C as seen from thecentral star θ (cid:39) ◦ . The fraction of the 4 π solid angle around the central star that is interceptedby component C is (1 / − cos ( θ/ (cid:39) .
25, consistent with the required value of 0.3.De Pree et al. (2004) report the local standard of rest (LSR) H52 α radial velocity of 15 ultracompactHII regions in the zone. Excluding component C, the mean and rms values of these velocities are11.9 ± − . Component C has v LSR = -8.0 ± − , blueshifted by more than 3 σ awayfrom the mean value. This confirms that component C has different motions than the rest of thestellar sources and the dense-gas clump and provides strong supporting evidence of the anomalouskinematics of component C.The potential well of the W49N region is given by the sum of the molecular gas mass within a 3pc radius M gas ≈ × M (cid:12) and the more uncertain embedded stellar mass M (cid:63) ∼ M (cid:12) (Galv´an-Madrid et al. 2013). Thus, the escape velocity from the embedded cluster is v esc ≈ . − . The3D velocity of source C with respect to the stellar cluster is v C , ≈ . − . Therefore, sourceC is clearly able to escape from the W49N clump region and constitutes a massive, runaway youngstar. If its velocity remains unchanged, source C could escape from the W49N clump within the next4 to 5 × yr, which is an order of magnitude smaller than the (cid:46) > ∼ SUMMARYWe analyze high-angular resolution 7-mm VLA data of the Welch ring of UC HIIs in W49N obtainedin epochs 1998.28 and 2016.74. We find that source C has large proper motions of 1 . ± .
11 masyr − toward the NW (P.A.= − ◦ ± ◦ ), equivalent to a velocity in the plane of the sky of 76 ± − . Together with its line-of-sight velocity ≈ −
20 km s − blueshifted from the mean cluster velocityas previously measured from recombination lines, we infer that source C is a young, massive runawaystar, which was ejected from the Welch ring ∼ D ec li n a t i on ( J2000 )
09 06 12.412.212.011.811.611.411.211.0 D ec li n a t i on ( J2000 )
09 06 12.412.212.011.811.611.411.211.0 D ec li n a t i on ( J2000 ) Right Ascension (J2000)19 10 12.94 12.92 12.90 12.88 12.86 12.8409 06 12.412.212.011.811.611.411.211.0
Figure 1.
Images of the components A1 and A2 for 1998.28 (top), 2016.74 (middle) and difference image(bottom). The beam (0. (cid:48)(cid:48) × (cid:48)(cid:48)
1) is shown in the bottom left corner. The contours are -4, 4, 5, 6, 8, 10, 12,15, 20, 30, and 40 times 1.6, 1.6 and 1.8 mJy beam − , respectively.. The residuals in the difference imageare not considered significant. D ec li n a t i on ( J2000 )
09 06 19.219.018.818.618.418.218.0 D ec li n a t i on ( J2000 )
09 06 13.213.012.812.612.412.212.0 D ec li n a t i on ( J2000 )
09 06 13.613.413.213.012.812.612.412.212.0 D ec li n a t i on ( J2000 )
09 06 19.219.018.818.618.418.218.0 D ec li n a t i on ( J2000 )
09 06 13.213.012.812.612.412.212.0 D ec li n a t i on ( J2000 )
09 06 13.613.413.213.012.812.612.412.212.0 D ec li n a t i on ( J2000 ) Right Ascension (J2000)19 10 13.20 13.19 13.18 13.17 13.16 13.15 13.14 13.13 13.1209 06 19.219.018.818.618.418.218.0 D ec li n a t i on ( J2000 ) Right Ascension (J2000)19 10 13.18 13.17 13.16 13.15 13.14 13.13 13.12 13.11 13.1009 06 13.213.012.812.612.412.212.0 D ec li n a t i on ( J2000 ) Right Ascension (J2000)19 10 13.48 13.46 13.44 13.42 13.40 13.38 13.36 13.3409 06 13.613.413.213.012.812.612.412.212.0
Figure 2.
Left column: images of the component C for 1998.28 (top), 2016.74 (middle) and difference image(bottom). The contours are -10, -8, -6, -5, -4, 4, 5, 6, 8, 10, 12, 15, 20, 30, and 40 times 1.6, 1.6 and 1.4 mJybeam − , respectively. The residuals in the difference image are interpreted to imply a motion of the source. Center column: images of the components B1, B2 and B3 1998.28 (top), 2016.74 (middle) and differenceimage (bottom). The contours are -5 -4, 4, 5, 6, 8, 10, 12, 15, 20, 30, 40 times and 60 times 1.6, 1.6 and 1.8mJy beam − , respectively. The residuals in the difference image are interpreted to imply an increase in thepeak flux density of components B1 and B2. Right column: images of the component G1, G2a, G2b andG2c for 1998.28 (top), 2016.74 (middle) and difference image (bottom). The contours are -4, 4, 5, 6, 8, 10,12, 15, 20, 30, 40, 60 and 80 times 1.6, 1.6 and 2.2 mJy beam − , respectively. The residuals in the differenceimage are interpreted to imply an increase in the peak flux density of components G2a and G2c. The beam(0. (cid:48)(cid:48) × (cid:48)(cid:48)
1) is shown in the bottom left corner of each panel. D ec li n a t i on ( J2000 )
09 06 19.219.119.018.918.818.718.618.518.418.318.2
DATA D ec li n a t i on ( J2000 )
09 06 19.219.119.018.918.818.718.618.518.418.318.2
MODEL D ec li n a t i on ( J2000 ) Right Ascension (J2000)19 10 13.19 13.18 13.17 13.16 13.15 13.14 13.13 13.1209 06 19.219.119.018.918.818.718.618.518.418.318.2
MODEL - DATA
Figure 3.
Image difference for the component C (top), model for the component C made as described inthe text (middle) and difference of the model minus data (bottom). The beam (0. (cid:48)(cid:48) × (cid:48)(cid:48)
1) is shown in thebottom left corner. The contours are -10, -8, -6, -5, -4, 4, 5, 6, 8, and 10 times 1.4 mJy beam − . Figure 4.
Image of W49 North for epoch 1998.28. The beam (0. (cid:48)(cid:48) × (cid:48)(cid:48)
1) is shown in the bottom leftcorner. The contours are 4, 6, 10, 20, 40, 60, 80 and 100 times 1.6 mJy beam − . The arrow indicates thedisplacement of source C for a time span of 10 years. Figure 5.
Left:
Histogram of the best-fit mass-loss rate. The vertical red-dotted lines indicate the 16%,50% and 84% of the percentiles in the distribution. The mass loss rate is given in units of M (cid:12) yr − . Right:
Image of the bow shock with the extracted profile (red-line); the best-fit model of the shape (white-line);the position of the central source (white star) and; the position of the bow-shock apex (dotted-white arrow).The angular scale of the image and its orientation are displayed in the panel.
ACKNOWLEDGMENTSWe thank an anonymous referee for suggestions that improved the clarity of the paper. LFRacknowledges the financial support of PAPIIT-UNAM and of CONACyT (M´exico). RGM acknowl-edges support from UNAM-PAPIIT project IN104319. CGD acknowledges support from NSF RUIaward 1615311. This research has made use of the SIMBAD database, operated at CDS, Strasbourg,France.
Facility:
VLA
Software:
AIPS (van Moorsel et al. 1996); CASA (McMullin et al. 2007); emcee (Foreman-Mackeyet al. 2013). REFERENCES
Arthur, S. J., & Hoare, M. G. 2006, ApJS, 165,283Canto, J., Raga, A. C., & Wilkin, F. P. 1996,ApJ, 469, 729Comer´on, F., & Pasquali, A. 2007, A&A, 467, L23 De Pree, C. G., Mehringer, D. M., & Goss, W. M.1997, ApJ, 482, 307De Pree, C. G., Wilner, D. J., Mercer, A. J., et al.2004, ApJ, 600, 286De Pree, C. G., Galv´an-Madrid, R., Goss, W. M.,et al. 2018, ApJL, 863, L91