Disks and Outflows in the Intermediate-mass Star-forming Region NGC 2071 IR

Among these sources, HOPS-361-A (hereafter IRS1) and HOPS-361-C (hereafter IRS3) have the strongest 1.3 mm continuum, with peak intensities of 152 and 114 mJy beam⁻¹, respectively. As stated in Section 1, IRS1 and IRS3 are the dominant mid/far-IR luminosity contributors and also presumed driving sources of the large-scale outflow (e.g., Torrelles et al. 1998; Eislöffel 2000). As shown in Figure 2 (see also Figure 3), IRS1 is partly resolved at 1.3 mm, and no clear elongation is apparent. At 0.87 mm IRS1 appears better resolved. The inner brighter part of IRS1 (i.e., flux intensity above 60 mJy beam⁻¹) has a bar-like shape, with an elongation at P.A. of about 25°. This elongated structure is further embedded in low-level extended emission (below 60 mJy beam⁻¹ but still above 20σ = 11 mJy beam⁻¹). This weaker component is approximately elliptical, and its major axis extends about 03 along the NW-SE direction, i.e., ∼110° offset in orientation from the inner bright component. There are some hints of spiral-like bending features at the interface between the two components and may further connect with larger-scale spiral-like features extending to ∼1'' (see Figure 3). These features have added to the complexity on inferring the configuration of the protostellar disk of IRS1. The kinematic information, which is discussed in the following sections, is more supportive of a disk major axis oriented in the NW-SE direction, i.e., consistent with the extended component. In this scenario, the inner bright component revealed in 0.87 mm appears as an unusual substructure of the protostellar disk. The deconvolved FWHM from Gaussian fit to the 0.87 mm continuum emission, which is dominated by the inner component, is 025 × 015 (108 au × 65 au).

Zoom In Zoom Out Reset image size

The VLA 9 mm continuum, which traces both free–free emission and thermal dust emission, shows a distinctly different morphology with respect to the ALMA images. The 9 mm continuum emission of IRS1 appears as a marginally resolved condensation, which coincides well with the position of the 0.87 mm flux peak. The 9 mm emission has a T-shaped elongation, i.e., with the brighter part extending along the NE-SW direction (P.A. ∼ 25°), and a fainter part extending slightly to the east (Figure 3). The NE-SW extension has a direction similar to the brighter part seen in 0.87 mm. In addition, some low-level (∼5σ) diffuse emission is also seen to the east and west of the central source, which extends as far as 05 (see Figure 3(c)). Trinidad et al. (2009) reported radio knots ejected from IRS1 along the E-W direction (IRS1E, IRS1W) based on the VLA 1.3 cm continuum. Comparing with their detections, the weak diffuse emission in the 9 mm could also trace a radio jet in the E-W direction. It is likely that some of the previously detected radio knots, like IRS1W, have dissipated most of their energy and can no longer be detected, thus absent in our map, though the different surface brightness sensitivities in the two observations may hinder a conclusive interpretation.

On the other hand, IRS3 shows a clear disk at 0.87 mm and 1.3 mm. Gaussian fits to the continuum in both bands give a similar deconvolved size of 0.48'' × 0.19'' with a position angle of 130°. Assuming that the Gaussian semimajor axis corresponds to the disk radius, it has a radius of 103 au, and the inclination can be estimated to be 67° by assuming that it is a geometrically thin disk and then calculating the inverse cosine of the minor axis divided by the major axis. The flux distribution in the 0.87 mm continuum appears asymmetric, with the emission peak offset by 0.16''(∼69 au) to the northwest compared with the geometric center. A similar asymmetry could be present at 1.3 mm, but this is less clear due to the lower spatial resolution. Interestingly, in the center of the disk, the 9 mm continuum further reveals a binary system separated by 01 (∼43 au). The two components (IRS3A, IRS3B) have a flux ratio of ∼10 at 9 mm, and the more luminous component, IRS3A, is coincident with the geometric center of the disk structure seen at 0.87 mm and 1.3 mm. IRS3B is located to the northwest of IRS3A and closer to the emission peak at 0.87 mm. IRS3B could be (at least partly) contributing to the asymmetric flux distribution seen at 0.87 mm via enhanced heating toward the surrounding dust/gas. While the detection of IRS3B is a point source, IRS3A is resolved and extends along the NE-SW direction to about 024 (∼100 au) on both sides, with a position angle of ∼15°. Extension from IRS3 in this direction has been reported in earlier VLA 1.3 cm and 3.6 cm observations, albeit with a lower resolution (Carrasco-González et al. 2012; Trinidad et al. 2009), and interpreted as a radio jet.

