HOPS 361-C's Jet Decelerating and Precessing through NGC 2071 IR

NGC 2071 was observed twice with HST's Wide Field Camera 3/Infrared (WFC3/IR). The first epoch is centered on NGC 2071 IR and the protostellar sources HOPS 361 (A)–(F). These images were taken on 2010 March 8 from HST General Observer (GO) proposal 11548 (Kounkel et al. 2016; Habel et al. 2021). The second epoch of images were taken between 2021 September 16 and November 3 as part of a joint HST/SOFIA GO proposal 16493. That program observed three fields to cover the entire HOPS 361-C jet, the HOPS 361-A outflow cavity, and the region surrounding these features.

The image frames were acquired with HST/WFC3 in star tracker-guided mode and using a standard 4-point dither. WFC3's field of view is 136'' × 123'' at NIR wavelengths. The images were captured using the following narrowband, NIR filters (and relevant spectral features): F126N (1.26 μm [Fe ii]), F128N (1.28 μm Paβ), F130N (1.30 μm continuum), F164N (1.64 μm [Fe ii]), F167N (1.67 μm continuum). Images are standard data products outputted by the calwf3 pipeline, retrieved from the Barbara A. Mikulski Archive for Space Telescopes (MAST) with suffix "FLT," and have details in Table 1.

We aligned images using the Python drizzlepac package with default parameters. The images in all filters were simultaneously aligned using the TweakReg routine version 1.4.7 to give a subpixel accuracy. We estimate errors in the image shifts to be at most 0.05 pixels or ∼0006 based on comparing centroids of stars in images at the two different epochs.

The Astrodrizzle algorithm (Hoffmann et al. 2021) version 3.3.0 was used with default parameters to resample the images to the same pixel size of 012825 × 012825 and remove defects (e.g., cosmic rays, bad pixels). The same routine was used to mosaic sets of second epoch narrowband images with the same filter. We then used the IPAC software Montage to shift, regrid, and update the World Coordinate System to the same frame size and orientation. This gave a final spatial (angular) resolution by 2D Gaussian FWHM of about 2 pixels (02565) for epoch (1) and about 3 pixels (038475) for epoch (2).

To convert e-/s to units of intensity, we used the flux density and bandwidth calibration values in the FITS file headers produced by the HST calibration pipeline (keywords PHOTFLAM and PHOTBW). The rms noise in the central frame for all images and within circular areas ranging from 1 to 15 arcsec² (used through this work) is 1.55 ×10⁻¹⁸ erg cm⁻² s⁻¹ per pixel as measured in regions devoid of stars and nebulosity.

For the precise coalignment needed to compute proper motions, the F160W frames were separately reprocessed through with drizzlepac, Astrodrizzle, and Montage as described with the corresponding second epoch F164N and F167N frames covering the closest overlapping area. Our images from both epochs centered on HOPS 361 are shown in the top left and top right of Figure 1. HOPS 361-C has an R.A. of 5^h47^m4631 and decl. of +0^d21^m4782 according to ALMA data in Tobin et al. (2020). This diverges from the coordinates found by Walther & Geballe (2019), which are based on Two Micron All Sky Survey K-band images, but we choose to use those from Tobin et al. (2020) because they were found with higher spatial resolution.

Zoom In Zoom Out Reset image size

For the second epoch images, the closest continuum filter was subtracted from each line image to construct narrowband images in [Fe ii] 1.26 μm and 1.64 μm as well as Paβ 1.28 μm. The appropriate continuum filters were F130N for F126N and F128N, and F167N for F164N. The continuum-subtracted narrowband images from the second epoch are shown in Figure 2.

Zoom In Zoom Out Reset image size

SOFIA observations taken on 2021 February 9–10 use the 4GREAT channel to observe the [O i] 63 μm line. The 4GREAT footprints overlap the HOPS 361 region, and [O i] emission traces the shocked material to give a jet's radial velocity and calibrate shock models, as shown with observations of other emission lines (e.g., Hartmann & Raymond 1984; Carr 1993; Hartigan et al. 2000; Graham et al. 2003) and theoretically for jets exhibiting bow shocks (e.g., Hartigan et al. 1987). The observations have a sensitivity of 0.04 K with 2 km s⁻¹ bins, resulting in a signal-to-noise ratio (S/N) of ∼6.

There is only one pointing that overlaps the base of the HOPS 361-C jet (see the circular region centered on the R.A. and decl. of 05^h47^m04723, 00^d21^m5784 in Figure 1). The footprint is shown with a diffraction limited beam size with a radius of 315. After removing the second-order baseline, the resulting [O i] 63 μm line spectrum is shown in Figure 3 with a Gaussian profile that is fit by least squares minimization. The fit gives a mean velocity of 3.3 ± 0.5 km s⁻¹ and a full width at half the maximum value of 15.8 ± 2 km s⁻¹. Comparing the mean velocity with the HOPS 361-C protostar's systemic velocity of 9.5 km s⁻¹, Cheng et al. (2022) yields an average radial velocity at the base of the jet of approximately −6.2 km s⁻¹ and a velocity spread with FWHM of about 16 km s⁻¹.

Zoom In Zoom Out Reset image size

We need spatially varying extinction measurements in order to interpret line ratios in terms of shock models and to estimate how deeply a protostar and its jet are embedded within the molecular cloud. We estimate extinction with the second epoch [Fe ii] 1.26 μm and [Fe ii] 1.64 μm spectral line images. These lines share a common upper state, so their intrinsic or unextinguished ratio can be computed from the wavelengths and Einstein A-coefficients for the transitions. The primary uncertainty is that the Einstein A-coefficients are difficult to determine theoretically.

We mitigate uncertainties in A-coefficients by adopting an observed zero-reddening value of 2.6 for the [Fe ii] line intensity ratio based on observations of bright HH objects in NGC 1333 (Rubinstein 2021). Next, we use the R_V = 5.5 model by Weingartner & Draine (2001) to find an optical depth from the ratio of the observed and theoretical [Fe ii] line intensity ratios. The resulting map of A_V for the HOPS 361 region is shown in Figure 5. An S/N of 2 is used to clip any signal in either HST line image below twice the noise of 1.55 × 10⁻¹⁸ erg cm⁻² s⁻¹ per pixel from Section 2.1.

