NGC 7538 IRS1—an O Star Driving an Ionized Jet and Giant N–S Outflow

Göran Sandell; Melvyn Wright; Rolf Güsten; Helmut Wiesemeyer; Nicolas Reyes; Bhaswati Mookerjea; Stuartt Corder

doi:10.3847/1538-4357/abbf5b

1. Introduction

The molecular cloud south of the H ii region NGC 7538 harbors several young massive stars, of which at least three, IRS 1, IRS 9, and NGC 7538 S, are centers of young clusters. At a distance of 2.65 kpc (Moscadelli et al. 2008) NGC 7538 is by no means the closest high-mass star-forming region, but the extreme youth of these stars makes them key targets for our understanding of high-mass star formation. NGC 7538 IRS 1 was first detected as an OH maser by Hardebeck (1971), indicating that there is another H ii region SE of the optically visible H ii region NGC 7538. It was first detected in the radio at 5 GHz by Martin (1973), who found three compact radio sources at the SE edge of the compact (∼4') H ii region NGC 7538. The brightest of the three, i.e., Martin's source B, is commonly referred to as NGC 7538 IRS 1, based on near/mid-IR high spatial resolution mapping and photometry by Wynn-Williams et al. (1974). In addition to the radio sources IRS 1–IRS 3, Qiu et al. (2011) also found eight cores at 1.3 mm using the Submillimeter Array (SMA) within a projected distance of 0.35 pc from IRS 1, many of which are associated with H₂O maser emission or outflow activity, suggesting that they are young pre-main-sequence stars. Frau et al. (2014) increased the number of cores to 13 by more sensitive observations with the SMA at 880 μm. IRS 1 is the most massive star in a cluster with more than 150 young stellar objects (YSOs) and protostars (Mallick et al. 2014; Sharma et al. 2017).

The radio emission from IRS 1 was first resolved with the Karl G. Jansky Very Large Array (VLA) at 14.9 GHz by Campbell (1984), who showed that the radio emission consists of a compact (∼02) bipolar N–S core with faint extended fan-shaped lobes, suggesting an ionized outflow. High angular resolution observations of radio recombination lines (Gaume et al. 1995; Sewilo et al. 2004; Keto et al. 2008) revealed extremely broad lines, indicating substantial mass motion of the ionized gas. Keto et al. (2008) modeled the spectral energy distribution (SED) of IRS 1 as a combination of two H ii regions plus dust, which could adequately fit the data available at the time. However, Sandell et al. (2009), who analyzed high angular resolution VLA archive data from 4.86 to 43.4 GHz, found that the bipolar core had a spectral index of ∼0.8 and that the size of the core shrank as a function of frequency with a power-law index of ∼−0.8, which is what one would expect from an ionized jet (e.g., Reynolds 1986), but not from an optically thick H ii region. It is therefore clear that IRS 1 drives a bipolar ionized jet.

IRS 1 is heavily accreting with an accretion rate of ∼10⁻³ to 2 × 10⁻⁴ ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Sandell et al. 2009; Klaassen et al. 2011; Qiu et al. 2011; Beuther et al. 2012; Zhu et al. 2013), which quenches the expansion of an H ii region. IRS 1 itself is an extremely rich maser source, with not only the normal masers found in high-mass stars, like OH, H₂O, and CH₃OH, but also rare masers like NH₃, H₂CO, and the 23.1 GHz class II CH₃OH maser. The latter has to date been only detected in three sources; see Galván-Madrid et al. (2010) and references therein. These masers are most likely pumped by the radio emission from the jet. IRS 1 is almost certainly surrounded by an accretion disk, but so far there is no firm evidence for a Keplerian accretion disk surrounding IRS 1. The orientation of the linear polarization of the 157 GHz CH3OH class II maser toward IRS 1 agrees with that of the submillimeter dust polarization, suggesting a magnetic field oriented along the outflow axis (Wiesemeyer et al. 2004). The same holds for the 133 GHz class I maser, residing 40'' south of IRS 1 and possibly associated with a bow shock in the redshifted outflow lobe. While it is unclear whether the 157 GHz maser originates from an accretion disk or from the outflow, there have been many claims that IRS 1 is surrounded by an accretion disk in Keplerian rotation (Minier et al. 1998, 2000; Pestalozzi et al. 2004, 2009; Surcis et al. 2011; Moscadelli & Goddi 2014) based on linear methanol maser features. Moscadelli & Goddi (2014) identified three maser groups and suggested that IRS 1 is a triple system, where each high-mass YSO in the system is surrounded by an accretion disk and probably drives an outflow. Beuther et al. (2017), who imaged IRS 1 at 23 GHz with the VLA in the A configuration, identified two of these maser features being close to peaks in the free–free emission and argued that IRS 1 is a binary, powered by two hypercompact H ii regions, each with a spectral index of 1.5–2. This does not agree with the observed spectral index of IRS 1 and with high angular (CARMA⁶ ) A-array data analyzed in this paper. Here we will therefore reexamine whether IRS 1 is a single star or a multiple system.

Scoville et al. (1986), using the Owens Valley Radio Observatory (OVRO) millimeter array at 2.7 mm, was the first to positively identify the high-velocity CO (1–0) emission seen toward NGC 7538 as a bipolar outflow powered by IRS 1 with the outflow at a position angle (PA) of 145° (counterclockwise against north), being redshifted in the SE and blueshifted to the NW with a total extent of ∼3', i.e., 2.3 pc. This agrees well with other studies (Fischer et al. 1985; Kameya et al. 1989; Davis et al. 1998). However, if this is the case, the ionized jet and the molecular outflow are misaligned. This led Kraus et al. (2006) to interpret mid-IR speckle images as a precessing jet where IRS 1 would undergo non-coplanar tidal interaction with an undiscovered close companion within the circumbinary protostellar disk. Sandell et al. (2012) presented single-dish James Clerk Maxwell Telescope (JCMT) CO (3–2) data, which show blueshifted high-velocity gas ∼80'' north of IRS 1, which suggested that the molecular outflow is quite large and well aligned with the ionized outflow. Here we present an extended CO (3–2) map from JCMT, as well as [C ii] 158 μm and high-J CO data obtained with GREAT on the Stratospheric Observatory for Infrared Astronomy (SOFIA). Analysis of these new data sets conclusively shows that IRS 1 drives a large N–S molecular and ionized outflow.

2. Observations

2.1. GREAT Observations

The NGC 7538 IRS 1 outflow was observed in the GREAT⁷ L1/L2 configuration (Heyminck et al. 2012) on SOFIA (Young et al. 2012) on 2015 January 23 as part of the open time project 02_0038. These observations were done on a 44-minute leg at 41,300 feet. The L2 mixer was tuned to the [C ii] ${}^{2}{{\rm{P}}}_{3/2}{\to }^{2}{{\rm{P}}}_{1/2}$ transition at 1.9005369 THz, while the L1 mixer was tuned to CO (11–10) at a rest frequency of 1.26701449 THz. The main-beam coupling efficiency η_mb = 0.69 for L2 and 0.65 for L1. The half-power beam width (HPBW) for L2 was measured to be 141 for [C ii] and 211 for CO (11–10). These observations were complemented with observations on 2015 May 20 and 21 during the Low Frequency Array (LFA) commissioning. The LFA array is a hexagonal array with two polarizations co-aligned around a central pixel, providing a 2 × 7 pixel array with the pixels separated by two beam widths. For a more complete description of the instrument, see Risacher et al. (2016). The LFA receiver was tuned to [C ii], while the second channel, L1, was tuned to CO (11–10). The measured HPBW for LFA was 151, and the main-beam coupling efficiency η_mb = 0.69, while the HPBW for L1 was 211 and η_mb = 0.64. On May 20, we had a 50-minute leg at 45,000 feet. On May 21 the leg duration was 52 minutes at 43,000 feet. The precipitable water vapor was ≤9 μm on both nights. In 2015 January (L1/L2) we did an OTF TP map with a size of 825 × 975 scanning in R.A. with 75 sampling and an integration time of 2 s per point. The off position was at +600'', 0'' relative to IRS 1. We also took a 10-minute pointed integration on IRS 1. In 2015 May we did much larger maps with upGREAT. These were done in classic on-the-fly (OTF) total power (TP) mode, with scanning in R.A. with a 7'' sampling and an integration time of 2 s per point using the same off position. The final map size in [C ii] is ∼160'' × 240''. For CO (11–1) the map size is 105'' × 196''. All the data were added together in the postprocessing.

The back ends for both GREAT channels are the latest generation of fast Fourier transform spectrometers (FFTS; Klein et al. 2012). These have 4 GHz bandwidth and 16,384 channels, which provide a channel separation of 244.1 kHz (0.0385 km s⁻¹ for [C ii]). The data were reduced and calibrated by the GREAT team. The post-processing was done using CLASS.⁸ We removed linear baselines, threw away a few damaged spectra, and created maps with rms weighting.

2.2. CARMA Observations