Figure 4 presents the integrated intensity map of a selected number of spectral lines toward IRS3. The line emission integrated over two velocity intervals (1.5 km s⁻¹ < ∣v − v_sys∣ <6.5 km s⁻¹, relative to a systemic velocity v_sys ≈ 9.5 km s⁻¹) is shown in blue and red, respectively. Almost all lines exhibit a clear velocity gradient along the major axis of the millimeter continuum, with emission transitioning from blueshifted in the northwest to redshifted in the southeast. This monotonic velocity transition and its correspondence with the dust continuum are strongly indicative of a Keplerian rotating disk. Figure 4 also reveals the difference in spatial distribution of line emission from different molecules. Species including C¹⁸O, ¹³CO, H₂CO, and SO exhibit strong emission beyond the disk boundary defined by 1.3 mm/0.87 mm dust continuum. The line emission from ¹³CH₃OH, CH₃OH, SO)₂, and other organic molecules are more spatially compact, i.e., within 025 (∼108 au) from the center. Hereafter we refer to the two groups of lines with distinct morphology as group A and group B, i.e., group A lines include those from C¹⁸O, ¹³CO, H₂CO, and SO, while group B lines include transitions from CH₃OH, SO₂, ¹³CH₃OH, as well as other complex organic molecules (CH₃OCHO is shown here as an example).

Zoom In Zoom Out Reset image size

These two categories are better illustrated in the PV diagram extracted along the major axis of the dust continuum, as shown in Figure 5. The group A lines, i.e., C¹⁸O, H₂CO, SO, etc, show bright emission peaks in the first and third quadrant. For ¹³CO and H₂CO, the detection in first quadrant (i.e., blueshifted emission) is stronger. For all species in group A the detection close to the systemic velocity is relatively weak, possible due to self-absorption of cold gas along line-of-sight and/or spatial filtering of the extended emission from the ambient cloud in interferometric observations, especially for C¹⁸O and ¹³CO. In contrast with group A, lines in group B appear as a continuous linear feature crossing the first and third quadrant, and no low-velocity emission extending beyond 05 (∼215 au) is apparent. The linear feature is consistent with a velocity gradient around 20 km s⁻¹ arcsec⁻¹, or 0.046 km s⁻¹ au⁻¹, and the intensity distribution across it is relatively uniform. The spatial extents of these lines line up well with the disk boundary inferred from the 1.3 mm/0.87 mm continuum. Nevertheless, the outer emission edges on the PV diagram have a convex shape, in contrast with the expectation for a Keplerian disk, but we will show in Section 5 that this is mainly due to the limited spatial resolution. Our band 6 observation has a resolution of ∼025 (108 au), comparable with half of the major axis of dust continuum, so the detailed PV structures have been smoothed out in this plot, especially for group B lines. The different spatial and kinematic distributions in the two groups are likely to be reflecting different excitation conditions required for different transitions, e.g., group A and B lines have systematically different upper energy levels (see Table 1).

Zoom In Zoom Out Reset image size

The kinematics of IRS1 are more difficult to infer as we do not have clear knowledge about the disk orientation from the dust continuum. Figure 6 presents the integrated intensity map of spectral lines toward IRS1. There is an extended blueshifted structure, about 06 to the west of IRS1, seen in C¹⁸O, ¹³CO, H₂CO, SO, and CH₃OH. This feature is not associated with the inner disk of IRS1 and should be tracing the extended emission adjacent to it (see Figure 1 or Figure 3). Interestingly, for most species there appear to be one blueshifted clump (cB1) and two redshifted clumps (cR1, cR2) associated with IRS1. This is most clear for H₂CO, SO, and CH₃OH. For an individual rotating disk, one would expect a monotonic velocity gradient along the major axis, as observed for IRS3. Here we think cB1 and cR1 are tracing the protostellar disk, while the third gas clump cR2 is a separate structure that is not associated with IRS1 disk, for the following reasons. First, the position of cR2 is more spatially offset from the emission peak of dust continuum compared with cB1 and cR1. If cB1 and cR2 are tracing the gas rotation on both sides of a disk, then the inferred position of a disk will disagree with that traced by the dust continuum. Second, not all the lines exhibit clear detection at the position of cR2. The cR2 clump is absent in organic molecules like CH₃OCHO, and for SO₂ and ¹³CH₃OH only some weak extension from cR1 toward cR2 is seen. Again, this is in contrast with the expectation that cR2 is tracing one side of the disk as similar chemical/excitation properties are generally expected for both sides of a disk. Therefore the IRS1 disk, as traced by cB1 and cR1, is oriented in the NW-SE direction. The other peak cR2 is likely a shock knot excited by the ejection from the protostar.

Zoom In Zoom Out Reset image size

Figure 7 illustrates the PV diagram extracted along the inferred disk orientation of IRS1, i.e., following the direction of P.A. ∼ 135°. Similar to IRS3, there are two categories of molecular lines: C¹⁸O, ¹³CO, H₂CO, and SO show low-velocity emission extending beyond 05, while ¹³CH₃OH, CH₃OH, SO₂, and other organic molecules are exclusively tracing the inner disk. However, unlike the prototypical case of IRS3, for the lines in group B, the linear feature is composed of two separate emission peaks in the first and third quadrants, instead of a more continuous distribution. In addition, the PV diagram of IRS1 appears more asymmetric against the origin, in terms of both the intensity and the shape. This deviation from symmetric kinematics may arise from an imperfect determination of the disk position/orientation, or confusion by other mechanisms, like ejection, or hidden multiplicity inside IRS1.

Zoom In Zoom Out Reset image size

In light of this we add an envelope component starting from R_disk and extending to an outer boundary R_out. We assume the envelope has a flat geometry that is similar to the disk, i.e., density ρ = 0 for h > 0.2 × r. The envelope starts from the centrifugal barrier at radius R_disk = R_cb and extends to an outer boundary R_out. We also assume the density changes smoothly from the disk to the envelope, and the envelope follows a density distribution a ρ ∝ r^−1.5. This corresponds to the typical density profile of an infalling cloud (e.g., Shu 1977; Harvey et al. 2003). The envelope has the motion described by

$\begin{eqnarray}&&{v}_{r}=-{v}_{\mathrm{cb}}\displaystyle \frac{\sqrt{{R}_{\mathrm{cb}}(r-{R}_{\mathrm{cb}})}}{r},\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray}&&{v}_{\phi }={v}_{\mathrm{cb}}\displaystyle \frac{{R}_{\mathrm{cb}}}{r},\end{eqnarray} \tag{ 9 }$

where v_cb = $\sqrt{2{{Gm}}_{* }/{R}_{\mathrm{cb}}}$ is the rotation velocity at the centrifugal barrier. Such motion conserves both angular momentum and mechanical energy (see, e.g., Sakai et al. 2014). An example of the adopted velocity profiles is shown in Figure 15.

Zoom In Zoom Out Reset image size

The same MCMC routine is run (with one more free parameter R_out). The resultant best-fit parameters are shown in Table 6 (disk + envelope (A)). The PV diagram of the best-fit model, which gives m_* = 1.24 M_⊙, R_disk = 113 au, and R_out = 394 au, is shown in panel (b) of Figure 14. This model has a better performance in fitting the low-velocity emission, and the best-fit disk radius R_disk also agrees well with the measurements in Section 5.2. This agrees with the explanation that group A lines (at least for SO) are tracing both the disk and the inner envelope components.

Interestingly, when an envelope component is added, the best-fit stellar mass, 1.24 ± 0.04 M_⊙ is smaller than the best-fit masses from modeling group B lines, i.e., 1.4–1.5 M_⊙. This underestimation of the dynamical mass could be partly attributed to the oversimplified envelope kinematics in the calculation. For example, the presented envelope model has ignored the motions of the infalling material in the z direction. In particular, the tangential velocity profile is not continuous at the location of R_disk (the radius of the centrifugal barrier), i.e., there is a jump by a factor of $\sqrt{2}$ when transitioning from the disk to the envelope. Such velocity jump may not be physically realistic and has not been suggested in the observations. Mathematically this discontinuity in the velocity profile could have led to the underestimation in the central mass, as a smaller mass is needed to compensate for the large tangential velocity beyond R_disk.