Zoom In Zoom Out Reset image size

We use the regions we define in Section 3 to estimate the mean extinction to each knot, both of which are listed in Table 2. For uncertainties, we use the "error" image or uncertainty plane directly outputted by drizzlepac for each image. The uncertainty for each pixel consists of exposure times taking account of Poisson noise (e.g., shot noise), read noise, and other sources of uncertainty. We propagate through the steps of calculating line ratios and extinction. The uncertainties show the extinction is well estimated in the vicinity of HOPS 361-C and the arced jet where we detect proper motions, until the end of the jet where the extinction drops off.

The emission-line knots along the HOPS 361-C outflow are HH objects: ionic-line emission follows the passage of a dissociative and ionizing shock produced by supersonic variations in outflow speed or the outflow supersonically interacting with the ambient medium. At each knot, an outflow, initially launched at speed v_j and pre-shock hydrogen number density n₀, enters a shock at relative speed or shock speed v_S . It decelerates and compresses in a (Rankine–Hugoniot) jump to speed v_S /4 and density 4n₀, then to a lower speed and higher density determined by momentum conservation and radiative cooling.

Post-shock velocity, the vector sum of proper motion v_P and radial velocity v_R , differs from outflow velocity v_outflow by v_S . Post-shock ionic-line-emission intensities are sensitive to v_S and n₀, so v_P and v_R measurements enable us to determine outflow speed v_outflow ≈ v_S + v_P , mass density ρ₀ = μ n₀ (μ as the mean molecular weight), mass flux ρ₀ v = μ n₀ v, momentum flux ρ₀ v², and kinetic energy flux ρ₀ v³ lost by the outflow at every knot.

We determine v_S and n₀ from observed line emission with the 1D shock and photoionization code MAPPINGS V (Sutherland et al. 2018). For [Fe ii], we substitute the network of [Fe ii] radiative and collisional rate coefficients calculated by Tayal & Zatsarinny (2018), but otherwise use the default atomic–physical data from the CHIANTI database (Del Zanna et al. 2021). The typical protostellar shock speeds are small enough that little sputtering of dust is expected, so we assume the same gas-phase abundances as Watson et al. (2016). We generate a grid of model line-emission spectra with v_S = 20–62 km s⁻¹, and n₀ = 0.1–320,000 cm⁻³. At each grid point, we assume a magnetic field oriented perpendicular to the outflow and frozen into the medium at typical interstellar strength, $B\,=\,5\,\mu G\sqrt{{n}_{0}/100\,{\mathrm{cm}}^{-3}}$. Each calculation terminates at gas temperature T = 100 K. While MAPPINGS is 1D, HH objects are well-resolved structures observed close to edge-on, so the fits to the observed intensity can still estimate v_S and n₀.

Our model-grid results, derived parameters, and sensitivity to observations are in Figure 6. The intensity ratios are from our HST images. Blue, nearly vertical contour lines can be a ratio of H i Pa β 1.28 μm to either NIR [Fe ii] line. This line ratio is sensitive to v_S and relatively insensitive to n₀. Green, nearly horizontal contours show the ratio of mid-infrared [Si ii] 34.8 μm and [Fe ii] 26.0 μm emission lines. This ratio is sensitive to n₀ and insensitive to v_S . The northern half of the jet, imaged with the Spitzer Infrared Spectrograph at 15'' spatial resolution, includes the spectral lines [Fe ii] 26.0 μm and [Si ii] 34.8 μm (Melnick et al. 2008).

Zoom In Zoom Out Reset image size

Three bright regions in the Spitzer images (labeled P1, P2, and P3 in Figure 6 of Melnick et al. 2008) are chosen to measure the [Si ii] 34.8 μm to [Fe ii] 26.0 μm line ratio of ∼0.45. Peaks in the HST images in the same regions are used to measure the ratio of [Fe ii] 1.26 μm to Paβ, about −0.15. These two line ratios are placed on the MAPPINGS shock model in Figure 6 with a pink circle whose size approximates the scatter from differences of ratios at the different peaks. From the location of the pink dot in Figure 6, we estimate n₀ ≈ 3.2 × 10⁴ cm⁻³, and v_S = 40–50 km s⁻¹.

HST images may show morphology not evident in relatively lower-resolution Spitzer images. We use the [Fe ii] 1.26 μm to Paβ line ratio, computed from the HST images alone, to create an image that shows the shock speed derived from the MAPPINGS model grid. Taking advantage of the density insensitivity for the NIR [Fe ii]/H i Paβ ratio, we fit a third-order polynomial to our grid points to estimate the shock speed, which is accurate to a few kilometers per second:

$\begin{eqnarray}\begin{array}{rcl}{v}_{S}(R) & = & -12.617\,{{\mathrm{log}}_{10}R}^{3}+30.684\,{{\mathrm{log}}_{10}R}^{2}\\ & & -42.044\,{\mathrm{log}}_{10}R+33.222,\end{array}\end{eqnarray} \tag{ 1 }$

where R is the [Fe ii] 1.26 μm/Paβ 1.28 μm intensity ratio. We use this relation to generate the shock-speed image in the bottom panel of Figure 7 with each line image at an S/N above 1. We neglect correcting the [Fe ii] 1.26 μm and Paβ lines for extinction. The two lines are close in wavelength, so the effect of extinction is negligible.

Zoom In Zoom Out Reset image size

Unlike proper motions, shock speeds do not appear to depend strongly on distance along the outflow from HOPS 361-C. The shock-speed map in Figure 7 shows the path of knots according to our proper motions may extend farther than about 0.2 pc from HOPS 361-C, following bright trails of molecular hydrogen emission (Walther & Geballe 2019; Eislöffel 2000) and extended line emission in Figure 2. Making conclusions about neighboring frames is difficult without detecting proper motions and bright knots to the southwest of the HOPS 361 clump. There appears to be a strip of shocked material neighboring HOPS 361-E, but the source is unclear in the first epoch. According to our extinction map in Figure 5, this part of the HOPS 361 clump has a higher degree of extinction. Ambiguity remains as the trail may consist of a single arc or form a longer wave as the result of previous periodic ejections in a precessing jet.

Using shock speeds (estimated via Equation (1)) and proper motions, we estimate the total outflow speed through each knot. We then examine how outflow speed depends on distance from the HOPS 361-C source. A decrease in flow velocity as a function of distance from the source may be consistent with a jet slowed or decelerated by the host molecular cloud. Such a jet may inject its energy and momentum into its host molecular cloud and locally affect stellar feedback.

To estimate total flow speed, following Coffey et al. (2004), Bally et al. (2002), we sum the 3D velocities of post-shock gas, as measured by the proper motions and radial velocity added in quadrature, with the velocity jump through the shock from Section 4.2. We ignore radial velocity because the majority of the speed is due to the 100–350 km s⁻¹ tangential motions as opposed to the radial velocity of <5 km s⁻¹, so these velocities indicate that the jet is oriented close to the plane of the sky (for examples exploring protostellar jets, their resultant shocks, and inclination effects, see Hartigan et al. 2000; Graham et al. 2003; Jhan et al. 2022). According to the 1D shock model described in Section 4.2, half of the velocity spread of shocked material measured from [O i] in Figure 3 should be between the speed of turbulent gas measured from C¹⁸O (1 km s⁻¹; see Iwata et al. 1988; Stanke et al. 2022) and roughly a quarter of the shock speed (10 km s⁻¹). With a velocity spread of 6.2 km s⁻¹, our assumed radial velocity and low inclination are consistent with the model.

In Figure 8, we show quantities for each knot listed in Table 3 as a function of distance from the likely protostellar source, HOPS 361-C. In the top panel, blue triangles show tangential velocity directly found from the proper motions. Yellow x's show the velocity of the shocked gas, estimated from the [Fe ii]/Paβ line ratio. The green dots show the sum of these two quantities, which is an estimate for each knot's total outflow velocity. The red diamond shows the radial velocity for the jet from our SOFIA spectrum in Figure 3.

Zoom In Zoom Out Reset image size

The top panel of Figure 8 shows that the outflow velocity is decelerating as a function of distance from the source. Because proper motions dominate the speed, their uncertainty values are plotted by applying deviations of 2 pixel, from the angular resolution or PSF's FWHM of epoch (1), to the offsets listed in Table 3, which corresponds to a speed of about 44 km s⁻¹. The velocity drops from about 350 to about 100 km s⁻¹ over a distance of 0.2 pc, which may be because of collisions in a dense environment (e.g., Velázquez et al. 2013).

The amount of deceleration is relatively larger than that in more linear jets suspected to also impact into dense gas. For example, HH 1 is noted to have almost no deceleration (<50 km s⁻¹) over a distance of approximately 0.15 pc despite interactions within a dense clump (Raga et al. 2017; Castellanos-Ramírez et al. 2018). HH 2 is similarly noted to not have much or any deceleration over a distance of 0.02 pc after finding low dispersion in over 70 yr of data (Raga et al. 2016).

Jets have knot velocities typically aligned with the direction of their source (e.g., Raga et al. 1993; Heathcote et al. 1996; Schwartz & Greene 1999). In the middle panel of Figure 8, the blue points show the direction of motion for each knot relative to celestial north, based on the proper motions. Black squares show the position angle of the vector between each knot and HOPS 361-C. The uncertainties plotted for the proper motions apply the same worst case scenario using 2 pixel deviations in our measured offsets. The observed velocity vectors have position angles that are about 4° higher than what is predicted if the knots were shot out ballistically (i.e., a straight line on the sky) from HOPS 361-C. That said, our uncertainties cannot rule out ballistic motion on average for the majority of knots except the outlier at approximately 0.13 pc.

In case the discrepancy between the proper motion's direction and direction to source is real, the potential explanations are that the knots are launched by a different protostellar source than HOPS 361-C, a projection effect (e.g., more extreme radial velocities, another protostar launching the knot), or nonballistic motion. We searched for an alternative source but did not find a candidate among the known HOPS 361 protostars (A)–(E) and the nearby protostars HOPS 335 and HOPS 366. The deviation from ballistic motion for the outlier at 0.13 pc could be caused by deflection due to dense clumps of gas in the ambient molecular cloud material, which may change the direction of a knot (Raga & Canto 1995).

To match the positions and proper motions of our knots of ionized gas, we construct a model similar to those used to model wiggles observed in jets. For background and diagrams of precessing jets, please see Masciadri & Raga (2002), Anglada et al. (2007).

We adapt the analytical model from Raga et al. (2009) and invoke constant jet ejection velocity (Raga et al. 1993) but with an additional damping term (i.e., time-dependent exponential decay) to account for the observed decreasing speeds in Figure 8. The constant density jets with a decaying, time-dependent velocity profile are not without precedent in jet models (e.g., Kofman & Raga 1992; Velázquez et al. 2013), modeling observed radial velocities and proper motions for HH objects (e.g., HH 34 by Cabrit & Raga 2000; HH 80/81/80N by Masqué et al. 2015; HH 223 by López et al. 2015), and modeling knots of CO gas for a moving protostellar source (e.g., PV Cep by Goodman & Arce 2004).

We consider a binary comprised of two stars with masses m₁ and m₂ in a circular orbit of radius a_B . The jet is assumed to be emitted from m₁, and the binary orbit lies in the x, y plane. Due to the motion of m₁ about the center of mass of the binary system, the initial position from the jet source as a function of ejection time t_e is

$\begin{eqnarray}\begin{array}{ccc}x({t}_{e}) & = & {r}_{0}\cos ({\omega }_{B}{t}_{e});\\ y({t}_{e}) & = & {r}_{0}\sin ({\omega }_{B}{t}_{e});\end{array}\end{eqnarray} \tag{ 2 }$

where the binary system's mean angular motion is

$\begin{eqnarray*}&&{\omega }_{B}=\sqrt{\displaystyle \frac{G({m}_{1}+{m}_{2})}{{a}_{B}^{3}}}.\end{eqnarray*}$

