Proper motions of young stellar outflows in the mid-infrared with Spitzer II HH 377/Cep E

We have used multiple mid-infrared observations at 4.5 μm obtained with the infrared array camera, of the compact ( ∼ 1.4 ′ ) ?> young stellar bipolar outflow Cep E to measure the proper motion of its brightest condensations. The images span a period of ∼ 6 ?> yr and have been reprocessed to achieve a higher angular resolution ( ∼ 0.8 ″ ) ?> than their normal beam ( ∼ 2 ″ ) ?> . We found that for a distance of 730 pc, the tangential velocities of the north and south outflow lobes are 62 ± 29 ?> and 94 ± 26 km s − 1 ?> respectively, and moving away from the central source roughly along the major axis of the flow. A simple 3D hydrodynamical simulation of the H2 gas in a precessing outflow supports this idea. Observations and models confirm that the molecular hydrogen gas, traced by the pure rotational transitions, moves at highly supersonic velocities without being dissociated. This suggests either a very efficient mechanism to reform H2 molecules along these shocks or the presence of some other mechanism (e.g. strong magnetic field) that shields the H2 gas.

94 26 km s 1 respectively, and moving away from the central source roughly along the major axis of the flow. A simple 3D hydrodynamical simulation of the H 2 gas in a precessing outflow supports this idea. Observations and models confirm that the molecular hydrogen gas, traced by the pure rotational transitions, moves at highly supersonic velocities without being dissociated. This suggests either a very efficient mechanism to reform H 2 molecules along these shocks or the presence of some other mechanism (e.g. strong magnetic field) that shields the H 2 gas.

Introduction
The measurement of proper motions for Herbig-Haro (HH) outflows and stellar jets has a long tradition and has played a fundamental role in our understanding of early phases of evolution of low mass stars, including their accretion rates, mass loss and disk dissipation (see e.g. McKee andOstriker 2007, Bally 2009) The original work by Herbig and Jones (1981) on the first and brightest HH objects 1, 2 and 3, using photographic plates over a 34 yr period, set up the framework of the mass loss process in proto-stars and their interaction with the surrounding medium. The 'knots' of HH 1 and 2 were found to have tangential velocities ranging from 100 to 350 − km s 1 , in disagreement with the spectroscopic measurements that indicated that their emission was due to shocks of at most 100 − km s 1 (see e.g. Raymond 1979). Modern observations and models have shown that HH 1 and 2 are the leading working surfaces or 'bowshocks' of a highly collimated jet/counter-jet system that arises from a deeply embedded protostellar source of class 0. And that a time dependent ejection can account for both larger proper motion and relative smaller shock velocities between their knots (see e.g. Raga et al 2011a for a review). At optical wavelengths, using narrow band images at some of the brightest collisionally excited emission lines (e.g. Hα and/or [SII]), it has been possible to measure proper motions for ∼50 HH flows, within a distance of roughly one kiloparsec, both from the ground and space (see e.g. Bally et al 2007, Garatti & Eislöffel 2009) and over relatively short periods of time ⩽ ( 10-20 yr). In the case of objects like HH 1/2, 34 and 46/47, such measurements have led to some spectacular time sequences of the outflows (see e.g. Hartigan et al 2011). Proper motions are still a fundamental tool to find and correlate outflows across the sky and over distance scales of parsecs (see e.g. Reipurth et al 2013). This method is particularly important when the outflow driving source has not been clearly identified. Radio observations have also been very successful in measuring proper motions of outflows for low and high mass protostars, taking advantage of the gas flow thermal (free-free) emission and high angular resolution measurements obtained by interferometric observations (see e.g. Rodriguez 2011).
At infrared wavelengths, however, it has been more difficult to measure proper motions because of the lack of large format arrays. It has not been until recently that near/mid-IR arrays have had a wide enough FOV to include a reasonable number of reference stars and be able to cross-correlate multiple epoch observations to derive the proper motions. The large proper motions observed in the HH 1/2 system in the atomic ionic gas (Herbig and Jones 1981) have been measured in the NIR as well, using the molecular hydrogen shock excited emission vibrational transition (v = 1 − 0 S(1)) at 2.121 μm . It is now possible, using similar tracers, to measure the motions of multiple H 2 features in Cha II (Caratti o Garatti et al 2009) or ρ Oph clouds (Zhang et al 2013), with a similar range of transversal velocities (from 30 to 120 − km s 1 ). In the mid-IR, thanks to the stability and longevity of the infrared array camera (IRAC; Fazio et al 2004) on board the spitzer space telescope, it has been possible to measure the proper motions of several outflows in NGC1333 at an angular resolution of ″ 2 (Raga et al 2013). IRAC channel 2 at 4.5 μm is an excellent tracer of the H 2 rotational emission, since three of the brightest lines, 0-0 S(11) 4.18, 0-0 S(10) 4.40 and 0-0 S(9) 4.18 μm fall within its passband (Noriega-Crespo et al 2004a, Noriega-Crespo et al 2004b, Ybarra and Lada 2009, Raga et al 2011a. Using images obtained at 4.5 μm over a period of ∼7 yr, Raga et al (2013) obtained tangential velocities ranging from ∼10 to 100 − km s 1 (for a 220 pc distance) for eight outflows in NGC 1333 cloud. For the bright HH 7-11 system, that lies at the center of the cloud, the H 2 tangential velocities of ∼10-15 km s −1 are very close to those obtained at optical wavelengths using atomic gas line tracers, such as [SII] 6717/31 Å or Hα Jones 1983, Noriega-Crespo andGarnavich 2001).
In summary, proper motions, and the corresponding tangential velocities, are essential for determining the dynamics of the outflows, plus their momentum and energy transfer into the surrounding interstellar medium (see e.g. Padoan et al 2009, Plunkett et al 2013, Quillen et al 2005. The striking morphological similarity between the atomic/ionic gas emission (obtained from optical or near-IR observations) and that of the molecular hydrogen (obtained either from near/mid-IR observations), suggests that the kinematics of the protostellar outflows allow this relatively fragile molecule (H 2 ) either to survive the shocks or to regenerate itself rapidly in the dense postshock regions (see e.g. Le Bourlot et al 2002, Panoglou et al 2012. These issues can be partially addressed with studies of proper motions of the H 2 emission from stellar outflows. In this study we determine the proper motions of the deeply embedded and compact molecular outflow Cep E, driven by a intermediate mass class 0 protostar , Lefloch et al 1996, Ladd and Hodapp 1997, Noriega-Crespo et al 1998, Hatchell et al 1999, Moro-Martín et al 2001, Noriega-Crespo et al 2004b. Cep E is considered an excellent prototype of its kind, and therefore, has prompted many recent observations at millimeter and sub-millimeter wavelengths to study its H 2 O and CO molecular content (Lefloch et  At a distance of ∼ 730 pc to Cep E, measuring proper motions with a time interval of ∼6 yr (covered by the available IRAC images) is a considerable challenge, since velocities of ∼100 − km s 1 would correspond to shifts of only ∼ ″ 0.17 . In order to achieve as high a resolution as possible, we have employed a high angular resolution enhancement of the IRAC images, reaching a resolution of ″ 0.6 -″ 0.8 (see Velusamy et al 2007,Velusamy et al 2008,Velusamy et al 2014. Such an enhancement has recently been successfully applied to IRAC images of Cep E (Velusamy et al 2011). Finally, given that observationally there is a tremendous morphological similarity between the mid-IR and NIR emission (Noriega-Crespo et al 2004b), we expand the study to include some ground based H 2 2.12 μm NIR data that allows us to extend the time baseline of the observations to ∼16 yr (table 1).
The paper is organized as follows. The observations and their high angular resolution reprocessing are described in section 2. The determination of the proper motions is described in section 3. Finally, section 4 presents a summary of the results.

Observations and high angular resolution reprocessing
The Cep E outflow was one of the first young stellar outflows observed with the Spitzer space telescope as part of their early release observations (PID 1063, P.I. Noriega-Crespo) because of its strong brightness at mid-IR wavelengths. Its emission in the mid-IR (5-17 μm) is due mostly to bright H 2 rotational lines clearly detected already by the infrared space observatory with the infrared camera using its circular variable mode (see e.g. Boulanger et al 2005) that provided a low spectral ∼ R ( 45) and angular (FWHM ∼ ″ 6 ) resolution spectral map of the region (Moro-Martín et al 2001). The outflow was later observed by two ambitious programs, one to map the Cepheus OB3 molecular cloud to study its star formation (PID 20403, P.I. Pipher), and more recently by the GLIPMSE360 survey during the warm Spitzer phase as one of the large exploration programs (PID 60020, P.I. Whitney). The data, consisting of the basic calibrated frames or BCDs, have been recovered from the Spitzer legacy archive, version S18.18.0 (Cryo) and S19.0.0 (Warm). In all cases, the data was collected using the high-dynamic-range mode with a 12 s integration time for the 'long' frames (10.4 s on target) and 0.6 s for the 'short' ones. A summary of the observations is presented in table 1.
The BCDs were then reprocessed with the HiREs deconvolution software AWAIC 7 (A WISE Astronomical Image Co-Adder), developed by the Wide Field Infrared Survey Explorer (WISE) for the creation of their Atlas images (see e.g. Masci andFowler 2009, Jarrett et al 2012). The AWAIC software optimizes the coaddition of individual frames by making use of the point response function as an interpolation kernel, to avoid flux losses in undersampled arrays like those of IRAC, and also allows a resolution enhancement (HiRes) of the final image, by removing its effect from the data in the deconvolution process. We have used this method quite successfully in the HH 1/2 outflow (Noriega-Crespo and Raga 2012), and as mentioned above, a similar method has been used on Cep E (Velusamy et al 2011) and HH 46/47 (Velusamy et al 2007). On IRAC images, the HiRes enhances the angular resolution from the standard ∼ ″ 2 to ∼ ″ 0.6 -″ 0.8 (Velusamy et al 2008, Noriega-Crespo and. The combination of being deeply embedded and its youth, ∼5000 yr (Ladd and Hodapp 1997), perhaps makes Cep E one of the outflows where the morphology of the vibrational H 2 is nearly identical to the H 2 rotational emission observed with IRAC at 4.5 μm. This similarity has encouraged us to introduce an earlier H 2 v = 1 − 0 2.12 μm image from 1996 and one recently obtained in 2012, in the analysis of the proper motions, providing a ∼ 16 yr time baseline. The NIR 1996 image was obtained with the 3.5 m telescope at the Apache Point Observatory with a 256 × 256 array at f/5 and a ″ 0. 482 pixel −1 scale using a 2.12 μm filter (1% width) and 2.22 μm (4% width) to subtract the continuum. The complete analysis of these data was already presented by Ayala et al (2000). The 2012 image was obtained at Palomar Observatory with the wide field infrared camera (WIRC) mounted in the 200 in prime focus using a 2048 x 2048 Hawaii-II HgCdTe detector. WIRC has a field-of-view of 8. ′ 7 and a 0.2487 arcsec/pixel scale (Zhang et al 2013). The observations were carried out on 29 August 2012 in the 2.12 μm and K-continuum (2.27 μm, 2%) filters, with a 25 min total integration time. Figure 1 shows a comparison of the Palomar 2.12 μm continuum subtracted image with that at 4.5 μm from Noriega-Crespo et al (2004b). Among some of the small obvious differences are the lack of H 2 vibrational emission on the same region where there is a 'wide angle' cone at 4.5 μm (Velusamy et al 2011), which appears to be scattered light by small dust particles; and the emission at 2.12 μm on the south lobe that fits within these 'cones' and reaches further into the IRAS 23011 + 6126 central source. Other than these differences, the knots that we have

Proper motion measurements
As can be seen in figure 2, the Cep E outflow shows bipolar 'cavity' structures which extend out to ∼10-″ 15 from the source (with an approximately NNE-SSW orientation). These cavities have two ridges that first open out from the outflow source, and then converge into compact emission structures. Further away (∼ ″ 20 from the source), we find two bow-like prolongations of the outflow lobes. As can be seen in figure 2, the emission along the tworidged cavities shows a complex time-dependence, with the northern cavity becoming fainter and the southern cavity brightening from 2003 to 2010 (i.e., the period covered by the IRAC images).
We have defined four boxes, including the regions of convergence of the cavities and the bow-like structures, which we show in figure 4. For each of these four boxes, we have carried out cross correlations between the emission in the 2003.89 frame (the first of the IRAC frames, see table 1) and the other four available frames. From paraboloidal fits to the peak of cross correlation functions we then determine the offsets of the emission (within the four boxes, see figure 4) with respect to the 2003.89 frame. For the two 2.12 μm images, we performed the cross correlation initially with the 2003.89 frame as well, and then between themselves. In this way the offsets at 2.12 and 4.5 μm are measured on a common reference frame.
The resulting RA and DEC displacements for the four selected boxes are plotted as a function of time in figure 3, where the stars and plus signs correspond to the 2.12 and 4.5 μm frames respectively. From this figure, one notices that boxes 1 and 4 (corresponding to the bow-like structures, see figure 4) show substantial N-S motions, while the 'cavity tips' (boxes 2 and 3) have considerably lower proper motions. We have then carried out linear fits to the time dependencies of the RA and DEC offsets of our four boxes, from which we obtain the proper motions (and their associated errors) given in table 2 (2.12 μm) and 3 (4.5 μm).
The proper motions obtained from the 2.12 μm and IRAC 4.5 μm do show similar trends and are within a factor of two (table 2). One interesting thing to notice is that by placing all displacements on the same frame of reference, one can also measure a shift between the positions of the vibrational (2.12 μm) and rotational (4.5 μm) H 2 emission in the outflow. These are of the order of ″ 0.25 -″ 0.30 , certainly larger than the ″ 0.05 -″ 0.10 positional uncertainty of the selected regions. The H 2 proper motions show evidence, for first time, that both the Figure 3. The offsets in arcseconds in RA and DEC between the five epochs of Cep E as measured in four different boxes; two each for the north and south lobes of the outflow (see figure 4). The symbols correspond to the 2.12 μm (stars) and 4.5 μm (plus signs) data, respectively. vibrational and rotational H 2 can share the same kinematics within an outflow; and furthermore, that the vibrational and rotational H 2 gas does trace slightly different regions within the outflow. Although this last statement may seem obvious within the now classical scenario of the acceleration of a molecular outflow, where the atomic/ionic supersonic flow drags and/or excites the surrounding molecular gas (see e.g. Raga et al 1995), this is (to our knowledge) one of the first pieces of observational evidence that this is the case for the different components of the H 2 gas. Recall that spectroscopically we do have some good examples where both near-IR  Box  (Rosenthal et al 2000) where the short wavelength spectrometer detected 56 H 2 transitions within its 2.5-45 μm spectral range, and where the rotational emission was tracing a H 2 gas with an excitation temperature of 600 K, while the vibrational emission was tracing one at 3200 K. In OMC-1 it was possible to explain the bulk of the emission with collisional excitation produced by a combination of C-type and J-type shocks (Rosenthal et al 2000). In Cep E, the spectroscopic evidence (Moro-Martín et al 2001) also suggests a mixture of excitations to explain the optical and near/mid/far IR observations. The resulting tangential velocities are shown (together with their error boxes) in figure 4. Because of the multiple epochs, the motions seem to be better defined at 4.5 μm than at 2.12 μm.

Numerical simulations
In one of the first studies of the H 2 NIR emission in Cep E  the authors suggested that the 'wiggles and sideways positional offsets' were due to precession, with a relatively small precession angle of 4°. In the same study they noticed the presence of a couple of H 2 knots emanating westward from the central source and nearly perpendicular to the main Cep E outflow (Eislöffel et al 1996 figure 3), suggesting a very close-by second protostar, and therefore, a possible mechanism to drive the precession. To further test this hypothesis and compare with the kinematical behavior of the H 2 gas derived from the proper motion measurements, we present in this section some relatively simple 3D hydrodynamical simulations. The YGUAZU-A code (Raga et al 2000, Raga et al 2003 was selected for this simulation. The code, in a nutshell, uses a binary adaptive grid and integrates the gas-dynamic equations with a second-order accurate scheme (in time and space) using a flux-vector splitting method (van Leer 1982). The code has been used over a decade to simulate the gas dynamical processes that take place in several astrophysical scenarios including YSO outflows (e.g. Raga  Table 3. Cep E proper motions and tangential velocities from IRAC 4.5 μm images a . Box  , and MIRAʼs turbulent wake (Raga et al 2008) among others. In this version of the code we include a H 2 gas component, although we cannot distinguish between its vibrational/rotational excitation. The goal of the simulation is not to model in detail Cep E, but to show that the observed proper motions can indeed be explained when taking into account precession. Because of the relatively small dynamical age of the Cep E outflow (∼3000 yr), we set the initial conditions of the model as that of a jet emerging from a compact dense cold core with the simplifying assumption that it is in thermal balance with the surrounding medium ( = n 1000 core cm −3 , = T 1 core K, = n 10 ISM cm −3 , = T 10 ISM 3 K). Based on our previous simulations , we assume that the jet has an initial radius of × 1.5 10 16 cm, a temperature of 1000 K and a velocity of 200 − km s 1 , with a time dependent velocity variation ('pulses') of a 50% over a 60 yr period. The jet is set to precess on a 10°angle with respect to its cylindrical symmetry axis over a 1200 yr period. Since we are interested on the bulk motion of the gas to compare with what is observed in the proper motions, we have chosen a medium resolution computational grid of × × 128 128 256 cells set to a scale of = X Y Z ( , , ) × (5, 5, 10) 10 17 cm, respectively. The time resolution is set to 200 yr per step, so once the jet plunges through the core, it takes 15 to 20 frames to reach a dynamical age close to that of Cep E.
The results of the simulation for the H 2 number density are shown for five time steps starting at the moment when the jet finally breaks free from the dense core ( figure 5; step 45). In all cases we present the XZ projection, i.e. perpendicular and along the flow. In the H 2 number density one already can see that the initial effect of precession has been to widen the path of the flow and the creation of two different density maxima upstream. At about 2000-3000 yr later (steps 56 and 59), the maxima have nearly merged and the flow is compact and asymmetric.
In an effort to better compare the numerical model with the observations we have integrated in the 3D grid the H 2 emission along the line of sight and projected it on a°30 angle. We have taken the difference in the projected emission for models 56 and 59 (i.e. a 600 yr interval), to mimic as much as possible the proper motion measurements. The result (figure 6) Figure 5. The H 2 number density (from 0 (black) to × 1.2 10 3 cm -3 (white)) for a precessing jet emerging from a cold core to simulate the Cep E outflow.
shows at least four 'knots' with tangential velocities ranging from 10 to 50 − km s 1 in slightly different directions, although the bulk of the motion is away from the center of the grid. A 50 − km s 1 tangential velocity is certainly consistent with the value of ∼62 − km s 1 of the north lobe, that plunges deeper into the cloud, and is a bit smaller than the 94 − km s 1 value of the south lobe, where one detects at optical wavelengths HH 377, i.e the observed proper motions seem to reflect the difference in physical environment between the two lobes.

Summary
We have derived proper motions based on six IR images of the Cep E outflow: two ground based images obtained in 1996 and 2012 in H 2 v = 1 − 0 at 2.12 μm, plus another four obtained at 4.5 μm with the Spitzerʼs IRAC camera over the 2003-2010 time period.
We have defined four cross-correlation boxes that included the more compact emission structures. Two boxes for the tips of the two-ridged cavities (boxes 2 and 3, respectively, in figures 3 and 4 and those for the northern and southern bowshocks (boxes 1 and 4, respectively). The proper motion of the tip of ridges is complex in both the rotational and vibrational emission and it may reflect time dependent variation in their illumination or excitation. On the other hand, the proper motions of the main bowshocks are well defined; they are moving away from the central source along the symmetry axis with tangential velocities of ± 62.6 29.5 and ± 94.0 26.5 − km s 1 , respectively. The southern bowshock is detected also at optical wavelengths and is known as HH 377, and its proper motion has been measured (Noriega-Crespo and Garnavich 2001) rendering a tangential velocity of ± (107 14) − km s 1 , directed approximately along the outflow axis. This motion is roughly consistent with the proper motion that we have obtained for the IR emission of this object. The proper motions based on the 2.12 μm emission are about a factor of two smaller than those in the mid-IR. With a time baseline of ∼ 16 yr and an angular resolution of ∼ ″ 1 , a priori one does not have any reason to believe that this difference in magnitude is not real. If this is the case, then the offset between the vibrational and rotational H 2 emission, plus the difference in velocity, suggests a different physical 'layer' in the outflow where the vibrational H 2 gas is excited. That not all the molecular tracers originate in the same place in young stellar outflows, including Cep E, has been nicely illustrated by a recent study of water using Herschel space telecope observations (Tafalla et al 2013). For Cep E, for instance, the H 2 O ( − 2 1 12 01 1670 GHz), CO (J = 2-1) and H 2 emission (from IRAC 3.6 μm channel) at the same angular resolution ( ″ 13 ), show a very different spatial distribution along the flow axis (Tafalla et al 2013, figure 4). In this case, the CO emission peaks closer to the source, while the H 2 O and H 2 share the same distribution farther away from the source. This means that gas at a temperature of tens of Kelvins (from CO) resides at a different place than gas at hundreds of Kelvins (from H 2 O and H 2 ). A similar process could be taking place in our case, where the H 2 rotational emission, as measured by IRAC at 4.5 μm (i.e. S(11), S(10) and S(9) lines) is tracing a higher kinetic temperature than the vibrational H 2 traced by the 2.12 μm emission (see e.g. Giannini et al 2011, Neufeld et al 2009) and it is moving at a higher velocity as well.
As mentioned in section 1, the observed emission at 4.5 μm from young stellar outflows is likely produced by three rotational transitions of H 2 . In the case of Cep E, in particular, one expects only a small contribution from the dust continuum emission at these wavelengths based on what is observed spectroscopically at 5 μm in its north lobe (Noriega-Crespo et al 2004b). Thus the infrared proper motions imply that the molecular hydrogen gas in Cep E is moving supersonically, just like its atomic/ionic counterpart. That the H 2 gas can have large tangential velocities in young stellar outflows was first noticed in the HH 1-2 system , where velocities as high as 400 − km s 1 , comparable to those from the optical tracers, were measured. Large H 2 flow velocities derived from proper motion measurements are now certainly not uncommon in young stellar outflows (Chrysostomou et al 2000, Zhang et al 2013,