Observations at 1.3 mm were obtained in the A configuration on 2011 December 6 using 16 500 MHz bandwidth spectral windows and on 2011 December 7 using 14 500 MHz windows and two 125 MHz windows that covered the H30α recombination line with ∼1 km s⁻¹ spectral resolution. The data were calibrated using the calibration sources J0102+584, 3C 454.3, and MWC 349 in a standard way using the Miriad software package (Sault et al. 1995). Due to hardware and software limitations in the second and third local oscillator phase tracking for the wide-band correlator during these commissioning observations in the A configuration, the phase tracking center was offset in each spectral window. Apart from the offset, the images in each spectral window agreed within the noise, and we used a shift-and-add procedure (using the Miriad tasks IMDIFF and AVMATHS) to average the 500 MHz spectral windows. The 2011 December 7 data were self-calibrated using the strong H30α emission. The 2011 December 6 data used the continuum emission in the same spectral window in a 500 MHz bandwidth. The continuum images at 227.4 GHz obtained from the two data sets are consistent.

CARMA A-array observations were also obtained at 111.1 GHz on 2011 December 10 and 12, and in the B-array configuration on 2012 January 3–4 and 5–6. B-array observations were also obtained at 222.2 GHz on 2012 January 2 and 3. All these observations used the same back-end configuration as described above. These data sets were reduced in Miriad following the same procedure described in Zhu et al. (2013).

We also retrieved some mosaicked CARMA ¹²CO (1–0) and ¹³CO (1–0) images, which allows us to get a better look of the outflow close to IRS 1. The observations and data reduction of these data sets are described in Corder (2009). In particular, we use a ¹²CO (1–0) image with 45 angular resolution and a velocity resolution of 1.27 km s⁻¹. This image has an rms sensitivity of 0.64 K per channel. The total velocity coverage of the ¹²CO is from −96 to 16 km s⁻¹. We also retrieved a ¹³CO (1–0) image made from the C- and D-array configurations, resulting in a resolution of 79 × 74 at P.A. = 55°. We corrected this image for missing zero spacing by adding the ¹³CO (1–0) OSO map to this image using the MIRIAD task IMMERGE. The velocity resolution of this map is 0.33 km s⁻¹, and it covers the velocity range from −65.8 to 45.8 km s⁻¹.

2.3. CO (1–0) and ¹³CO (1–0) Mapping with the 20 m Onsala Space Telescope

We observed the NGC 7538 molecular cloud in CO (1–0) and ¹³CO (1–0) with an SIS mixer on the 20 m Onsala Space telescope in Sweden from 2010 January 28 to February 4. The back end was a 1600-channel hybrid digital autocorrelator spectrometer with 50 kHz resolution. The beam size is 33'' and 34'' for ¹²CO and ¹³CO, respectively. The main-beam efficiency is elevation dependent, but for the bulk of the data it was ∼0.42 for CO and 0.45 for ¹³CO. All observations were done in Total Power position switch mode. We mapped the cloud on a regular grid with a 15' spacing using an integration time of 3 minutes per position. In total we mapped an area of 240'' × 300'' in CO (1–0) and 400'' × 380'' in ¹³CO (1–0). Both maps were centered on IRS 1. The weather was generally quite good, although some time was lost owing to heavy snowfall. The system temperature for the ¹³CO setting was typically between 250 and 350 K, going up to about twice that much for ¹²CO. All the data reduction was done in CLASS. We removed low-order baselines and occasionally standing waves. The spectra were then co-added with rms weighting. The final maps were exported as FITS cubes and imported into MIRIAD for further analysis.

2.4. JCMT Archive Data

We obtained a fully sampled CO (3–2) map from C. Fallscheer covering a 1° × 1° region of NGC 7538. This map was observed with the Heterodyne Array Receiver Program (HARP; Buckle et al. 2009) on the JCMT⁹ at Maunakea, Hawaii. HARP has a beam size of 146 at 345.8 GHz, and the main-beam efficiency is 0.63. The reduction of the CO (3–2) data is described in Fallscheer et al. (2013). The CO data cube covers a 1 GHz wide band with 2048 channels, providing a velocity resolution of 0.42 km s⁻¹ per channel. The 1σ rms per channel is ∼0.8 K on the T_mb scale. Here we only use part of the map covering the H ii region and the molecular cloud southwest of it.

We also got two long-integration HARP maps done in position-switched jiggle mode from Jan Wouterloot covering ∼100'' × 100'', one of ¹²CO (3–2) and one of ¹³CO (3–2). Both maps have a bandwidth of 1000 MHz, i.e., a frequency resolution of 488.3 kHz, which for ¹²CO corresponds to a velocity resolution of 0.42 km s⁻¹. The ¹²CO map is a co-add of 22 observations, which were all of good quality, while the ¹³CO is based on six independent observations. The final averages were regridded onto a 30'' grid. The rms noise per channel is 30 mK for CO (3–2) and 40 mK for ¹³CO (3–2).

2.5. FORCAST Archive Data

A single field around IRS 1 was observed with the Faint Object infraRed Camera (FORCAST) on SOFIA on two flights in 2013 September 10 and 13. These observations were part of the cycle 1 project 01_0034 (A. Tielens). We retrieved the fully calibrated level 3 data from the SOFIA archive. The level 3 data products are pipeline reduced and calibrated data, which require no further processing; see Herter et al. (2013) for details about data reduction and calibration. FORCAST is a dual-channel mid-infrared camera covering 5–40 μm with a suite of broad- and narrowband filters (Herter et al. 2012). Each channel consists of 256 × 256 pixels. The pixel size after pipeline processing is 0768. Both channels can be operated simultaneously by use of a dichroic mirror internal to FORCAST. FORCAST is diffraction limited at wavelengths longer than 15 μm. At 7.7 μm the beam size is 30; at 19.7 and 37.1 μm it is 29 and 36, respectively. The observations on September 10 were done on a 1 hr leg at 38,850 feet in dual-channel mode with filter combinations 19.7 μm/31.5 μm and 25.3 μm/37.1 μm in C2NC2 mode. Each filter combination was observed for eight cycles with 30 s on-time each cycle. The observations on September 13 were done on a half-hour leg at 39,000 feet. These observations were done in single-channel mode with the 7.7 μm filter. Most of the leg was lost owing to instrument problems, resulting in only three integration cycles of 30 s on-time each cycle. Due to field rotation and boresight acquisition, some images partially capture IRS 4 at the western edge of the image, while some see IRS 11 and South at the southern border of the images. These sources are cut out in the co-added images, which only adds data common to each cycle. Even though IRS 1 and IRS 2 were completely saturated in the Spitzer IRAC images, the FORCAST images reveal no additional sources close to IRS 1 and IRS 2. A three-color image in Figure 1 shows that IRS 1–IRS 3 completely dominate the emission and are surrounded by much fainter diffuse extended emission. The hot photodissociation region (PDR) between the NGC 7538 H ii extends to the west and curves up to IRS 4, only partially seen in this image. Aperture photometry is presented in Table 1. Because IRS 1 overlaps with IRS 2, we used a small aperture and empirical aperture corrections. IRS 2 has a size of ∼10'', while IRS 1 appears marginally extended. All other sources are unresolved. The errors for IRS 4 may be underestimated, because it is at the edge of the field, and the same is true for IRS 2, which is extended and surrounded by bright emission from the surrounding cloud.

**Figure 1.** Three-color image of the IRS 1–IRS 3 region generated by combining the FORCAST 19 (blue), 25 (green), and 37 μm (red) images in logarithmic scale. IRS 1 is the strong yellow source at offset 0'', 0''. IRS 2, at 07, +96, is the extended bluish source with a size of ∼10''. The third strong source, IRS 3 is at −13'', +26. IRS 4 is at the northwestern part of the image and not fully covered in this image. IRS 1 E is the faint source east of IRS 1, and IRS 1 SE is not really visible in this color image.
Download figure:
Standard image High-resolution image

Table 1. Positions and Flux Densities of FORCAST Mid-Infrared Sources