Here, the distance from the binary's center of mass position to the jet source, m₁, is

$\begin{eqnarray}&&{r}_{0}=\displaystyle \frac{{m}_{2}}{{m}_{1}+{m}_{2}}{a}_{B}.\end{eqnarray} \tag{ 3 }$

The jet precesses if misaligned with the binary orbit's normal (Terquem et al. 1999). The precessing jet then has an initial velocity vector originating from m₁

$\begin{eqnarray}&&{{\boldsymbol{v}}}_{j}({t}_{e})={v}_{j}(\sin \beta \cos ({{\rm{\Omega }}}_{p}{t}_{e}-{\phi }_{0}),\sin \beta \sin ({{\rm{\Omega }}}_{p}{t}_{e}-{\phi }_{0}),\cos \beta ),\end{eqnarray} \tag{ 4 }$

where ϕ₀ is a phase with respect to the binary orbit, the angle β sets the opening angle of the jet with respect to the orbital plane, and v_j is the jet velocity. The jet is assumed to be ballistic, so the fluid parcels representing knots preserve their velocity at the launch. Here Ω_p is the precession rate of the jet. Taking into account the orbital motion of the binary, the initial velocity vector of jet material emitted at t_e is

$\begin{eqnarray}\begin{array}{rcl}{v}_{x}({t}_{e}) & = & -{r}_{0}{\omega }_{B}\sin ({\omega }_{B}{t}_{e})+{v}_{j}\sin \beta \cos ({{\rm{\Omega }}}_{p}{t}_{e}-{\phi }_{0});\\ {v}_{y}({t}_{e}) & = & {r}_{0}{\omega }_{B}\cos ({\omega }_{B}{t}_{e})+{v}_{j}\sin \beta \sin ({{\rm{\Omega }}}_{p}{t}_{e}-{\phi }_{0});\\ {v}_{z}({t}_{e}) & = & {v}_{j}\cos \beta .\end{array}\end{eqnarray} \tag{ 5 }$

We allow the velocity of emitted material to drop as a function of travel time with

$\begin{eqnarray}&&{\boldsymbol{v}}(t,{t}_{e})={\boldsymbol{v}}({t}_{e}){e}^{-\alpha (t-{t}_{e})},\end{eqnarray} \tag{ 6 }$

where α describes the decay rate. Here v (t_e ) is given in the previous equation. After a knot is emitted, its direction of motion does not change, so the jet propagates ballistically. However, the velocity of emitted material decelerates due to the α parameter. Without damping (in the limit of α → 0), the model reduces to that by Masciadri & Raga (2002) for orbital motion alone. What is more, if the precession rate is set to zero, or that for precession alone, then binary motion is neglected.

We integrate v (t_e ) in Equation (6) to find the position of ejected material at a later time t

$\begin{eqnarray}&&{\boldsymbol{x}}(t,{t}_{e})={{\boldsymbol{x}}}_{0}({t}_{e})+\displaystyle \frac{{\boldsymbol{v}}({t}_{e})}{\alpha }(1-{e}^{-\alpha (t-{t}_{e})});\end{eqnarray} \tag{ 7 }$

or specifically

$\begin{eqnarray}\begin{array}{rcl}x(t,{t}_{e}) & = & {v}_{x}({t}_{e})\frac{1}{\alpha }(1-{e}^{-\alpha (t-{t}_{e})})+{r}_{0}\cos ({\omega }_{B}{t}_{e});\\ y(t,{t}_{e}) & = & {v}_{y}({t}_{e})\frac{1}{\alpha }(1-{e}^{-\alpha (t-{t}_{e})})+{r}_{0}\sin ({\omega }_{B}{t}_{e});\\ z(t,{t}_{e}) & = & {v}_{z}({t}_{e})\frac{1}{\alpha }(1-{e}^{-\alpha (t-{t}_{e})}).\end{array}\end{eqnarray} \tag{ 8 }$

Here the ejection time is t_e , and the present time of observation is t. By setting the present time t = 0, the angle ϕ₀ sets the jet orientation at the present time. We can use Equation (8) to compute the position of ejected material as a function of ejection time and the velocity of this same material using Equation (6). This gives the position and velocity of ejected material in a coordinate system associated with the binary star.

We rotate along the x-axis to correct for the inclination i_B of the binary orbit normal vector. Here i_B = 0 corresponds to an edge-on binary orbit plane. Then we rotate the resulting coordinate system along the line of sight by angle ξ to correct for the position angle of the binary orbit's normal on the sky.

$\begin{eqnarray}&&\left(\begin{array}{c}{x}_{s}\\ {y}_{s}\\ {z}_{s}\end{array}\right)=\left(\begin{array}{ccc}\cos \xi & 0 & -\sin \xi \\ 0 & 1 & 0\\ \sin \xi & 0 & \cos \xi \end{array}\right)\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & \cos {i}_{B} & \sin {i}_{B}\\ 0 & -\sin {i}_{B} & \cos {i}_{B}\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right).\end{eqnarray} \tag{ 9 }$

Coordinates on the sky are x_s , z_s , and y_s , with y_s increasing away from the viewer along the line of sight. Here z_s is positive to the north, and x_s is positive to the west. This is a right-hand coordinate system with origin at the location of the source, HOPS 361-C. The rotations given in Equation (9) are also used to transform the velocity vectors. This gives v_xs , v_zs for motions in the plane of the sky (tangential motions), and v_r = v_ys corresponding to motion along the line of sight (the radial component). The resulting model has the following free parameters: v_j , Ω_p , α, β, i_B , ξ, ϕ₀, m₁, m₂, a_B .

We explore the parameter space by hand to see what ranges of parameters match the knot locations, the decrease in outflow velocity as a function of distance from the source protostar and the radial velocity of the first knot. We do not attempt to fit the semimajor axis (a_B ) or masses of the orbit's constituents (m₁, m₂) because these parameters only change the jet trajectory on scales smaller than observed, but we still allow the binary to change the jet's behavior. Lower inclinations relative to the plane of the sky (0°–5°; see Section 5.1) narrowed our search through the parameter space and better matched positions on the sky, proper motions, and opening angle (Figure 4). The first knot's low radial velocity (Section 2.2) and our measured proper motions were better matched with low initial phase angles (10 < ϕ₀ < 20).

We could not find a model that results in the observed deceleration (see Section 5.1) via projection effects alone. Models with an extreme opening angle, wider than what we measure in Figure 4 (33°), can exhibit a drop in proper motion as a function of distance from the source protostar. However, if our observed deceleration was caused by knots moving away from the observer, then the model knot positions predicted on the sky tend to be closer to the source than observed. We find that precessing jet models only match knot locations and velocities when the ejected knots rapidly decelerate, as described with our exponential decay parameter α. The models that match the positions and proper motions have half opening angle of β ≈ 15° and require moderate α values within about 20% of that given in the Table 4. Such a model is shown in Figures 8 and 9 and has parameters listed in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Precessing Jet Model Parameters

Parameter	Symbol	Fiducial Values
Jet ejection velocity	vj	−325 km s−1
Precession period	2π/∣Ωp ∣	2066 yr
Jet decay rate	α	1560 Myr−1
Half jet opening angle	β	145
Binary orbit inclination	iB	0°
Position angle of precession axis	ξ	235 + 180°
Initial phase	ϕ0	1275
Mass of jet source	m1	1 M⊙
Mass of binary companion	m2	1 M⊙
Binary semimajor axis	aB	20 au

Download table as:  ASCIITypeset image

Figure 9 shows knot positions on the sky from Table 2 as blue ellipses in the left panel. The black curve represents our model jet trajectory from HOPS 361-C. The plot axes x_s and z_s give coordinates on the sky in parsecs with origin at the jet source, which is shown with a fuchsia star. Celestial north is aligned with z_s , so it is facing up on this panel. A knot location can be used to estimate the time a knot was ejected from the source. The ejection times, estimated from this model, are shown on the top axis of Figure 8.

The molecular hydrogen knots from Walther & Geballe (2019) are shown by an approximate dashed curve, and the low-velocity CO cavity from Cheng et al. (2022) is shown with gray arrows. The helix of molecular hydrogen may be accelerated by the jet that we detect but cannot presently be analyzed in detail because of imaging differences (e.g., other protostars could heat gas, differences in brightness and epochs). We caution that past observations strictly correspond to any of the densest and warmest parts of the gas. The low-velocity CO gas encapsulates high-velocity CO gas and a jet of radio continuum emission from Trinidad et al. (2009), Carrasco-González et al. (2012) aligning with the molecular hydrogen knots and lying on top of the track of our detected [Fe ii] knots.

The right panel in Figure 9 shows the y_s coordinate of the model jet trajectory as a function of distance on the sky from the source (d_sky) as in Figure 8. Along the line of sight, the increasing y_s corresponds to greater depths into the cloud. With our model, the jet is moving toward then away from the observer as a function of d_sky.

HOPS 361-C was only recently resolved as a close binary system in VLA 9 mm continuum images with a model binary separation of about 01 or 43 au (Cheng et al. 2022). Following Terquem et al. (1999), the binary system's orbital period is related to the jet's precessional period, assuming a rigid disk of gas within the binary system that precesses, uniform disk surface density, and Keplerian orbits,

$\begin{eqnarray}&&{\tau }_{\mathrm{orbit}}={\tau }_{\mathrm{prec}}\displaystyle \frac{15}{32}\displaystyle \frac{q}{{(1+q)}^{1/2}}{\sigma }^{3/2}\cos (\beta ),\end{eqnarray} \tag{ 10 }$

where τ_orbit is the binary orbital period, q = m₂/m₁ is the ratio of the secondary mass divided by the primary, which is assumed to be the jet source. Here σ is the ratio of m₁'s disk radius to the binary orbital semimajor axis. Typically, σ is taken to be 1/3, assuming tidal truncation; though it could range from 1/4 to 1/2.

We found precession models match positions and proper-motion data with periods of 2π/∣Ω_p ∣ ranging from about 1500 to 2500 yr. With a 2000 yr precession period, a binary star comprised of two 1 M_⊙ stars (q = 1) has an orbital period of 130 yr and a semimajor axis of 33 au, which is approximately consistent with the separation inferred from the VLA 9 mm images. However, kinematic modeling of molecular line emission gave an inner disk radius that is smaller, 10–25 au (Cheng et al. 2022). If this inner radius encompasses the binary, then the binary semimajor axis would necessarily have to be smaller than the inner disk radius.

We note the precession model shown in Figures 8 and 9 is not sensitive to the binary separation or the mass of its constituents; though the small amplitude wiggles in the model curves are due to the binary orbital motion. With a larger binary separation, these wiggles would have lower amplitude and longer wavelength.

Via kinematic modeling, Cheng et al. (2022) estimate a binary mass for HOPS 361-C of m₁ + m₂ ∼ 1.5M_⊙; though modeling the spectral energy distribution gave a protostellar mass of ∼2M_⊙ neglecting binarity. These estimates for the protostellar masses are similar to our assumed values of m₁ = m₂, and m₁ + m₂ = 2M_⊙. Based on an estimate for the location of the kinematic center, Cheng et al. (2022) estimate a mass ratio m₁/m₂ ≈ 1.3 −2.5, which is within a factor of 2 of our assumed value of m₁/m₂ = 1. Using our current measurements and uncertainties, the jet precession model is not inconsistent with what is known about the HOPS 361-C binary system.

The direction of precession is predicted to be in the opposite direction as the rotation of the disk from which the jet originates (Ω_p is negative in Equation (1) by Terquem et al. 1999). Since the precession rate is computed by averaging the tidal force over the binary's orbit, it is independent of the direction of rotation in the binary orbit. We take Ω_p < 0 so that position angle and inclination are given with respect to the binary orbit normal. We assume the direction of rotation for the disk driving the jet is similar to the direction of rotation of the binary orbit (for which we chose β < π/2). The position angle of the precession axis, ξ, is also the position angle of the binary orbit normal. This axis points to the south, consistent with the direction of rotation seen in the disk resolved by Cheng et al. (2022).

The CO J = 2−1 data for HOPS 361-C also show a high-velocity, bipolar jet with a position angle of 22°–32°. This differs from the radio jet axis seen at 9 mm that has a position angle of ∼15° (Cheng et al. 2022). Near the source, our precession model has the jet oriented near its maximum opening angle, giving a jet position angle of about 10° and approximately consistent with the radio jet axis.

The jet opening angle β is sensitive to the solid angle extended by the arc with respect to our origin, the protostellar source. The model is less sensitive to the binary inclination i_B , but we found models that matched both proper motions and knot positions with i_B < 25°. The radial velocity from the [O i] spectrum restricted i_B to a few degrees. Cheng et al. (2022) fit kinematic models to their ALMA molecular line data for the HOPS 361-C disk–jet system, and they estimated a disk inclination of 63°–77° and a position angle of around 15° for the radio jet. If we assume their disk inclination is that of a circumbinary disk in a plane containing the binary, this inclination corresponds to i_B ∼ 90° − 70° = 20° for the inclination of the binary orbit normal vector using our variable, which has zero inclination for an edge-on binary orbit. Our precession model is consistent with the estimates for the disk inclination by Cheng et al. (2022), assuming their model pertains to a circumbinary disk.

The disk resolved by Cheng et al. (2022) has a blueshifted component on the western side and redshifted component on the eastern side, which suggests that the binary orbit's angular momentum vector is pointed to the south. To compare, if the binary orbit and the disk that drives the jet have the same direction of rotation (opening angle β < π/2), then the jet's arc is part of a counter-jet that propagates in the direction opposite to the binary orbit's normal. The position of the counter-jet can be predicted using Equation (4) with v_j < 0, giving the expected S-shaped symmetry about the origin between jet and counter-jet. This is why we have chosen negative v_j for our model in Table 4.

The direction of precession in our model is consistent with a disk-driven jet having rotation direction similar to the presumed circumbinary disk observed in various molecular lines by Cheng et al. (2022). We concur with Cheng et al. (2022), who concluded that jet precession and associated interactions between jet and environmental material at different times are the most likely scenario accounting for the observed differences in the jet position angles seen with different tracers.

We achieve a quasi-3D positioning of the knots by comparing foreground extinction to each knot in Figure 8 (or see Table 2) with the extinction to material at the back of the molecular cloud. Their ratio gives an approximate value for how deeply embedded each knot is within the cloud. But the relative distances may change due to radial distance, column density, number density, or opacity. Therefore, they may not directly correspond to the radial positions predicted by our precession model (the right panel of Figure 9).

The extinction to the back of the cloud can be determined using molecular gas, like optically thin CO isotopologues, that traces relatively slow motions (e.g., Dickman 1975, 1978; Goldsmith et al. 1992; Schwartz et al. 1983; Gutermuth et al. 2008; Stojimirović et al. 2008). We use the correlation between extinction and C¹⁸O (J = 1–0) molecular line emission (Alves et al. 1999):

$\begin{eqnarray}&&{A}_{V}=\displaystyle \frac{I({{\rm{C}}}^{18}{\rm{O}})+0.4\pm 0.1}{0.21\pm 0.01}.\end{eqnarray} \tag{ 11 }$

The region of the integrated C¹⁸ O map from Iwata et al. (1988) that overlaps our NIR images gives a brightness temperature of 2.52–2.8 K km s⁻¹, which corresponds to an A_V of 13.9–15.2 mag using Equation (11). We only find a single value for the region because the radio observations have a nearly 1' spatial resolution, nearly the scale of our jet. We could use C¹⁸ O data in Schwartz et al. (1983), but their column density maps have fewer contour lines to use for measurements.

We present the fraction of each knot's A_V relative to total A_V through the cloud in the bottom panel of Figure 8. Translating the contours from Iwata et al. (1988) to an A_V value with the correlation from Alves et al. (1999) introduces an additional uncertainty of approximately ±1.6 mag to the uncertainties in the extinction, and we have used this uncertainty to estimate the error bars for the points plotted in the bottom panel of Figure 8.

The fractional depths in the bottom panel of Figure 8 indicate knots oscillate between being more and less embedded, potentially moving closer and farther from the viewer. We take caution making firm constraints with our precession model, since the cloud is unlikely to be perfectly uniform (e.g., Schwartz et al. 1983). But precession may impact jet morphology and dynamics more than variations in a jet's gas density (e.g., Castellanos-Ramírez et al. 2018). Quantitatively, the knots are at a fractional depth of 1/5–4/5 into the molecular cloud, where they may disrupt the HOPS 361 region.

We search for effects that connect jet properties relevant to stellar feedback, our precession model values, and HOPS 361 cloud clump properties. Assuming a 1D, incompressible fluid flow and mass conservation across each knot (ignoring changes in viscosity, plasma density, and mass flow rate) as in the MAPPINGS model, the jet has a mass outflow rate or jet output ($\dot{M}$) of

$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{{\rm{\Delta }}M}{{\rm{\Delta }}t}={\rho }_{j}{{Av}}_{\mathrm{outflow}},\end{eqnarray} \tag{ 12 }$

and the momentum flux ($\dot{P}$) or ram pressure through each knot for constant $\dot{M}$ as

$\begin{eqnarray}&&{\rm{\Delta }}\dot{P}=\displaystyle \frac{\dot{M}}{A}{\rm{\Delta }}v=\displaystyle \frac{\dot{M}}{A}({v}_{\mathrm{outflow}}-{v}_{P})=\displaystyle \frac{\dot{M}}{A}{v}_{S},\end{eqnarray} \tag{ 13 }$

and a flux of kinetic energy ($\dot{K}$) or energy density dissipated in the shock in each knot with constant $\dot{M}$ of

$\begin{eqnarray}&&{\rm{\Delta }}\dot{K}=\displaystyle \frac{1}{2}\,\displaystyle \frac{\dot{M}}{A}{\rm{\Delta }}({v}^{2})=\displaystyle \frac{1}{2}\,\displaystyle \frac{\dot{M}}{A}({{v}_{\mathrm{outflow}}}^{2}-{{v}_{P}}^{2}).\end{eqnarray} \tag{ 14 }$

Here ρ_jet is the jet's mass density, v_outflow is total flow speed through each knot, v_P is knot proper motion, v_S is shock speed, and A is the jet's cross-sectional area. We compute $A=\pi {{R}_{\mathrm{knot}}}^{2}$, taking R_knot to be the semiminor axis from Table 2 for an ellipse encircling each knot. As in Section 4.2, we estimate v_outflow from the sum of the knot proper motion and shock speed. Assuming the majority of the jet's mass consists of neutral atomic hydrogen, we find the mass density by multiplying the jet's number density (n₀ ≈ 3.2 × 10⁴ cm⁻³, or see Section 4.2) by the mass of a hydrogen atom and using a mean molecular weight of 1. We estimate the total mass density in the jet is ρ_jet = 5.34 × 10⁻²⁰ g cm⁻³ = 789 M_⊙ pc⁻³. We list each knot's $\dot{M}$, $\dot{P}$, and $\dot{K}$ multiplied by area A for ease of comparison in Table 5.

Table 5. Derived Feedback Properties for Each Knot

Knot Identifier	A	$\dot{M}$	Δ$\dot{P}$ × A	${\rm{\Delta }}\dot{K}\times A$
...	pc2	M⊙ yr−1	M⊙ yr−1 km s−1	L⊙
1	7.69 × 10−6	2.12 × 10−6	1.10 × 10−4	5.69
2	6.62 × 10−7	1.32 × 10−7	6.61 × 10−6	0.240
3	1.29 × 10−5	2.93 × 10−6	1.38 × 10−4	5.85
4	2.77 × 10−6	4.29 × 10−7	2.19 × 10−5	0.598
5	3.28 × 10−6	4.62 × 10−7	2.27 × 10−5	0.558
6	6.90 × 10−6	9.72 × 10−7	5.08 × 10−5	1.24
7	5.78 × 10−6	7.56 × 10−7	3.97 × 10−5	0.885
8	4.29 × 10−6	4.62 × 10−7	2.21 × 10−5	0.397
9	1.16 × 10−5	1.33 × 10−6	6.57 × 10−5	1.27
10	2.17 × 10−5	1.79 × 10−6	8.91 × 10−5	1.13
11	1.54 × 10−5	1.24 × 10−6	5.99 × 10−5	0.740

Note. Jet properties used to evaluate how the HOPS 361-C jet affects its host cloud. The cross-sectional area derives from squaring the smaller, minor axis of the ellipse centered on each knot and using values from Table 2. The flow speeds needed for other columns are listed in Table 3.

Download table as:  ASCIITypeset image

According to Watson et al. (2016), Sperling et al. (2021), the jet output efficiency is the ratio of the mean mass outflow rate through the jet to the accretion rate. Mass accretion rate can be estimated for class 0 and class I protostars using bolometric luminosity (L_bol) and stellar properties (mass, M_*, and luminosity, L_*) by assuming that the bolometric luminosity is equal to the total luminosity for the binary system as well as the accretion luminosity; then

$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}=\eta \displaystyle \frac{{L}_{\mathrm{bol}}{R}_{* }}{{{GM}}_{* }}\end{eqnarray} \tag{ 15 }$

(Enoch et al. 2009; Evans et al. 2009), where η is an efficiency factor often assumed to be on the order of 1 (Meyer et al. 1997; Calvet & Gullbring 1998; Muzerolle et al. 2003; Fischer et al. 2017). Using the methods from Furlan et al. (2016) to estimate bolometric luminosity, Cheng et al. (2022) estimated a mass accretion rate of 9 × 10⁻⁶ M_⊙ yr⁻¹ for HOPS 361-C. Using, the mean $\dot{M}$ from Table 5 of 1.15 ×10⁻⁶ M_⊙ yr⁻¹ and the above accretion rate estimate, we find an efficiency of 0.128. Efficiency is related to the location of the jet launching site, with larger efficiencies corresponding to launch at a smaller radius (Watson et al. 2016; Sperling et al. 2021). Based on Watson et al. (2016), Sperling et al. (2021), the values around 0.1 cannot rule out any jet launching models. If the bolometric luminosity is split between the components of the binary system, then the efficiency may increase, and the jet launching models that place the footpoint closer to the protostar may better explain why our value exceeds 0.1.

We determine the momentum injection rate through the jet imparted into the surrounding cloud, compare with escape velocity, and find the amount of mass the jet ejects from the HOPS 361 clump (as in Matzner 2007; Arce et al. 2010; Li et al. 2020). We sum $\dot{P}\times A$ values for each knot listed in Table 5. Because the interaction with the surroundings may not be perfectly efficient or uniform, then the total force is ≤6.26 ×10⁻⁴ M_⊙ yr⁻¹ km s⁻¹. Taking the clump as a uniform, nonrotating, nonmagnetic, spherical mass of M_c = 88.9 M_⊙, at a radius of R_c = 0.2 pc (using the Herschel map of Orion B described in Stutz & Gould 2016; Stutz p. com.; and Furlan et al. 2016; and the correlation between N(H) and extinction found by Pillitteri et al. 2013), and then the clump's escape velocity (v_esc) is 2 km s⁻¹. The ratio of the sum of momenta from the jet's knots to this escape velocity gives $\dot{M}$ _esc for material escaping the cloud equal to 3 ×10⁻⁴ M_⊙ yr⁻¹, assuming constant escape velocity and momentum output. Over a typical protostellar lifetime or a freefall timescale of 100,000 yr, then the jet would launch approximately 30 M_⊙ out of the system, or about a third of the clump's mass. Therefore, a few precessing jets from intermediate-mass protostars may have enough momenta to eject the majority of the clump's mass at this scale.

To find the effect of the jet's kinetic energy, we sum power ($\dot{K}\times A$) emitted by each knot in Table 5 to get a total of 18.6 L_⊙. The gravitational potential energy of an extended mass or cloud, M_c , with radius, R_c , is generally

$\begin{eqnarray}&&U=-C\displaystyle \frac{G{{M}_{c}}^{2}}{{R}_{c}},\end{eqnarray} \tag{ 16 }$

where C is a constant dependent on the cloud's geometry and uniformity. Assuming the mass and radius to compute our escape velocity, Equation (16) gives 1.21 × 10⁴⁵ erg. If the cloud is virialized and only has a total of half of this energy, then the ratio of half of U to $\dot{K}\times A$ is a 480 yr timescale. This timescale may suggest a precessing jet contributes a significant portion of a cloud clump's gravitational energy budget during a protostellar lifetime. We caution this is the fastest possible case with our present assumptions (perfectly efficient energy transfer, uniform cloud, constant output).

The jet's power may contribute to turbulence within the molecular cloud. The dissipation rate due to turbulence at a velocity of v_turb in the molecular cloud is

$\begin{eqnarray}&&{L}_{\mathrm{turb}}\sim \displaystyle \frac{{M}_{c}{v}_{\mathrm{turb}}^{3}}{{l}_{\mathrm{eddy}}},\end{eqnarray} \tag{ 17 }$

where l_eddy is the size of the largest eddies within the cloud (e.g., Stone et al. 1998; Elmegreen & Scalo 2004; Quillen et al. 2005). Taking the NGC 2071 cloud mass M_c ∼ 2000 M_⊙and l_eddy ∼ 1 pc (from ¹²CO and ¹³CO maps; see Stojimirović et al. 2008), v_turb ∼ 1.9 km s⁻¹ (based on C¹⁸O line widths by Iwata et al. 1988; and Stanke et al. 2022; also a velocity FWHM of 1 km s⁻¹ from an NH₃ line map using the Green Bank Telescope, by priv. comm. with J. Di Francesco), and then we estimate a turbulent dissipation luminosity of about 2.3 L_⊙. Because the knots deliver an average power of about 1.69 L_⊙ into the cloud, the energy dissipation rate via turbulence in the cloud is similar to the average power dissipated locally by the jet. The jet itself may have enough energy to locally drive the molecular cloud turbulence.

The HOPS 361-C jet may have enough momentum and energy to disrupt the HOPS 361 protostar-forming clump, but that depends on how efficiently the jet outputs material over time. The jet dissipates over a distance of ∼0.2 pc (Figures 8 and 9). A knot launched at 325 km s⁻¹ traveling 0.2 pc is a dynamical timescale of 600 yr, which is similar to our damping timescale, the reciprocal of α, of 640 yr. The damping or dynamical timescale for knots can occur about 3 times within our precession timescale of 2000 yr, so the jet may extend farther as seen in Figure 2. However, the knots outside this 0.2 pc distance may be due to other outflows or may be associated with lower-velocity shocks <20 km s⁻¹, so future observations are needed to determine their source. Our study of the arced jet associated with HOPS 361-C suggests the precessing jets with wide opening angles are locally damped within 0.2 pc of their driving source.

What condition allows subsequently emitted clumps in a jet to dissipate or remain within the cavity opened by previously emitted clumps? For small opening angles, the clumps in protostellar outflows lie on a similar path as the jet itself, and these jets open a cavity within a molecular cloud (Quillen et al. 2005; Cunningham et al. 2009a; Frank et al. 2014; Fendt & Yardimci 2022). Long jets with small wiggles can puncture their host molecular cloud (e.g., simulations by Velázquez et al. 2013) and extend more than 1 pc from their source (Bally 1982, 2016). In contrast, instead of passing through the molecular cloud, the HOPS 361-C jet seems to entirely dissipate locally within it.

We estimate the speed that a cavity is opened by a jet and compare it to its advance speed. The radial distance and radial velocity from the jet's central axis of symmetry are

$\begin{eqnarray}&&r(t)={v}_{j}\sin (\beta )t;\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray}&&\dot{r}(t)={v}_{j}\sin (\beta )\sim \mathrm{constant};\end{eqnarray} \tag{ 19 }$

where v_j is the jet speed, and β is the half opening angle. A jet cavity that would not fill in at the speed of sound, c_s , in the ambient molecular cloud would satisfy

$\begin{eqnarray}&&\dot{r}\lesssim {c}_{s}.\end{eqnarray} \tag{ 20 }$

This corresponds to a jet cone that moves perpendicular to its central axis slower than the sound speed. Combining this equation with Equation (19), a jet cavity that does not fill in would satisfy

$\begin{eqnarray}&&\sin (\beta )\lesssim \displaystyle \frac{{c}_{s}}{{v}_{j}}.\end{eqnarray} \tag{ 21 }$

Using a jet velocity of v_j =400 km s⁻¹ and a sound speed for the cloud's interstellar medium in the of c_s ∼ 1 km s⁻¹, we estimate the half opening angle β ≲ 1° for a jet that can maintain a cavity. In this scenario, any jet traveling at a few 100 km s⁻¹ within an opening angle of a degree requires motions above the local sound speed to fill in the associated cavity. If the opening angle is larger than a degree, then the cavity would fill in. Opening a cavity requires energy. A jet with a wide opening angle may then expend more energy to propagate. This suggests wider jets would dissipate more rapidly, and long jets with large opening angles should be rare.

This scenario ignores turbulence and magnetic fields. But note that the supersonic, turbulent speed of gas in the interstellar medium is approximately Mach 3, where the Mach number is the ratio of motion in the cloud to the local sound speed (Ballesteros-Paredes et al. 2007). To take into account turbulence in the molecular cloud, we multiply the right-hand side of Equation (21) by a factor of 3, and this would relax the restriction on the opening angle for a jet that can puncture a molecular cloud.

2D magnetohydrodynamic simulations have investigated how a precessing jet nozzle affects the propagation of a high-speed jet (Fendt & Yardimci 2022). The jets with wide opening angles (>20°) fail to propagate to large distances as they are dissipated locally (see their Figure 6) but suggest the condition allowing a jet to propagate is less restrictive than that by Equation (21). Future simulations could relate protostellar jet observations (e.g., knot proper motions) to other cloud, jet, and binary properties and elucidate connections between the jet opening angle and jet dissipation. Simultaneously observing molecular, atomic, and ionized components of the HOPS 361-C jet, including potential follow-up for proper motions, may also confirm H₂ knots accelerated and carved out by HOPS 361-C.