To test this we run the MCMC routine with a modified velocity profile, in which we assume the envelope to disk transition takes place at the centrifugal radius, i.e., R_disk = R_c =2R_cb to prevent the rotation velocity from increasing inside R_c and exceeding the Keplerian velocity. In this way the transition of the rotation velocity profile becomes continuous (see Figure 15). The best-fit results and parameters are presented in Figure 14 and Table 6 (disk + envelope (B)). As can be seen, this model returns m_* = 1.44 ± 0.09 M_⊙, which is consistent with the measurements using compact disk tracers in Section 5.2.

In order to further explore how the fitting is affected by different envelope models, we also try to fit the observations with Ulrich (1976; hereafter Ulrich76) envelope model. It is based on the solution of the collapse of a spherically symmetric cloud in uniform solid-body rotation when the pressure forces are negligible, in which material infalls following ballistic trajectories. With the assumption of solid-body rotation, particles falling near the rotational axis have smaller angular momentum and will fall in close to the central star, while particles falling in from regions near θ ∼ π/2 will fall in at a maximum centrifugal radius R_c , which is determined by the specific angular momentum measured around the rotation axis. The material arriving at the midplane will collide with material arriving at the same position from the opposite z direction, thus producing a flat disk (after the kinematic energy dissipates in shocks). For simplicity we ignore the detailed processes of disk formation and assume a thin disk undergoing pure Keplerian rotating motion inside R_c . Therefore, for r < R_disk (=R_c ), we keep the same geometry (h = 0.2 × r), density, and kinematic setup for the disk component as in other models, and for r > R_disk, we adopt the ballistic solution as in Ulrich76 to set the density and velocity for each model grid (the density is evaluated by assuming that the mass infall rate is steady). However, there is a singularity in the density solution at r = R_c and z = 0 (on the midplane). This will not cause numerical issues as long as no model grid is centered on this point, but we cannot assume a continuous density transition from the disk to the envelope as other models above. Instead we add another parameter, the mass ratio between the envelope and disk m_env/m_disk, to properly characterize the density contrast between the two components.

The results are also shown in Figure 14 and Table 6. The best-fit model gives m_* = 1.43 M_⊙, R_disk = 153 au, and R_out =347 au. Unlike the envelope model above, which shows a concave PV feature (i.e., at the second/fourth quadrant in the PV diagram), the Ulrich76 envelope model seems to predict a more compact and diamond-shaped PV diagram, in contrast with the observation. This appearance could partly arise from the difference in the assumed envelope geometry: in the simplified envelope case we have assumed a flat geometry (i.e., density ρ = 0 for h > 0.2 × r), which is consistent with a relatively late evolutionary stage when the surrounding material has been partly evacuated; in the Ulrich76 model there is still material distributed in smaller polar angles, so the infalling material close to the rotational axis could contribute more low-velocity emission, seen in the second/fourth quadrant of the PV diagram.

In summary, we attempt to model the SO 6₅ − 5₄ line with different setups, including a pure Keplerian disk and a disk plus an inner envelope, and different envelope models are also investigated. In particular, the disk + simple envelope (A) model, in which we adopt the radius of the centrifugal barrier as the transition point of the disk and envelope, appears to lead to an underestimation of the central stellar mass, mainly due to a jump in the rotation velocity profile. This suggests that the simplified thin envelope model based on the centrifugal barrier assumption may be oversimplified in describing the kinematic transition between the disk and envelope.

Overall the model will have better performance in reproducing the low-velocity extended emission when an envelope component is included, though the mass estimation seems to rely on the specific envelope kinematics, especially at the disk-to-envelope transition point. In this case, group A lines like SO arise from both the disk and part of the envelope, while group B lines exclusively trace the disk. However, the disk-only model also gives a reasonable fit, which can be potentially improved if a different disk geometry or different density profiles are to be used. In this case, group A lines and group B lines are both tracing the disk, but group A lines like SO 6₅ − 5₄ have a more widespread distribution, likely due to their relatively low upper energy level (see Table 1). It is difficult to distinguish the two different scenarios based on current observations. Higher-resolution observations and further exploration of the envelope/disk properties, such as the geometry and density contrast between the disk and envelope, are needed to clarify the issue. Given the systematic uncertainties involved in modeling the group A lines, we rely on the dynamical mass estimated for group B lines (compact disk tracers) in Section 5.2 for related discussions.

For IRS3 our kinematic modeling gives a central mass estimation of 1.43–1.52 M_⊙ from fits to different molecules. The variation is probably reflecting systematic differences of molecules in the spatial distribution within the disk, due to the variation in abundances and excitation conditions. More detailed modeling involving physically realistic disk properties and astrochemical evolution is required to better understand the variations in different molecules but is beyond the scope of this paper. Here we take the range 1.43–1.52 M_⊙ as a reasonable estimation for possible central masses of IRS3. This estimation is further consolidated by our SED fitting in Section 4, which favors a central mass around 2 M_⊙. Note that our VLA 9 mm observations have revealed the multiplicity in IRS3. In principle, if the molecular lines are tracing Keplerian orbiting motions around the two components (IRS3A, IRS3B), then the estimated dynamical mass should be treated as the sum of two.

Our kinematic modeling also has the ability to constrain the mass ratio of IRS3A/IRS3B. In Section 5 we attempted the modeling of IRS3 with the position offset as a free parameter to allow for a precise measurement of the "kinematic center", which corresponds to the binary barycenter as the disk kinematics are regulated by the center of mass. In the fixed-i case, the best-fit offsets range from −13.60 to −19.54 au for different lines, indicating that the kinematic center is about 13.60–19.54 au to the NW direction compared with the reference point, as illustrated in Figure 16. Therefore the kinematic center is located roughly in the middle of IRS3A/B, and one can further estimate a IRS3A/IRS3B mass ratio of 1.3–2.5 by comparing their relative distances from the kinematic center.

Zoom In Zoom Out Reset image size

It is also interesting to note that IRS3A appears to coincide with the "geometric center" of the disk. Here the "geometric center" can be defined as the center of symmetry of the elliptical contours in the 1.3 mm/0.87 mm continuum at relatively low levels, where the emission is not obviously skewed to the NW direction. In Figure 16 we marked the position of the intensity-weighted center for the outer disk (defined by the region between 100σ–400σ isophotal contours in 1.3 mm). If the disk is in a steady state with its material distribution regulated by the central mass, one would expect its geometric center to be close to the barycenter of the binary system. However, in our case there seems to be a nonnegligible deviation between the two.

However, we note that our spatial resolution is only 023 (∼100 au) in band 6, so it is challenging to infer the position of the kinematic center (or the relative position between ALMA/VLA detections) with a few au precision. One potential issue is the limitation in absolute positional accuracy, which hinders precise determination of the relative location between the VLA/9 mm and ALMA/1.3 mm detections. The theoretical accuracy limit imposed by the S/N can be estimated by δ p = θ_beam/(S/N)/0.9 (e.g., Cortes et al. 2021), which is around 0.3 mas for a FWHM beam size of θ_beam = 023 and S/N = 890 for IRS3 in band 6. However, the atmospheric phase fluctuations limit the δ p to about 0.05θ_beam (Cortes et al. 2021), which is around 002, or 5 au. Second, the determination of the kinematic center relies on the assumption of symmetrically distributed line emission on both sides of the disk. Nevertheless, we have observed an asymmetric intensity distribution in the high-resolution 0.87 mm continuum, which is heavily skewed to the NW direction. If the asymmetry arises from the local enhanced temperature or density (likely due to the existence of IRS3B), then the molecular line emission could be stronger toward the NW side as well, and thus leading to a deviation of the derived kinematic center to the NW direction. There are also some hints of asymmetry in the line intensity distribution in some lines (see, e.g., Figure 5). Follow-up observations in the future are required to clarify this.

The mass of IRS1 is less well constrained mainly due to a lack of knowledge of the inclination angle. In the case of a fixed inclination of 60°, the kinematic modeling gives a central mass in the range of 4.40–5.46 M_⊙ with scattering from different molecules. Systematically larger (smaller) masses can be obtained if smaller (larger) inclination angles are assumed. The free-i cases favor a modestly large inclination of 56°–70°.

One can also argue against a very small inclination of i ≲ 30°(close to face-on configuration) for IRS1 based on the observed morphology/kinematics of outflows driven by IRS1. The H₂ 1-0 S(1) map reveals a strong outflow associated with IRS1 roughly in the E-W direction, with the eastern lobe extending as far as 30'' (∼0.06 pc; outflow II in Eislöffel 2000; see also Walther & Geballe 2019). In our CO 2–1 observations the outflow associated with IRS1 mainly appears as a single blueshifted lobe located to the southwest of the source, in contrast with the face-on configuration for which one would expect more spatially overlapping blueshifted and redshifted line emission. Therefore the IRS1 disk should have a moderate or higher inclination, although we cannot be more certain about the precise range.

The SED analysis in Section 4 provides a good constraint on the upper limit of the central mass of IRS1. The best five fit models of IRS1 have a stellar mass m_* ≲ 4 M_⊙. A larger stellar mass will typically result in larger fluxes in mid-/far-IR and bolometric luminosities compared with the observation. The best-fit model with the largest stellar mass (i.e., m_* = 8 M_⊙) has χ² ≳ 10, significantly worse compared with that of the minimum χ² (∼5), and hence is highly unlikely. Nevertheless, the ZT model grid is rather sparse for smaller m_*, and it is difficult to derive a more precise upper limit from the SED fitting. Combining the information from both kinematics and SED, we speculate the stellar mass of IRS1 is about 3–5 M_⊙. This puts constraints on the inclination angle of IRS1, i.e., i ≳ 60° according to Table 5. A caveat is that in the SED analysis, we have implicitly assumed that IRS1 is dominated by a single protostar while multiple components could exist inside IRS1.

The CO 2–1 data reveal a high-velocity bipolar jet associated with IRS3, with a position angle in a range of 22–32° seen in different velocities. This further confirms the association of the large-scale NE-SW H₂ outflow with IRS3. Interestingly, there is a clear misalignment between the jet/outflow seen in different tracers: the H₂ outflow IA/IB has an average position angle of ∼45°, with its lateral extents covering a wide range of around 10°–60°. The molecular jet shows a position angle in a range around 22–32°. The radio jet in our 9 mm map has a position angle ∼15° (which is consistent with measurements in the literature; see Torrelles et al. 1998; Seth et al. 2002; Trinidad et al. 2009; Carrasco-González et al. 2012). Furthermore, the disk has a position angle of ∼130°, approximately perpendicular to the high-velocity CO jet (but not the radio jet). The most likely mechanism to account for these observational phenomena is a precessing jet wiggling over a range of jet position angles, and thus the various orientations seen in different tracers represent the interaction of the jet and environment material at different phases. Jet precession has long been known in both low-mass and high-mass YSOs (see reviews in Frank et al. 2014; Lee 2020), and large axis changes of up to ∼45° are also observed in a few sources (e.g., Cunningham et al. 2009).

The current radio jet seen in the centimeter continuum indicates the most recent ejection event from the protostar. It is also suggested in Carrasco-González et al. (2012) that the jet is precessing as they found variations in the jet orientation (a few degrees) of the IRS3 jet over a few years (1995–1999). However, their observations have relatively low resolutions (∼0.3'' in 3.6 cm), and hence the measurements may suffer from more uncertainties due to the imperfect Gaussian fitting, or contaminated by the unresolved binary component. Our high-resolution 9 mm data give a more unambiguous measurement of the jet position angle, i.e., ∼15°, which is close to the value measured for the 1998/1999 epochs in Carrasco-González et al. (2012). There is no smoking gun evidence for the changes in the radio jet orientation in a timescale of a few years, and further observations with similar spatial resolutions are needed to confirm/clarify it and better constrain the precession period.

In some velocities the CO outflow appears as misaligned segments, and these features are roughly symmetrically distributed in the blueshifted and redshifted lobes (e.g., at ∣v − v_sys∣ = 30 km s⁻¹; see Figure 8). Such point-symmetric (i.e., S-shaped) wiggles are expected in the precession of accretion disks because of tidal interactions in noncoplanar binary systems (see, e.g., Raga & Biro 1993; Terquem et al. 1999). Similar misaligned segments and S-shaped wiggles are also tentatively detected in the H₂ emission (see figures in Eislöffel 2000; Walther & Geballe 2019), but the spatial distribution of H₂ is more complex and may involve other physical processes including deflection or blocking of part of an outflow. At lower velocities we detected a wide-angle component of the CO outflow. Interestingly, the edges of this component line up well with the spatial extent of the H₂ emission. This further consolidates a unified origin of the CO and H₂ emission and suggests that this wide outflow cavity is at least partly shaped by the jet precession.

Combining the information in the literature and this work, IRS1 exhibits a similar complexity on the inferred ejection directions. Trinidad et al. (2009) detected a few ejected condensations (IRS1E, 1W) from IRS1 along the E-W direction (P.A. ∼ 100°) based on VLA 1.3 cm and 3.6 cm observations (see also Carrasco-González et al. 2012). In the H₂ emission, IRS1 seems to drive a east-oriented outflow (P.A. ∼ 70°). This outflow contains several shock knots, and some interspread diffuse emission and are collectively called outflow IIA in Eislöffel (2000), but there could be contribution from other protostars like IRS2. We have also detected some blueshifted CO emission that is likely associated with outflow IIA. The western lobe (IIB) is much fainter in H₂. Note that, given the disk orientation (P.A. ∼ 135°) of IRS1, it it likely that the large-scale outflow in the NE-SW direction seen in single-dish CO observations (e.g., Stojimirović et al. 2008) is partly contributed by IRS1. On the other hand, a smaller P.A. (∼45°) is required to account for the low-velocity CO outflow, especially the blueshifted wide-angle component to the southwest. There are some relatively high-velocity CO knots distributed close to the E-W orientation as well. The situation is further complicated by the nonideal disk appearance seen in the high-resolution 0.87 mm, 1.3 mm, and 9 mm continuum. The relatively diffuse part of the 0.87 mm emission appears consistent with a NE-SW oriented disk (P.A.∼135°) but it contains a bright inner part elongated at around P.A. ∼ 25°. Extension along a similar direction is also seen in the 9 mm continuum, as well as a protuberance to the east. This east-oriented extension could be due to an E-W jet, as suggested in Trinidad et al. (2009) based on ejected radio knots, whereas the extension with P.A. ∼ 25° is less clear without auxiliary information.

In summary, the observed radio condensations or high-velocity CO suggest a jet along or close to the E-W direction, while the disk kinematics and CO outflow agree with a NE-SW-oriented ejection direction. The H₂ emission has a spatial content with a P.A. range approximately in between but may have contributions from other YSOs. Similar to the case of IRS3, one might consider a precessing jet to account for the observed ejection events from IRS1. Carrasco-González et al. (2012) reported that the direction of the IRS1 jet appears to be changing slightly over a few years based on the 3.6 cm continuum morphology. In particular, the jet slightly bends to the northwest around 05 to the west of the protostar. Carrasco-González et al. (2012) argued that it could be due to either the superposition of a binary jet or a single jet interacting with the ambient medium. The latter scenario is reminiscent of our 1.3 mm continuum map, where we have also observed a dust clump around 08 to the west of IRS1. This dense clump could be responsible for the bending feature in the radio jet. Optionally, we note that a precessing jet can also naturally explain the variations in the jet orientation and morphology.

The possibility of unresolved multiplicity cannot be ruled out with current data. Carrasco-González et al. (2012) suggests that the extension in their 0.09'' resolution 1.3 cm map can be interpreted as a marginally resolved close binary system with a separation of ∼40 au. Our 9 mm map of IRS1 exhibits a similar protuberance to the east. A secondary component could be responsible for the origin of precession. Alternatively, one could assign the observed ejections in different directions to different components in the binary system, e.g., one along E-W and the other along NE-SW. However, more aggressive weighting of the long baselines in the imaging process does not present any definitive evidence of multiplicity at present. In either case, IRS1 is an interesting target to follow up with higher-resolution interferometer observations to resolve possible multiplicity and to investigate the accretion and jet ejection process associated with intermediate-mass protostars.