Source	α(2000.0)	δ(2000.0)	S(7.7 μm)	S(19.7 μm)	S(25.3 μm)	S(31.5 μm)	S(37.0 μm)
	(^h ^m ^s)	(° ' '')	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)
IRS 1	23 13 45.37	+61 28 10.5	170.7 ± 3.3	240 ± 8.1	787 ± 30	1320 ± 28.2	1710 ± 50
IRS 2	23 13 45.44	+61 28 20.1	9.4 ± 5.0	508 ± 16.2	360 ± 70	550 ± 30.0	540 ± 90
IRS 3	23 13 43.62	+61 28.14.2	2.9 ± 0.2	85 ± 6.0	166 ± 17	223 ± 12.5	249 ± 25
IRS 4	23 13 32.39	+61 29 06.2	5.2 ± 0.3	⋯	⋯	41.8 ± 0.6	32 ± 0.8
IRS 1E	23 13 48.64	+61 28 05.1	0.46 ± 0.10	1.49 ± 0.20	5.4 ± 0.6	5.3 ± 0.5	6.6 ± 2.1
IRS 1SE	23 13 48.40	+61 27 40.3	<0.5	1.11 ± 0.10	3.3 ± 0.5	3.6 ± 0.5	1.2 ± 0.8

Note. Errors are 1σ statistical errors.

Download table as: ASCII Typeset image

2.6. Herschel Archive Data

We have retrieved Herschel PACS and SPIRE data from the Herschel data archive. There were two different observing projects, which covered NGC 7538 in PacsSpire parallel mode. F. Motte's HOBYS Key Programme (OBS ID 1342188088 and 1342188089) observed NGC 7538 on 2009 December 14 with slow scan speed (20'' s^–1) and the PACS blue channel set to 70 μm, i.e., PACS observed 70 and 160 μm, while SPIRE always covers all three bands, PSW (250 μm), PMW (350 μm), and PLW (500 μm), simultaneously. Preliminary analysis of these data has been presented by Fallscheer et al. (2013). The second project, S. Molinari's OT3 project (OBS ID 1342249088 and 1342249089), was observed on 2012 August 5. It had the same instrument setup but used fast scanning (60'' s^–1) and covered a much larger area. However, since NGC 7538İRS 1 was severely saturated in the HOBYS SPIRE 250 μm image, the region was reobserved as an OT2 project (P. André OBS ID 1342239268) on 2012 February 13. This map was done in small-map cross-scan mode with a scan rate of 30'' s^–1 for high dynamic range. Here we only make use of the PACS data from Motte's program (i.e., AORs 1342188088 and 1342188089) and André's high dynamic range SPIRE data. Photometry using a combination of point-spread function fitting and aperture photometry is presented in Table 2. We also integrated the emission over the cloud core surrounding IRS 1 using the STARLINK application GAIA. At 70 μm the core has an equivalent radius of ∼20'', while the emission is much more extended at 160 μm, where the radius is ∼25''. The size of the core is similar in the SPIRE images. The integrated PACS and SPIRE flux densities are given in Table 2.

Table 2. Far-infrared Photometry of IRS 1 (PACS and SPIRE) and Integrated Flux Densities of the IRS 1 Core Derived from Maps (see Section 2.6)

Filter	Flux Density	Integrated Flux
(μm)	(Jy)	(Jy)
70	6930 ± 600	7890 ± 790
160	3236 ± 300	5550 ± 500
250	1083 ± 110	2760 ± 140
350	552 ± 55	1236 ± 60
500	190 ± 20	407 ± 20

Download table as: ASCII Typeset image

3. The Nature of NGC 7538 IRS 1

Sandell et al. (2009) proposed a model where the free–free emission from IRS 1 is dominated by a collimated ionized wind driving an ionized north–south jet. This is the only model that can explain all the observed characteristics of IRS 1. The model satisfactorily explains the observed morphology of IRS 1, including the shrinking in size of the hypercompact H ii region at higher frequencies due to the outer parts of the jet becoming optically thin with increasing frequency. It also explains the observed spectral index with no turnover at high frequencies, as expected in a traditional H ii region. Further, the extremely broad recombination lines seen toward IRS 1 (Gaume et al. 1995; Keto et al. 2008) and the observed variations in flux density of timescales of 10 yr or less (Franco-Hernández & Rodríguez 2004) are consistent with the proposed model. Since there have been quite a few observations of IRS 1 since 2009, we recomputed the spectral index by fitting the observed total flux densities from 4.8 to 115 GHz. At these frequencies the free–free emission completely dominates and there is no contribution from dust. If there was a contribution from dust emission at 100 GHz, it would be quite noticeable at 220 GHz, since dust emission has a spectral index of 3–5. We now derive α = 0.87 ± 0.03, which is slightly steeper than Sandell et al. (2009) derived. In Figure 2 we have also included all recent flux densities from interferometric observations at millimeter wavelengths. Our analysis indicates that we only pick up a clear excess due to dust emission at frequencies higher than 300 GHz. We have also redone our fit to the size (length) of the compact core as a function of frequency including the A-array data presented in this paper. We find that the size of the core falls off as ν^{−0.92 ± 0.02}, which is consistent with a jet (Reynolds 1986). If we assume that the jet is launched from a radius of 10–20 au, this fit predicts that the free–free emission will not become completely optically thin until somewhere between 700 and 1500 GHz. After that, the SED would follow the typical −0.1 law for optically thin free–free emission.

**Figure 2.** Top: least-squares power-law fit to the length of the compact core surrounding IRS 1 as a function of frequency. The length of the core falls off as ν^{−0.92 ± 0.02}, where ν is the frequency in GHz. Bottom: least-squares power-law fit to total flux densities from published observations of IRS 1 from 4.8–115 GHz. The fit gives a spectral index of 0.87 ± 0.03. The flux densities in the 200 and 300 GHz windows (open squares) were excluded from the fit. Clear excess from dust emission is only seen at frequencies above 300 GHz.
Download figure:
Standard image High-resolution image

Figure 3 shows the 43 GHz image from Sandell et al. (2009), together with high-resolution CARMA A+B-array images at 111.1 GHz and a CARMA B-array image at 222.2 GHz. All three images show a compact core centered on IRS 1 surrounded by faint extended emission to the south, i.e., what Gaume et al. (1995) called "south spherical." At 111.1 GHz the morphology of the emission looks very similar to the 43 GHz VLA image, suggesting that all the emission is free–free emission. At 222.2 GHz there is a hint of a more east–west structure, as well as a more pronounced tail to the southwest. It is therefore likely that we start to pick up some dust emission in the interface region between the free–free outflow and the surrounding dense infalling envelope. It is unclear whether any of this dust emission is due to the accretion disk, which should be aligned east–west, i.e., orthogonal to the outflow. This faint emission is more extended in the E–W direction and there is no sign of the northwestern extension. In the south the emission is less extended and curves more to the southwest. At the highest angular resolution, i.e., the CARMA A-array images at 227.4 GHz, the faint extended emission is filtered out, and we only see a compact core with a tail toward the southwest; see Figure 4.

**Figure 4.** High spatial resolution CARMA A-array wide-band images at 227.4 GHz, both formed by shift and add and reduced independently; see Section 2.2. Both images look essentially identical. Left: December 7 continuum image, which used H30α for selfcal. The synthesized beam is 015 × 013 at P.A. = 89°. Here we also marked the positions of the Beuther et al. (2017) continuum sources cm1 and cm2. Right: December 6 continuum image, which used a 500 MHz wide continuum channel for selfcal. The synthesized beam size is 016 × 013 at P.A. = 86°.
Download figure:
Standard image High-resolution image

To determine the size of the compact structure seen at 227.4 GHz, we created images with three different weightings, robust = −2, 0.5, and 2, with the full "uv-range (40–1400 kλ), as well as restricting the short spacings to 500, 600, and 800 kλ. The IRS 1 core shrinks in size with improved spatial resolution, although the position angle stays the same, i.e., roughly north–south. The core is barely resolved for the uv-range 800–1400 kλ, where the beam size is ∼01; see Table 3. At this resolution the flux density of the compact core is 1.18 ± 0.05 Jy, or about a third of the total flux density (Zhu et al. 2013; Frau et al. 2014). Table 3 also gives the source sizes for 111 GHz and the 222 GHz B-array image. For these we determined the deconvolved source size by fitting a double Gaussian, one for the extended free–free jet and one for the compact core.

Table 3. Gaussian Fits to the IRS 1 Core

Freq.	Array	uv-range	Synthesized Beam	Size
227.4	A^a	600–1400 kλ	016 × 013 P.A. = +86°	0101 ± 0003 × 0080 ± 0003 P.A. = +5° ± 6°
227.4	A^b	600–1400 kλ	015 × 013 P.A. = +89°	0099 ± 0004 × 0090 ± 0004 P.A. = −6° ± 13°
227.4	A^b	800–1400 kλ	01 × 01 P.A. = 0°	0043 ± 0007 × 0026 ± 0007 P.A. = −15° ± 15°
222.2	B	full	029 × 023 P.A. = −85°	014 ± 001 × 009 ± 001 P.A. = +6.1° ± 8°
111.1	A	full	032 × 028 P.A. = −83°	015 ± 001 × 005 ± 001 P.A. = +0.9° ± 5°

Notes.

^aDecember 5.^bDecember 6.

Download table as: ASCII Typeset image

We also made images of the H30α emission with ∼1 and ∼10 km s⁻¹ spectral resolution. The H30α emission (line minus continuum) is more compact than the continuum emission with a deconvolved size 72 × 67 mas extended in a P.A. of −50° ± 10°. Figure 5 shows moment 0 and moment 1 images made from spectral channels above a 3σ threshold. A velocity gradient in a NE direction (P.A. ∼ 30° ± 10°) is seen in the moment 1 image at both 1 and 10 km s⁻¹ resolution. The gradient was confirmed by Gaussian fits to the velocity channels. In order to check whether the observed velocity gradient might be due to some instrumental effect, we made a spectral image for the same wide-band spectral window in the 2011 December 6 observations. No spectral gradient was seen in these data; Gaussian fits to the spectral channels in the 500 MHz spectral window agree in position within 0001. The velocity gradient of the H30α emission agrees quite well in P.A. with the velocity gradient that Zhu et al. (2013) found, 43°, in subarcsecond imaging of OCS (19–18), CH₃CN (12–11), and ¹³CO (2–1).

**Figure 5.** Moment 0 and moment 1 images of the the H30α emission made from spectral channels above a 3σ threshold. The moment 0 emission is plotted with contours on top of the moment 1 map in color. The velocity scale is shown in the bar to the right. The images were made with a spectral resolution of 10 km s⁻¹. A velocity gradient in a NE direction (P.A. ∼ 30° ± 10°) is seen in the moment 1 image.
Download figure:
Standard image High-resolution image

The distribution of the integrated line emission is orthogonal to the velocity gradient. The structure of the H30α emission is unclear, because it is not aligned with the large-scale ionized jet. The elongation and velocity gradient could arise from a prolate rotating structure or be due to infall onto or outflow from an inclined disk. The fitted major and minor axes are consistent with an inclination of 21° (zero = face-on), although this is not consistent with the outflow, which suggests that the outflow is almost in the plane of the sky. Further evidence for an asymmetric H30α distribution comes from the H30α spectrum shown in Figure 6, where a three-component Gaussian fit is shown. Images of the red- and blueshifted spectral components are separated by 0005 in R.A. and 0009 in decl. Gaussian fits show that the redshifted emission (−50 to 0 km s⁻¹) is offset to the NE in P.A. = 30° from the blueshifted emission (−110 to −60 km s⁻¹). The redshifted component is also somewhat larger (81 × 74 mas) than the blueshifted component (60 × 54 mas). The dip in the spectrum at ∼−56 km s⁻¹ might be produced by a colder envelope. Even at larger angular scales the hydrogen recombination line emission appears rather chaotic. Gaume et al. (1995) imaged H66α with the VLA at 300 au resolution. They found extremely wide lines, 250 km s⁻¹ toward the IRS 1 core, but no clear velocity gradient along the ionized jet. With the upgraded Karl Jansky VLA it would now be possible to image recombination lines of the ionized jet farther out, where the ionized gas would not be affected by accretion and where the velocity field might be more well ordered.

**Figure 6.** H30α spectrum fitted with a three-component Gaussian. The red- and blueshifted spectral components are separated by 0005 in R.A. and 0009 in decl.
Download figure:
Standard image High-resolution image

There is no sign of the three linear maser features that Goddi et al. (2015) argued to be individual H ii regions. Therefore, these are more likely to be associated with the outflow (De Buizer 2003) or photoionized preexistent clumps of molecular material (De Buizer & Minier 2005) rather than self-luminous objects. The 227.4 GHz observations with an angular resolution of 01 show only a single compact, barely resolved source, and there is no evidence for IRS 1 being a binary with components cm1 and cm2 as suggested by Beuther et al. (2017). The position of the compact core is consistent with source cm1 within positional errors and also with the single compact source detected at 219 GHz (resolution 033 × 032) using NOEMA by Beuther et al. (2018). IRS 1 could still be a binary, but if it is, the separation between the two stars must be less than 30 au if the binary components have roughly equal brightness.

IRS 1 has a total luminosity of (1.1–1.4) × 10⁵ ${L}_{\odot }$ (corrected to a distance of 2.65 kpc), as determined from extinction-corrected near/mid-IR observations (Willner 1976; Hackwell et al. 1982) and modeling of the observed dust emission surrounding IRS 1 (Akabane & Kuno 2005). We have independently checked the bolometric luminosity of IRS 1 by integrating over the observed SED, which includes accurate high angular resolution far-infrared SPIRE, PACS, and FORCAST images from 500 to 19 μm. Using published photometry from 2.2 μm to 1.3 mm, including the photometry derived in this paper (Table 1 and 2), we get a bolometric luminosity of 1.0 × 10⁵ ${L}_{\odot }$ . Here we used the PACS and SPIRE integrated flux densities in Table 2 to include heating from IRS 1 into the surrounding core, although it does not really make much of a difference, because the bulk of the luminosity comes from shorter wavelengths. This is a lower limit to the total bolometric luminosity, because we know that a significant amount of far-UV (FUV) light escapes into the outflow without heating up the cloud envelope surrounding IRS 1.

The observed luminosity corresponds to an O7 star, if the central source is a single star. The spectral type determined from free–free emission has resulted in somewhat later spectral types, O8–B0 ZAMS (see, e.g., Hackwell et al. 1982; Akabane & Kuno 2005; Zhu et al. 2013), because the ionizing flux has been estimated at a wavelength where the free–free emission is still significantly optically thick. S. Hui (2015, personal communication) derives a spectral type of ∼O5.5 from modeling high spatial resolution VLA and millimeter-interferometry data, which is more consistent with the observed bolometric luminosity.

We have refined the mass estimate for the IRS1 envelope using the published single-dish data complemented with the FORCAST and Herschel data in Tables 1 and 2. Here we used a two-component graybody fit to the observed flux densities from 1.3 mm to 20 μm, thereby avoiding contribution from the hot dust that dominates the emission at near-IR wavelength. For the cold extended envelope surrounding IRS 1 our fit predicts a size of ∼20'', a dust temperature of ∼40 K, and a dust emissivity index, β = 1.80, resulting in a total mass (gas + dust) of 490 ${M}_{\odot }$ . The fit predicts a size of 57 and a temperature of 81 K for the warm inner envelope surrounding IRS 1. The dust emissivity index for the warm component was fixed at β = 1.5. We performed a similar fit for the extended 50''–60'' cloud core surrounding the IRS 1–IRS 3 cluster and obtained similar temperatures and sizes for the warm dust component. The temperature is about the same for the cold cloud core, ∼42 K, while the dust emissivity index is slightly lower, β = 1.6, resulting in a total mass of 1300 ${M}_{\odot }$ .

4. The IRS 1 Molecular Outflow

4.1. Background

It has been widely accepted that IRS 1 drives a bipolar molecular outflow at a P.A. of ∼145° (Scoville et al. 1986; Kameya et al. 1989; Davis et al. 1998; Qiu et al. 2011), although other outflows from the young cluster surrounding IRS 1 may also contribute (Qiu et al. 2011). Kraus et al. (2006) suggested that the outflow from IRS 1 is precessing and has therefore created a wide-angled outflow. However, it is hard to think of a mechanism where the ionized outflow and the molecular outflow would be completely misaligned. Molecular outflows are secondary phenomena, which are shaped by their interaction between the wind and the surrounding cloud. It is therefore possible that the high-velocity gas close to IRS 1 may give a distorted view of the molecular outflow, especially because of the extreme accretion rate toward IRS 1. Sandell et al. (2012) noticed a large fan-shaped outflow protruding out of the molecular cloud south of IRS 1 on Spitzer IRAC images of NGC 7538. This fan-shaped nebula points back toward IRS 1. In the north there is a large, elongated cavity, also pointing back toward IRS 1. The northernmost part of this cavity could not have been excavated by the NGC 7538 H ii region, because all the O stars illuminating it are too far west of it as seen in projection, so the real separation could be even larger. The only star luminous enough to create such a cavity is IRS 1. In the 8 μm IRAC color image (Figure 7) one can see the cavity extending up to ∼250'', i.e., 3.2 pc, from IRS 1. In the south a fan shape nebulosity emerges out of the cloud and curving to the west. The tip of the nebulosity is quite bright and almost certainly a bow shock, since there is no infrared source embedded in it. This bow shock, seen at −100'', −180'' in Figure 7, is ∼200'' (2.6 pc) south of IRS 1.

4.2. The Large N–S Outflow as Seen in CO (3–2)

In Figure 7(b) we have plotted blue- and redshifted CO (3–2) overlaid in contours on the same 8 μm image seen in Figure 7(a). The strongest high-velocity emission is close to IRS 1 with the redshifted (SE) and the blueshifted (NW) emission peaks separated by ∼30'' at a P.A. of 125°–135°, which agrees with earlier studies. However, there is still strong blueshifted emission extending up to ∼80'' north of IRS 1, i.e., approximately where the large cavity becomes apparent in the 8 μm image. At the northern end of this cavity blueshifted emission again becomes visible, where the outflow interacts with the surrounding molecular cloud. The extent of the southern outflow is less clear in CO (3–2). South of IRS 1 there is both blue- and redshifted emission, suggesting that the dense molecular cloud has affected the outflow, making it largely flow in the plane of the sky, i.e., perpendicular to us. We also see the "compact" outflow from NGC 7538 South, which is located ∼80'' south of IRS 1 and within the expected path of the southern outflow. A further complication in the south is the rather strong emission from the −50 km s⁻¹ cloud (see Figure A1), which also appears to be connected to the NGC 7538 H ii region. The strong NW−SE outflow seen at low to moderate velocities in low-J CO lines is somewhat puzzling. As we will see later on, this high-velocity emission is not seen in high CO or in [C ii]. The most likely explanation is that most of this outflow emission is from young embedded low- to intermediate-mass protostars in the IRS 1 core. Especially the northwestern blueshifted emission lobe appears to point toward the protostar MM 4, as shown in the high-resolution SMA CO (2–1) image by Qiu et al. (2011).

The position–velocity plot of CO (3–2) (Figure 8) shows the northern blueshifted high-velocity gas quite well. The CO emission is very broad at the leading edge of the high-velocity cloudlets at 60''–80'' north of IRS 1, where the outflow becomes invisible, because there is no molecular gas to interact with. The outflow becomes visible again where it plunges into the surrounding cloud again at ∼240'' north of IRS 1. Here the outflow velocity is more modest, ∼10 km s⁻¹, most likely because there is still considerable mass in the surrounding cloud. To the south the outflow appears both blueshifted and redshifted. Although the position–velocity plot picks up the outflow from NGC 7538 South, this outflow cannot explain the high-velocity emission seen to the south, because the outflow from South is very compact (Sandell & Wright 2010). It is therefore likely that most of this high-velocity emission is driven by IRS 1. There is a "bow shock" seen in blueshifted emission ∼140'' south of IRS 1 where the outflow emerges out of the large NGC 7538 cloud core, i.e., the dense, massive core at ∼−58 km s⁻¹.

4.3. The IRS 1 Outflow in [C ii] 158 μm Emission

The [C ii] fine-structure line is not an outflow tracer. It is a PDR tracer (Tielens & Hollenbach 1985; Hollenbach & Tielens 1999). Although some [C ii] emission can be produced in shocks (Flower & Pineau des Forêts 2010), it is two to several orders of magnitude fainter than the PDR emission. Yet [C ii] emission is seen toward many outflows (van Kempen et al. 2010; Podio et al. 2012; Alonso-Martínes et al. 2017). In most cases the [C ii] originates from PDR emission in the surrounding cloud, but there are also some cases where one might see [C ii] emission from the cavity walls of the outflow or from UV-irradiated shocks in the outflow (Visser et al. 2012; Green et al. 2013; Alonso-Martínes et al. 2017). Figure 9 shows that the blueshifted high-velocity [C ii] emission fills the outflow lobe and looks similar to the blueshifted CO (3–2) at high velocities. We definitely see [C ii] emission from IRS 2, because it is a bright compact H ii region and the peak of the [C ii] emission is approximately centered on IRS 2. IRS 2, however, does not have an outflow. Observations of the [Ne ii] fine-structure line at 12.8 μm shows that it has ionized gas in the velocity range from −82 to 51 km s⁻¹ (Zhu et al. 2008). There is no sign of low-velocity [C ii] tracing the cavity walls of the outflow. It appears that the strong FUV radiation from IRS 1 ionizes some of the carbon in the fast-moving CO cloudlets, resulting in strong PDR emission tracing the surface layers of the cloudlets. The emission knot at 53'', −32'' seen in redshifted [C ii] (Figure 9) is not an outflow. The [C ii] channel maps (Figure A2) show that it is only seen in the channels centered at −51 and −47.5 km s⁻¹. It is unclear what star illuminates this PDR. There is not much high-velocity redshifted [C ii] emission south of IRS 1, although we see strong [C ii] emission in the velocity range from ∼−60 to −50 km s⁻¹ south and northwest of IRS 1 (see the channel maps in Figure A2). Most of this emission is likely to be PDR emission illuminated by the NGC 7538 H ii. There is no [C ii] emission associated with the compact NGC 7538 South outflow, which is very prominent in CO (11–10); see Section 4.4. However, in the the channel maps we see a compact, ∼10'' knot of emission in the channel centered on −58 km s⁻¹ about 70'' south of IRS 1 (Figure A2). The position of this knot is within 1'' of IRS 11, a well-known young Herbig Be star near NGC 7538 South. The deeply embedded IRS 11 apparently illuminates a small nebula, which we see in [C ii]. This PDR emission is also seen in the [C ii] position–velocity diagram at ∼70'' south of IRS 1 (Figure 10).

**Figure 9.** High-velocity [C ii] emission in blue and red contours overlaid on a smaller IRAC 8 μm image than in Figure 7 and plotted with a logarithmic stretch. IRS 1 is at 0'', 0'' and marked with a red star symbol. The blueshifted emission is averaged from −74 to −64 km s⁻¹ and plotted with eight linear contours going from 0.75 to 3.5 K. The redshifted emission is averaged from −52 to −42 km s⁻¹ and plotted with eight linear contours from 0.75 to 2.5 K. The redshifted peak at ∼52'', −32'' is a bright peak in the −50 km s⁻¹ cloud and not associated with any outflow. The redshifted emission NE of IRS 1 at P.A. ∼ 45° is true high-velocity emission. IRS 1, IRS 4, IRS 11, and (NGC 7538) South are labeled.
Download figure:
Standard image High-resolution image

The blueshifted emission north of IRS 1 has the same velocity structure as CO (3–2) (Figure 10), with the leading bow shock at 80'' north of IRS1 having the highest velocities. The highest velocity seen in [C ii] is −74 km s⁻¹, whereas it is ∼−85 km s⁻¹ in CO (3–2). The difference is most likely due to the lower sensitivity (shorter integration times) in [C ii]. Using observations with longer integration times, we find that the high-velocity emission has very much the same extent in both [C ii] and CO (3–2), toward IRS 1 (Figure 11). Both [C ii] and CO (3–2) show high-velocity wings of ±40 km s⁻¹ or more.

**Figure 11.** Long-integration [C ii] (red), CO (11–10) (blue), and CO (3–2) (black) spectra in native angular resolution toward NGC 7538 IRS 1. The gray vertical line in the top panel marks the systemic velocity, −58.5 km s⁻¹. In the top panel the CO (3–2) brightness temperature is multiplied by a factor of 1.5; in the bottom one, which shows the high-velocity wings, it is scaled by a factor of 10. The faint absorption between −12 and 0 km s⁻¹ is due to emission or absorption from gas in the Local Arm. In [C ii] it is absorption against the strong continuum from IRS 1. In CO (3–2) it is due to emission in the reference position.
Download figure:
Standard image High-resolution image

There are also some clear differences between [C ii] and CO (3–2). The dominant redshifted (SE) and blueshifted (NW) outflow lobes near IRS 1 are far less obvious in [C ii] than in CO (3–2); see Figures 8 and 10.

4.4. CO (11–10)

CO (11–10) was observed simultaneously with [C ii] in the L1 channel of GREAT. Since L1 is a single-pixel receiver, the CO (11–10) map is not as extended as the [C ii] map. One would expect that the emission from CO (11–10), which has a critical density of 4 × 10⁵ cm⁻³ for a gas temperature of 100 K (Yang et al. 2010), would be quite compact. Yet CO (11–10) emission extends almost all the way to NGC 7538 South, 80'' south of IRS 1 (Figure 12), which is seen even more clearly in the channel maps (Figure A3). South of IRS 1 the outflow velocities are rather modest, with line wings ∼3 km s⁻¹ in both blue- and redshifted gas to ∼40'' south of IRS 1 (Figures 12 and 13). Farther south the line is very narrow until the emission runs into the outflow from NGC 7538 South, where one sees strong blue- and redshifted emission. The southern outflow lobe appears to behave very much the same way in CO (11–10) as it does in CO (3–2) and [C ii]. The outflow is at near-cloud velocities or perhaps even slightly blueshifted, suggesting that the southern outflow lobe is close to the plane of the sky. The emission in CO (11–10) is undoubtedly excited by the outflow, because CO (11–10) is only seen along the expected path of the outflow. Only hot gas in the outflow can excite CO to such high energy levels. The velocity of the northern outflow is low and only seen to ∼40'' north of IRS 1. The high-velocity cloudlets, which stand out prominently in CO (3–2) and [C ii], are too diffuse to excite CO (11–10). There is no sign of the SE–NW outflow lobes, which dominate the outflow near IRS 1 in low-J CO lines and other low-excitation tracers, like HCO⁺.

**Figure 12.** High-velocity CO (11–10) emission in blue and red contours overlaid on an IRAC 8 μm image plotted with a logarithmic stretch. IRS 1 is at 0'', 0''. The blueshifted emission is averaged from −66 to −60.5 km s⁻¹ and plotted with eight logarithmic contours going from 0.45 to 4.8 K. The redshifted emission is averaged from −56 to −50.5 km s⁻¹ and plotted with six linear contours from 0.5 to 2.7 K. IRS 1, IRS 4, IRS 11, and (NGC 7538) South are labeled.
Download figure:
Standard image High-resolution image

**Figure 13.** Position–velocity diagram in gray scale overlaid with contours of CO (11–10) through IRS 1 at P.A. = 10°, i.e., the same as for CO (3–2) and [C ii]. IRS 1 is at offset 0''. The compact outflow from the young high-mass star NGC 7538 South at ∼81'' to the south is quite prominent in blueshifted emission. CO (11–10) emission still continues to the south. The contours are logarithmic, with 10 contours going from 0.8 to 29 K.
Download figure:
Standard image High-resolution image

4.5. The Morphology of the Outflow near IRS 1; CARMA CO (1–0) and ¹³CO (1–0) Imaging

Single-dish images do not have enough spatial resolution to show how the outflow behaves near IRS 1. We therefore also make use of some CARMA images of CO (1–0) and ¹³CO (1–0), the latter filled in for missing zero spacing with maps from the OSO 20 m telescope; see Section 2.3. The emission from the surrounding molecular cloud is very strong, and it is difficult to study the outflow at low velocities, especially in ¹²CO. In Figure 14 the lowest blue- and redshifted emission starts at ∼6 km s⁻¹ from the systemic velocity of the cloud. Some of the redshifted emission, however, is not from the IRS 1 outflow, but from gas associated with the expanding H ii NGC 7538 NW of IRS 1 (see, e.g., Davis et al. 1998), or from dense condensations in the "−50 km s⁻¹" cloud. The ¹³CO, which is more optically thin, can probe the outflow at lower velocities, but it is difficult to evaluate, especially for redshifted velocities, because of contamination from the extended "−50 km s⁻¹" cloud. This cloud has several velocity components between −52 and −45 km s⁻¹, some of which also appear to trace the PDR layer associated with NGC 7538. The cluster of submillimeter sources (Qiu et al. 2011) most likely powers molecular outflows as well. Some of the blueshifted emission NW of IRS 1 is likely from MM 4, which is shown as a filled black square in Figures 14 and 15. Some of the redshifted emission south of IRS 1 is almost certainly "contaminated" by emission from one or several of these submillimeter sources as well.

**Figure 14.** Contour plots of the low-, intermediate-, and high-velocity blue- and redshifted CO (1–0) emission from CARMA observations with 45 spatial resolution. Due to the primary-beam correction, the maps are noisier in the northernmost part of the map, because it is close to the HPBW of the CARMA primary beam. The L panel shows low-velocity emission, ±9 km s⁻¹ from the systemic velocity, −58.5 km s⁻¹, averaged over a 6.4 km s⁻¹ wide velocity interval. The I panel shows the outflow at intermediate velocities, ±15 km s⁻¹, also averaged over a 6.4 km s⁻¹ wide velocity interval, while the H panel shows high-velocity emission, ±25 km s⁻¹ from the systemic velocity averaged over a 14 km s⁻¹ wide window. At low and intermediate velocities we see blueshifted emission from the cloudlets in the blueshifted outflow lobes. At higher velocities the outflow is more jet-like and aligned with the ionized jet. The contour levels for the low and intermediate velocities are at 0.9, 3, 7, 12, 18, and 25 K, and at 0.5, 1.5, 3, 4.5, and 7 K for the high-velocity panel (right). IRS 1 is at 0'', 0'' and marked with a star symbol on the L panel. Also marked with a black square is the submillimeter and H₂O maser source MM 4 (Qiu et al. 2011), which almost certainly contributes to the blueshifted high-velocity emission NW of IRS 1. We also marked another H₂O maser with a black triangle at offset −31'', −47''. This source appears to be associated with a compact blueshifted high-velocity knot. We have also drawn the symmetry axis of the outflow and our best guess of the boundary between the outflow cavity and the surrounding cloud. These are plotted in black on each panel.
Download figure:
Standard image High-resolution image

At low velocities the ¹³CO emission is redshifted east and E–W of IRS 1. There is also some redshifted emission to the NW (Figure 15). One can also see some redshifted emission associated with the high-velocity cloudlets near the symmetry axis of the blueshifted outflow. ¹²CO (1–0) behaves the same way as ¹³CO at low and intermediate velocities, i.e., it strongly suggests that the outflow is rotating, possibly in a spiral pattern (Wright et al. 2014). What is noticeable is that at high spatial resolution the high-velocity jet or high-velocity cloudlets inside the blueshifted outflow lobe break up into compact clumps. At intermediate velocities the clumps become more centered in the middle of the outflow lobe; see panel I in Figure 14. At high velocities the blueshifted outflow is well collimated, almost N–S, like the free–free jet, and extending to ∼15'' N of IRS 1. At these velocities we no longer see the northwestern blueshifted outflow feature, which dominates at lower velocities.

4.6. What We Learned about the Outflow from CO (3–2), [C ii], and CO (11–10)

Here we summarize what we learned about the morphology of the large-scale IRS 1 outflow from examination of CO (3–2), [C ii], and CO (11–10). The large CO (3–2) map confirms that the large cavity north of IRS 1 seen in IRAC images has been created by the molecular outflow from IRS 1. We can trace the outflow in high-velocity blueshifted emission up to about 80'' north IRS 1, after which it becomes invisible owing to the low gas density of the surrounding cloud. At the tip of the outflow, ∼250'' to the north, the outflow again becomes visible as blueshifted CO (3–2) high-velocity emission where it hits the dense surrounding molecular cloud. The southern outflow from IRS 1 is probably deflected by the dense molecular cloud south of IRS 1 and is almost in the plane of the sky. Therefore, it is seen as mostly blue- and redshifted low-velocity emission. The southern outflow passes close to NGC 7538 South, another young high-mass star, and becomes invisible when the density of the surrounding cloud becomes too low to excite CO emission. The plume seen in IRAC images where the outflow emerges out of the cloud appears to contain no or very little molecular gas.

These findings are strongly supported by the smaller maps we have obtained in [C ii] and CO (11–10). The high-velocity clouds or cloudlets in the northern outflow lobe are also seen in [C ii], where they show the same morphology and velocity structure as seen in the CO (3–2) high-velocity emission. This suggests that they are illuminated by the strong UV radiation escaping from IRS 1 into the outflow. They are not seen in CO (11–10), most likely because they have too low gas density to excite CO (11–10). The outflow from IRS 1 also has modest blue- and redshifted velocities south of IRS 1 in [C ii] and CO (11–10), confirming that the southern outflow lobe is almost in the plane of the sky. However, since we trace CO (11–10) emission almost all the way to NGC 7538 South, this emission must be excited by IRS 1, since CO (11–10) has a critical density of 4 × 10⁵ cm⁻³ and an upper energy level of 304.2 K. No other mechanism can produce such warm high-density gas in this region, except radiation from IRS 1.

4.7. Outflow Properties

We use the CO (3–2) map to estimate the physical properties of the outflow, because it is the only outflow tracer we have, which covers the whole outflow. We restrict the analysis to the blueshifted outflow lobe, because it is clearly defined with very little contamination for the surrounding cloud. It is harder to extract information from the counterflow ("red"-shifted), because of the strong cloud emission and contamination from other outflows. It is not possible to directly estimate the inclination of the blueshifted outflow lobe, but it appears to be close to the plane of the sky. Here we assume an inclination of 70° ± 10°. If we use the total extent of the blueshifted outflow, 3.6 pc, a velocity of 10 km s⁻¹ (see Figure 8), we get t_dyn ∼1.3 × 10⁵ yr after correcting for inclination. This estimate has an uncertainty of about a factor of two owing to the uncertainty in the inclination. Usually dynamical timescales are underestimates (see, e.g., Parker et al. 1991). When an outflow first starts, it has to drill through dense gas, so it starts slowly. Then, once the density gets lower (and it does seek the path of least resistance), the outflow speeds up. How long this phase takes is hard to estimate. After that, it may move with relatively constant velocity as long as the gas density is relatively uniform, which is definitely not the case here. Nevertheless, here we will use the dynamical timescale derived above, i.e., 1.3 × 10⁵ yr. For estimating physical parameters of the outflow, like mass, momentum, and kinetic energy, we use the large JCMT CO (3–2) map, which is the only data set we have that covers the whole outflow. Here we only analyze the northern blueshifted outflow lobe, because it is well defined and only moderately affected by the surrounding molecular cloud, whereas the same is not true for the southern outflow, which is strongly affected by the dense molecular cloud south of IRS 1. There may be some contribution from outflows from embedded lower-mass protostars (Qiu et al. 2011), but the IRS 1 outflow will dominate overall. CO (3–2), however, is extremely optically thick, and we do not have ¹³CO (3–2) maps covering more than a small portion of the outflow. To get an idea of the ¹²CO optical depth, we therefore use the deep HARP CO (3–2) and ¹³CO (3–2) maps, which are centered on IRS 1 and cover ∼100'' × 100'' and 80'' × 80'' for ¹²CO (3–2) and ¹³CO (3–2), respectively. These maps show that ¹²CO is significantly optically thick over all outflow velocities seen in ¹³CO, or up to blueshifted velocities of >20 km s⁻¹ from the systemic velocity, here assumed to be −58.5 km s⁻¹. This is shown in Figure 16, where we show the ¹²CO (3–2) and ¹³CO (3–2) spectra in the top panel and the ratio of ¹²CO (3–2) to ¹³CO (3–2) in the bottom panel. This figure shows that ¹²CO is optically thick over the velocity range where we detect ¹³CO emission. Here we have assumed a ¹²C-to-¹³C ratio of 79 (Wilson & Rood 1994) and scaled the ¹³CO (3–2) by 7.9, i.e., the ratio would be 10 if ¹²CO (3–2) is optically thin. However, ¹²CO (3–2) is not only optically thick on IRS 1; similar ratios are seen over the whole area mapped in ¹³CO (3–2), indicating that ¹²CO is likely to be optically thick over a large part of the outflow.

**Figure 16.** Top panel: long-integration ¹²CO (3–2) and ¹³CO (3–2) spectra, the latter plotted in red, toward IRS 1. Bottom panel: ratio of T_mb(¹²CO) over T_mb(¹³CO), with ¹³CO multiplied by a factor of 7.9, i.e., a ratio of 10 corresponds to optically thin ¹²CO (3–2) emission. This shows that ¹²CO (3–2) is significantly optically thick over the whole velocity range, where we see the high-velocity gas in ¹³CO, i.e., to ∼−78 km s⁻¹, or ∼20 km s⁻¹ from the systemic velocity in the blueshifted outflow wing. At high outflow velocities the ratio is dominated by the noise in the spectra. The gray vertical line marks the systemic velocity, −58.5 km s⁻¹. There are two clouds in the direction of IRS 1, one at −58 km s⁻¹ and the other one, which is more diffuse, at ∼−50 km s⁻¹. The −58 km s⁻¹ cloud is self-absorbed in ¹²CO. Hence, the two clouds look about equal in ¹²CO, while the −50 km s⁻¹ cloud is fainter in ¹³CO (3–2).
Download figure:
Standard image High-resolution image

To estimate the outflow parameters, we therefore divided up the map into three parts. The first is a polygon encompassing the outflow around IRS 1 to 50'' north of the star, the second goes from there to ∼100'' from IRS 1, and the third one covers the tip of the outflow ∼250'' to the north. We integrate over velocity intervals of 2.1 km s⁻¹ for low velocities, then gradually increasing the velocity range to 3.8 km s⁻¹. The highest velocity range, 20.5 km s⁻¹ from the systemic velocity, is integrated over 6.4 km s⁻¹. For the area closest to IRS1 we apply the average opacity corrections derived from the ¹²CO-to-¹³CO ratio. Since we do not know what the optical depth is farther out, we take half of the correction, and at the tip of the outflow we assume that the CO is optically thin.

**Figure A1.** Channel maps of CO (3–2) averaged over 4.2 km s⁻¹ wide velocity intervals. The center velocity of each panel is indicated in the upper right corner of each panel. The emission is plotted with eight logarithmic contours from 0.35 to 26 K. The red stars in each panel are IRS 1, NGC 7538 South, and IRS 4. These are labeled in the bottom right panel (velocity = −26.3 km s⁻¹). At high blue- and redshifted velocities one only sees emission from IRS 1 and NGC 7538 South. IRS 4 is not associated with any outflow. At velocities around −60 km s⁻¹ we can see the boundary between the NGC 7538 H ii region and the molecular could south of it. Emission in the channel at −47 km s⁻¹ also appears to trace the boundary of the H ii region.
Download figure:
Standard image High-resolution image

**Figure A2.** Channel maps of [C ii] averaged over 3.5 km s⁻¹ wide velocity intervals. The center velocity of each velocity interval is indicated in the upper right corner of each panel. The emission is plotted with 10 logarithmic contours from 0.4 to 40 K. The red stars in each panel are IRS 1 and NGC 7538 South. These are labeled in the bottom right panel (velocity = −37.0 km s⁻¹). One can see widespread PDR emission from both the −58 and −50 km s⁻¹ clouds. NGC 7538 South is not associated with any [C ii] emission. The compact emission NW of NGC 7538 South appears centered on IRS 11, another young early B star in the cloud.
Download figure:
Standard image High-resolution image

**Figure A3.** Channel maps of CO (11–10) averaged over 2.5 km s⁻¹ wide velocity intervals. The emission is plotted with six logarithmic contours from 0.5 to 24 K. The center velocity of each velocity interval is indicated in the upper right corner of each panel. The red stars in each panel are IRS 1 and NGC 7538 South. These are labeled in the bottom right panel. There is strong CO (11–10) high-velocity emission from both IRS 1 and NGC 7538 South. At low velocities, panels −60.5 km s⁻¹ to −55.5 km s⁻¹, the CO (11–10) emission extends all the way from IRS 1 to NGC 7538 South. The emission extending to the west from IRS 1, which is seen in panels −60.5 km s⁻¹ and −58 km s⁻¹, is almost certainly associated with the NGC 7538 PDR rim.
Download figure:
Standard image High-resolution image

The IRS 1 molecular outflow is hot. Near IRS 1, Klaassen et al. (2009) and Zhu et al. (2013) find that the temperature of the molecular gas is 250 K or more from analysis of methylcyanide emission toward IRS 1. Since one sees CO (11–10) emission up to more than 40'' from IRS 1 and strong [C ii] emission to more than 80'' from IRS 1, the gas must be hot even farther away from IRS 1. It is therefore clear that strong UV radiation from IRS 1 shining into the outflow cavity and illuminating the gas in the outflow is capable of heating the gas over most of the outflow. Here we have assumed a gas temperature of 100 K for our column density estimates, which is the same temperature Klaassen et al. (2011) used in their study of outflows from high-mass stars. The error resulting from the uncertainty in temperature is relatively minor. A 50 K error in temperature only changes the column density by 30%–34%, whereas errors in the ¹²CO optical depth over the outflow is likely to result in an uncertainty of a factor of two or more in the mass of the outflow. We estimate the outflow mass for optically thin emission assuming a standard CO-to-H₂ abundance ratio of 10⁻⁴ (Mangum & Shirley 2015) and a mean molecular weight ${\mu }_{{{\rm{H}}}_{2}}$ of 2.8 per hydrogen molecule (Kauffman et al. 2008) and apply the opacity corrections as explained earlier in this section. By summing over all the velocity intervals from −61.4 to −82.2 km s⁻¹, we get a total mass of the blueshifted outflow lobe of 45.5 ${M}_{\odot }$ . More than 90% of the mass is the area extending up to 50'' from IRS 1. To calculate momentum, P, momentum flux, F, and kinetic energy, E_k, we use the formulae presented in Choi et al. (1993), although we correct for inclination assuming an outflow inclination of 70°. For the blueshifted outflow we get P = 290 ${M}_{\odot }$ km s⁻¹, E_k = 1.9 × 10⁴⁰ J, and F = 0.4 × 10⁻² ${M}_{\odot }$ km s⁻¹ yr⁻¹. If we assume that the outflow is momentum conserved, the mass in the redshifted outflow is likely to be similar, since the outflow velocities are about the same for both the blue- and redshifted outflows. The same is true for the other outflow parameters. If we compare our results to Qiu et al. (2011) over the same velocity interval as they used, our mass estimate, 91 ${M}_{\odot }$ , agrees quite well with their estimate, ∼100 ${M}_{\odot }$ . We can also estimate the mass of the outflow in [C ii] emission. To compute the [C ii] column density, we use Equation (26) in Goldsmith et al. (2012) assuming that the [C ii] emission is optically thin with an excitation temperature of 100 K. If we adopt 1.2 × 10⁻⁴ for the C⁺/H abundance ratio (Wakelam & Herbst 2008; see also Ossenkopf et al. 2015), we get a total mass of the [C ii] outflow of 18 ${M}_{\odot }$ , by integrating the blueshifted [C ii] emission from −74–61.5 km s⁻¹. Here we have corrected the outflow mass for contribution from IRS 2, which we know has strong blueshifted [C ii] emission (see Section 4.3), by placing a 20'' aperture centered on IRS 2. The mass of gas ionized by IRS 2 is ∼3 ${M}_{\odot }$ .

The total mass of the blueshifted outflow is therefore ∼64 ${M}_{\odot }$ , and if we assume that the mass in the redshifted outflow is the same, the total outflow mass could be as high as 130 ${M}_{\odot }$ , which is still less than the mass of the core of IRS,1, which we estimated to be ∼490 ${M}_{\odot }$ (Section 3). Compared to low-mass outflows, the IRS 1 outflow is much more massive, but the mass is similar to other high-mass star outflows (Beuther et al. 2002). There are examples of high-mass star outflows where the outflow mass even exceeds the core mass (Shepherd et al. 1998; Qiu et al. 2009). Earlier we estimated a dynamical timescale of 1.3 × 10⁵ yr, which would give us a mass outflow rate of 1.0 × 10⁻³ ${M}_{\odot }$ yr⁻¹, which agrees reasonably well with the observed mass accretion rate.

5. Discussion

The IRS 1–IRS 3 region, which borders the H ii region NGC 7538, is the center of a young star cluster with more than 150 young stars and protostars (Mallick et al. 2014; Sharma et al. 2017). The IRS 1–IRS 3 cluster is almost certainly a case of star formation triggered by the expanding NGC 7538 H ii region, although it has not been conclusively proven to be the case (Balog et al. 2004; Puga et al. 2010; Sharma et al. 2017). IRS 1 is the most massive star in the cluster, and it is still heavily accreting with an accretion rate in the range from 10⁻³ to a few times 10⁻⁴ M_⊙ yr⁻¹ (Sandell et al. 2009; Klaassen et al. 2011; Qiu et al. 2011; Beuther et al. 2012; Zhu et al. 2013), which quenches the expansion of an H ii region. It is almost certainly surrounded by an accretion disk, although there is no firm evidence of such a disk, because the emission from the ionized jet powered by IRS 1 is so strong in the millimeter and submillimeter regime that it completely hides the much fainter accretion disk.

In this paper we have analyzed high angular resolution (∼01) CARMA observations at 1.3 mm, which conclusively shows that none of the previously proposed companions to IRS 1 are high-mass stars. It could still be a binary, but then the secondary must be within 30 au, if the binary components have roughly equal brightness or if the secondary has a later spectral type. In a way this is somewhat surprising since most massive stars are close binaries or multiple systems (Sana et al. 2012). A single star is, however, consistent with the observed luminosity of the star, ≥1.0 × 10⁵ ${L}_{\odot }$ , for which the lower limit corresponds to the luminosity of a single O7 star. If the spectral type is as early as O5.5, the secondary would have to be a later spectral type. For a spectral type of O7–O5.5, the stellar mass would be in the range of 30–50 ${M}_{\odot }$ (Davies et al. 2011).

What makes IRS 1 so extraordinary is that it is an extremely young O star, which has such a high accretion rate that it quenches the formation of an H ii region. The FUV radiation by the star can only escape in the polar regions, where it excites a "jet-like" ionized outflow with velocities in excess of 250 km s⁻¹ (Gaume et al. 1995). In this paper we have also shown that the strong FUV radiation from IRS 1 heats molecular gas inside the outflow lobes and ionizes carbon, which is why we see high-velocity [C ii] 158 μm emission up to ∼100'' (1.3 pc) north of IRS 1 in the blueshifted outflow lobe. Most of this gas is still molecular and seen in CO (3–2), with similar morphology and velocity structure; see Figures 8 and 10.

We can use the jet model by Reynolds (1986) to estimate the turnover frequency and the mass-loss rate of the ionized jet, i.e., his Equations (18) and (19). In his paper Reynolds (1986) used IRS 1 as an example to illustrate of how his model could be applied. Since we now have more accurate data and knowledge of IRS 1, we can expect more robust results. From the VLA data in Sandell et al. (2009) we estimate a collimation angle of ∼15°. We further assume a wind velocity of 500 km s⁻¹, which is about what one would expect for a wind from an O star (see, e.g., Johnston et al. 2013; Sanna et al. 2016; Rosero et al. 2019), and that the jet is fully ionized. We determined the spectral index α = 0.87 and estimated the inclination of the jet to be ∼70°. Since we know the spectral index, we can estimate the rest of the jet parameters from Equations (15) and (17) in Reynolds (1986; see also Sanna et al. 2018), who analyzed radio jets from 33 radio jets. If we assume a launching radius of 20 au, i.e., the radius at which we estimated the jet to become optically thin (Section 3), we get a turnover frequency of 620 GHz, which agrees very well with our estimate from the fall-off in size, 700–1500 GHz. If the jet is launched closer to the star, the turnover frequency would be higher. Reynolds's Equation (19) now gives a mass-loss rate of 1.2 × 10⁻⁴ ${M}_{\odot }$ yr⁻¹. Although this mass-loss rate is 3–10 times higher than we found for any jet in the literature (see, e.g., Johnston et al. 2013; Zhang et al. 2019), it appears quite plausible and suggests that the jet now has enough momentum to drive the IRS 1 molecular outflow.

The core accretion models by Tanaka et al. (2016) predict that a high-mass star starts with a small thermal jet confined by the outflow from the protostar. The photoionized region is confined by the magnetohydrodynamically driven outflow. The models predict that almost the whole of the outflow is fully ionized in 10³–10⁴ yr after the initial H ii region is formed. The models do not describe the phase we see for IRS 1, i.e., when the star's ionizing flux is turned on and the accretion rate is still so high that the ionizing flux can only escape in the polar regions, i.e., into the outflow, creating the ionized jet we now see. Tanaka et al. (2017) do address the onset of ionization, but their models only go up to a stellar mass of 16 ${M}_{\odot }$ and are therefore not applicable to IRS 1, which is much more massive.

We have confirmed that IRS 1 drives a large parsec-scale outflow, which in the north extends to a projected distance of 3.6 pc from IRS 1, giving a dynamical timescale of the outflow of 1.3 × 10⁵ yr (see Section 4.7), which is probably an underestimate. In their review of massive star formation, Zinnecker & Yorke (2007) argued that high-mass star formation is not a scaled-up version of low-mass star formation. It requires a cluster environment with a much larger mass reservoir and higher pressure than in a low-mass star-forming environment. However, in the very earliest stages of the formation of a high-mass star it may be very difficult or impossible to judge whether the protostar will form a low-mass star or evolve into a high-mass star. Based on the dynamical timescale, IRS 1 most likely started as a much lower mass star. Although there is no way to estimate the accretion rate in the past, it probably started with a much lower accretion rate. The accretion rate increased in time, when the mass of the core in which it was embedded grew in size, i.e., essentially what is often called competitive accretion, but perhaps better called runaway accretion (Zinnecker 1982) or hierarchical accretion (Larson 1978). IRS 1 will continue to grow in mass until the accretion rate drops, so that it can no longer quench the H ii region. When this happens, the H ii will rapidly expand and disperse the material around it. There are even two other examples of high-mass stars in NGC 7538 that are going through the same process, IRS 9 (Sandell et al. 2005; Rosero et al. 2019) and NGC 7538 South (Sandell & Wright 2010). Both IRS 9 and NGC 7538 South are younger than IRS 1 and have similar accretion rates to IRS 1. Both have currently weak radio jets and luminosities suggesting that they are still early B stars. Given the current accretion rates, they will rapidly evolve into O stars like IRS 1.

6. Summary

NGC 7538 IRS 1 is a very young O star at a distance of 2.65 kpc, which is still so heavily accreting that the accretion prevents the formation of an H ii region. The new observational data presented in this paper provide new insight on IRS 1 and suggest that at least some O stars form the same way as low-mass stars, i.e., with an accretion disk and driving an outflow, but with much higher accretion rates.

The data presented in this paper confirm the model proposed by Sandell et al. (2009), who showed that the free–free emission from IRS 1 is dominated by a collimated ionized wind driving an ionized north–south jet. Reanalyzing with new data from the literature and from this paper, we find that the radio SED follows a power law with a spectral index α = 0.87 ± 0.03, which is very similar to what Sandell et al. (2009) found. Only at frequencies higher than 300 GHz can one start seeing some excess emission due to thermal dust emission. The size of the free–free core falls off as ν^{−0.92 ± 0.02}, which is consistent with a jet (Reynolds 1986). If we assume that the jet is launched from a radius of 10–20 au, this fit predicts that the free–free emission will not become completely optically thin until somewhere between 700 and 1500 GHz. Mapping of the outflow from IRS 1 in [C ii] and CO (11–10) confirms that strong FUV radiation leaks out in the polar regions of the IRS 1 core (or disk), ionizing carbon and heating up the molecular gas.

We find no evidence for IRS 1 being a binary or triple system as has been suggested in the past. Analysis of high angular resolution (∼01) CARMA observations at 1.3 mm shows a barely resolved single source. It could still be a binary, but if it is, the separation between the two stars must be less than 30 au if the binary components have roughly equal brightness. With our spatial resolution and sensitivity we would have easily detected the two binary components proposed by Beuther et al. (2017).

We find that IRS 1 drives a massive parsec-scale outflow, which in the north has sculpted a cavity extending to 3.6 pc from IRS 1. To the south it is more difficult to trace the outflow, because it penetrates the massive molecular cloud south of NGC 7538. IRAC images suggest that it terminates at ∼2.6 pc from IRS 1. At high velocities the northern outflow is well aligned with the free–free outflow. We trace high-velocity cloudlets inside the outflow cavity to at least 80'' north of IRS 1 in both CO (3–2) and [C ii]. Since CO and [C ii] show the same morphology and velocity structure, they originate in the same gas, i.e., [C ii] traces the outer layers of the gas ionized by the strong FUV field from IRS 1 shining into the cavity. We do not see these cloudlets in CO (11–10), which is a tracer of high-density hot gas, but we see CO (11–10) all the way down to NGC 7538 South. This gas can only be heated by FUV radiation from IRS 1. We estimate a dynamical timescale of the outflow of 1.3 × 10⁵ yr and a total outflow mass of 130 ${M}_{\odot }$ , which would give us a mass outflow rate of 1.0 × 10⁻³ ${M}_{\odot }$ yr⁻¹. This agrees quite well with the observed accretion rate, which is in the range from 10⁻³ to a few times 10⁻⁴ ${M}_{\odot }$ yr⁻¹, especially since the dynamical timescale is likely to be underestimated.

We find that the O star IRS 1 was formed in a similar way to low-mass stars, i.e., with an accretion disk and driving an outflow, but in a much denser environment and with a much higher accretion rate. It is likely that many other high-mass stars may form the same way, although other formation scenarios are also possible. In the NGC 7538 molecular cloud there are also two other high-mass stars: IRS 9 and NGC 7538 South, which both appear to be younger analogs of IRS 1. Therefore, at least some high-mass stars form in a similar fashion to low-mass stars.

Based in part on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. The development of GREAT (German REceiver for Astronomy at Terahertz frequencies) was financed by the participating institutes, by the Federal Ministry of Economics and Technology via the German Space Agency (DLR) under Grants 50 OK 1102, 50 OK 1103, and 50 OK 1104 and within the Collaborative Research Centre 956, subprojects D2 and D3, funded by the Deutsche Forschungsgemeinschaft (DFG). Support for CARMA construction was derived from the Gordon and Betty Moore Foundation, the Kenneth T. and Eileen L. Norris Foundation, the James S. McDonnell Foundation, the Associates of the California Institute of Technology, the University of Chicago, the states of California, Illinois, and Maryland, and the National Science Foundation. Ongoing CARMA development and operations were supported by the National Science Foundation under a cooperative agreement and by the CARMA partner universities. We acknowledge support from Onsala Space Observatory for the provisioning of its facilities/observational support. The Onsala Space Observatory national research infrastructure is funded through Swedish Research Council grant No. 2017-00648. We thank Dr. Hans Zinnecker for constructive criticism of the paper. We thank the referee for comments that improved this manuscript.

Appendix: Channel Maps

In Figures A1–A3 we present channel maps of the JCMT CO (3–2) map and the GREAT [C ii] and CO (11–10) maps discussed in the main body of the paper.

NGC 7538 IRS1—an O Star Driving an Ionized Jet and Giant N–S Outflow

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction