1. INTRODUCTION

The HH 1/2 outflow (discovered by ^{Herbig 1951} and ^{Haro 1952}) has played a fundamental role in the study of collimated flows from young stellar objects (YSOs), and the associated observational and theoretical work has been reviewed by ^{Raga et al. (2011)}. This system has two bright “heads”: HH 1 (to the NW) and HH 2 (to the SE), centered on the “VLA 1” radio continuum source (^{Pravdo et al. 1985}).

The VLA 1 source also has a jet/counterjet system visible at IR wavelengths (^{Noriega-Crespo & Raga 2012}) extending out towards HH 1 and 2. Optically, only the slightly blueshifted N jet (pointing to HH 1) is visible (^{Bohigas et al. 1985}; ^{Strom et al. 1985}), as shown in Figure 1. This optical feature has been called the “HH 1 jet”. Apart from the papers mentioned above, a limited number of papers have studied some of the characteristics of the HH 1 jet:

optical images and proper motions:

^{Reipurth et al. (1993)},^{Eislöffel et al. (1994)},^{Bally et al. (2002)},radio proper motions:

^{Rodríguez et al. (2000)},infrared images:

^{Davis et al. (2000)},^{Reipurth et al. (2000)},infrared spectra:

^{Eislöffel et al. (2000)},^{García López et al. (2008)}.

Some of the most striking characteristics of HH 1/2 are their proper motions (^{Herbig & Jones 1981}; ^{Eislöffel et al. 1994}; ^{Bally et al. 2002}; ^{Hartigan et al. 2011}) and time-variability (^{Herbig 1969}, ^{1973}; ^{Raga et al. 1990}; ^{Eislöffel et al. 1994}). The fact that there are now four epochs of HST images of HH 1/2, covering a time span of ≈ 20 yr (^{Raga et al. 2015a}, ^{b}, ^{c}; ^{2016a}, ^{b}, ^{c}) has allowed progress on both these issues.

^{Raga et al. (2016a}, ^{b)} used the 4 epochs of HST images to determine proper motions of HH 1 and 2, finding a small acceleration for the motion of HH 1 and a small braking for HH 2 (when comparing their proper motions to the ones of ^{Herbig & Jones 1981}). They also used the photometrically calibrated HST images (^{Raga et al. 2016c}) to evaluate the recent time-variability of the emission of HH 1 and 2 (comparing their line fluxes to the ones of ^{Brugel et al. 1981}).

For their study of HH 1/2 proper motions, ^{Raga et al. (2016a}, ^{b)} explored a new method for determining motions of angularly extended objects, based on a two-step process:

convolving the frames of the different epochs with wavelets of chosen widths,

spatially fitting the peaks in the (degraded angular resolution) convolved frames.

In the present paper, we apply this new method to the four available epochs of HH 1/2 HST [S II] images, in order to determine the proper motions and intensity variations of the knots along the HH 1 jet (which was not studied in the papers of ^{Raga et al. 2016a}, ^{b)}. We also present a detailed description of the method, and describe a pipeline (written in Python) developed for applying this method to observational or simulated emission map time-series.

The paper is organized as follows. § 2 reviews the methods that have been used to measure proper motions in CCD frames of HH outflows. § 3 presents the new method for deriving proper motions and intensities of extended structures, and describes the Python pipeline. § 4 describes the proper motions of the knots along the HH 1 jet, and § 5 the time-variability of the [S II] emission. § 6 describes the standard attempts at using the observed proper motions to reconstruct the history of the time-variability of the ejection and to predict the future evolution of the ejected material. Finally, the results are summarized in § 7.

2. PROPER MOTIONS AND TIME-VARIABILITES OF HH OBJECTS FROM CCD IMAGES

As far as we are aware, the first attempt at measuring positions and fluxes of condensations in CCD frames of HH objects was done by ^{Raga et al. (1990)}, who analyzed H*α* and [O III] 5007 images of HH 1/2. These authors found the then non-trivial result that even though they had only two stars in their CCD frames (and were therefore only able to compute a scaling, rotation and translation rather than a “real” astrometric calibration of the images) they obtained positions for the HH 1/2 condensations that coincided with the forward time-projection obtained with the photographic proper motions of ^{Herbig & Jones (1981)}.

^{Raga et al. (1990)} measured the positions (and peak intensities) of the HH 1/2 condensations by carrying out paraboloidal fits to the emission peaks seen in the images. This kind of “peak fitting” procedure (fitting mostly either a paraboloid or a Gaussian) has been extensively used for obtaining proper motions of HH outflows (see, e.g., ^{Eislöffel & Mundt 1992}, ^{1994}; ^{Eislöffel et al. 1994}).

^{Heathcote & Reipurth (1992)} tried a different method to obtain proper motions from CCD images of HH outflows. In their analysis of images of HH 34, they defined a box (including the emission of the HH 34 jet) within which they carried out cross-correlations between pairs of images. This method proved to be a major improvement in determining proper motions of HH outflows, since instead of relying on the positions of sometimes illdefined peaks, the proper motions were determined with the emission within a spatially more extended box. This process yielded a cross-correlation function with a much better signal-to-noise ratio (compared to the images themselves), the peak of which could be fitted to a many times surprising accuracy. This cross-correlation technique has become the standard method for determining proper motions of HH objects (see, e.g., ^{Curiel et al. 1997}; ^{Reipurth et al. 2002}; ^{Hartigan et al. 2005}; ^{Anglada et al. 2007}).

The main inconvenience of the cross-correlation method is the fact that one has to choose boxes of arbitrary shapes (mostly square boxes have been used), sizes and locations so as to include features that one judges to be well-defined “entities” within the images. This is of course inconvenient in images with complex structures of different sizes, and also somewhat problematic, since the determined proper motions clearly depend on the “cross correlation boxes” that have been chosen.

In a study of a planetary nebula, Szyskza et al. (2011) used the interesting method of covering the images with a regular array of cross-correlationboxes (which they call “tiles”). The shifts of the peaks of the correlation functions corresponding to these boxes then give a “proper motion map” of the whole field (actually, a low-intensity cut-off has to be imposed so as to avoid random motions in boxes with no visible emission structures). ^{Raga et al. (2012a}, ^{2013)} applied this method to HH objects (with the implementation of the method being presented in detail in the latter paper).

This method of cross-correlation “tiles” has the clear advantage that one only needs to define:

a size for the tiles,

a “beginning point” at which to begin to draw one of the tiles,

a “low intensity cutoff” necessary for the proper motions to be calculated.

There are of course many fewer free parameters than the ones involved in a “free choice” cross-correlation box scheme.

However, it is evident that there are complications in this method. Two of these are that:

the tiles sometimes include only part of an apparently coherent structure (algorithmical efforts to surmount this problem are described by

^{Raga et al. 2013}),identifiable features sometimes are shifted away from a tile into neighbouring tiles in the image pairs (so that a shift has to be applied to one of the images before applying the division into tiles, see

^{Raga et al. 2012a}).

3. MEASURING PROPER MOTIONS OF HH OBJECTS WITH A “WAVELET TECHNIQUE”: A PIPELINE

In order to try to avoid these problems, ^{Raga et al. (2016a}, ^{b)} proposed (and used) an alternative, two-step method:

convolving the images with a wavelet of a chosen size,

determining proper motions from spatial fits to the peaks in the convolved maps.

This method is, of course, a “peak fitting method”, but it also incorporates a spatial averaging (through the convolution with a wavelet function) such as is obtained with the “cross correlation method”. The only free parameter of this method basically is the half-width *σ* of the wavelet function (and of course, the choice of which peaks are identified as “pairs” in two different epochs!).

Convolving an image with a wavelet of half-width *σ* has three effects:

improving the signal-to-noise ratio at the expense of spatial resolution,

eliminating emitting structures with scales

*< σ*,eliminating structures with scales

*> σ*.

If one convolves images with functions similar to instrumental “point spread functions” (e.g., with a Gaussian), one eliminates small scale structures, but larger scale structures in the images still remain. It is, however, unlikely that proper motions determined on images convolved with Gaussians would be substantially different from proper motions measured on convolutions with wavelets. We prefer convolutions with wavelets basically because of the mathematical properties of wavelet decompositions, which allow partial rebuildings of images with arbitrary ranges of spatial scales (see, e.g., ^{Kajdic et al. 2012}). However, this feature is not used in the present proper motion determinations.

The choice of the particular form of the wavelet function does not affect the obtained results in a substantial way. In our implementation, we have chosen a “Mexican hat” wavelet:

where *σ* is the half-width of the central peak. This function has an approximately Gaussian central peak, surrounded by a negative ring (such that its spatial integral is zero). Together with the “French hat” wavelet, this is one of the standard “wavelet kernels”. For an astronomically oriented discussion of the properties of these wavelet kernels (along with graphic depictions) see, e.g., ^{Rauzy et al. (1993)}.

The convolved maps *I*
_{
σ
} are then calculated through the usual integral

where *I*(*x*
^{′}
*,y*
^{′}) is the original (i.e., not convolved) image, and (*x,y*) are the coordinates of the convolved image. The convolutions are carried out with a standard, “Fast Fourier Transform” method.

On the convolved image, we then carry out paraboloidal fits to intensity peaks, from which we determine the positions and intensities of the peaks. We then determine proper motions from the shifts of the positions between successive epochs.

We have developed a pipeline (written in Python) that:

reads an image;

convolves it with Mexican hat wavelets of thespecified

*σ*values;finds peaks (either chosen by the user, orsearches for all peaks above a given intensity threshold) and carries out paraboloidal fits;

has a “user confirmation and labelling” routine(with which the user can choose the relevant peaks);

identifies the same peaks in two or more imagesand calculates the proper motions (with linear least squares fits to the knot positions as a function of time).

Item number 5 allows for several possibilities:

the more straightforward one is to calculate proper motions for the knots identified by the user with the same label in the available epochs. This is of course appropriate for images with a small number of emitting knots;

to automatically associate the “nearest knots” detected in two successive epochs;

to search for the nearest knot (in the following epoch), but only in the general direction away from the outflow source.

It is also possible to use the wavelet spectrum of the individual knots in order to find the knot pairs that are morphologically closest to each other, and to then use the identified pairs to calculate proper motions. This kind of “morphological evaluation” using wavelet spectra has been studied in detail by Masciadri & Raga (2004, in the context of the search for exoplanets), but has not yet been implemented in our pipeline.

Finally, our Python pipeline has routines to produce appropriately labeled plots for publication. Figures 1-4 (see the following sections) were produced with these routines. After further testing and improvements to the user interface, the routine will be available to the community.

4. PROPER MOTIONS OF THE HH 1 JET

We have taken the four epochs of [S II] HST images of HH 1/2 described by ^{Raga et al. (2016a}, ^{b}, ^{c)} obtained in 1994.61, 1997.58, 2007.63 and 2014.63 (we have not analyzed the H*α* frames because the HH 1 jet is very faint in this line). Figure 1 shows a region of the 1997 frame including the position of the VLA 1 source, HH 1 and the HH 1 jet. The CohenSchwartz (CS) star, despite its strategic location, is apparently not associated with the outflow.

The analysis presented in this paper is restricted to the region around the HH 1 jet shown with a white box in Figure 1. The [S II] emission within this region in the four epochs is shown in Figure 2, with the knots labeled with identifications that correspond to the ones of ^{Reipurth et al. (2000)} and ^{Hartigan et al. (2011)}. Also, we have labeled knot B of the HH 501 jet (see ^{Bally et al. 2002}) with a lower case “b”. This outflow appears to have been ejected by another source in the vicinity of the HH 1/2 source (^{Bally et al. 2002}).

In Figure 3, we show the four [S II] frames after convolution with a *σ* = 4 pix wavelet (i.e., with a central peak with a full width of 0^{′′}
*.*8). In these convolved frames, the jet breaks up into knots with well defined peaks, to which we fit paraboloids (giving peak fluxes and the positions shown in Figure 3).

With the positions measured for the successive knots (some of them seen in all frames, but others in only two or three frames) we carry out linear least squares fits to determine their proper motion velocities. These velocities are given in Table 1 (for a distance of 400 pc to HH 1/2) and are shown in Figure 4.

Knot |
v
_{
x
}
^{a}[kms ^{−1}] |
v
_{
y
}
^{a}[kms ^{−1}] |

A | 128 (91) | 86 (45) |

B | 245 (12) | −6 (4) |

C | 258 (17) | 3 (4) |

D | 235 (12) | 5 (6) |

E | 271 (12) | 9 (2) |

F | 259 (12) | 18 (2) |

G3 | 286 (13) | 9 (5) |

G1 | 288 (16) | 17 (2) |

H | 247 (38) | 38 (19) |

b | 149 (12) | −24 (5) |

^{a}The values in parenthesis are the estimated errors.

The errors for the proper motion velocities given in Table 1 are calculated as follows. We estimate that the errors of the fits to the knot positions are at most 1 pixel (0^{′′}
*.*1). This estimated error is used to calculate the errors in the proper motion velocities of knots A and H, which have measured positions in only two frames. For all of the other knots, we calculate the error in the knot positions using the standard deviation of the measured positions with respect to the (straight line) least squares fits.

These errors (as is standard) are then propagated using the covariance matrix to calculate the errors in the proper motion velocities. Therefore, the quoted errors (see Table 1) correspond to estimates of the standard deviations.

It is clear that we do not obtain a significant result for knot A (as the errors are comparable to the determined proper motions). This is not surprising because we only see this knot in the two first epochs (which only have a time range of 3 years). For the remaining knots we do obtain proper motions along (*v*
_{
x
} ) and across (*v*
_{
y
} ) the outflow axis with reasonable errors (ranging from ≈ 10 to 40 km s^{−1}, see Table 1).

5. THE INTENSITIES OF THE KNOTS IN THE HH 1 JET

Figure 5 shows the peak [S II] intensities for knots A-I of the HH 1 jet in all epochs, obtained through paraboloidal fits to the peaks of the convolvedimages (see § 4). It is clear that for distances from the VLA 1 source larger than *x* ≈ 7^{′′} there is a general trend of decreasing intensities as a function of *x*. This trend approximately follows an *I* ∝ *x*
^{−3} power law (shown with a dashed line in Figure 5).

Our observations do not show in a conclusive way that individual knots have intensities that “slide down” the *x*
^{−3} slope as a function of time. This is because in the 2014 frame (intensities shown with open circles in Figure 5) we obtain systematically larger intensites for all knots than in the 2007 frame. This is a result of the fact that the HH 1 jet region has a relatively strong reflection nebula, with peak intensities aligned with the jet. This reflection nebula has a stronger contribution in the 2014 frame, which was obtained with the WFC3 camera (with a [S II] filter of 118 ˚A width). The first three epochs were obtained with the WFPC2 camera (with a [S II] filter of 47 ˚A width), and have a smaller contribution from the reflection continuum.

The observed *I*
_{
[SII]
} ∝ *x*
^{−3} dependence for large distances along the HH 1 jet (see Figure 5) is in remarkable agreement with the prediction of the analytic, “asymptotic regime” of periodic internal working surfaces of ^{Raga & Kofman (1992)}. These authors note that at large enough distances from the source, the decaying working surfaces should have an intensity *I* ∝ *x*
^{
−(κ+1)
} , where *κ* is the index of an assumed power law dependence of the line emission as a function of shock velocity (see equation 25 of ^{Raga & Kofman 1992}). If one takes the planeparallel shock models of ^{Hartigan et al. (1987}), from the lower range of the shock velocities of their models one obtains that the [S II] intensity has a scaling ^{Raga & Kofman (1992)} predicts a [S II] intensity ∝ *x*
^{−3}, in surprisingly good agreement with our observations of the HH 1 jet (see Figure 5).

6. THE PAST AND FUTURE EVOLUTION OF THE HH 1 JET

As can be seen in Table 1, along the HH 1 jet we see a general trend of decreasing velocities with increasing distances from the outflow source. Such a decreasing velocity trend could in principle be the result of drag due to entrainment of stationary, environmental material.

It is clear that in some of the “parsec scale HH jets” (e.g. in HH 34, see ^{Devine et al. 1997}) a progressive decrease in proper motions for the “heads” at larger distances are seen, and that this trend cannot be explained as a result of a secularly increasing ejection velocity from the outflow source (^{Cabrit & Raga 2000}). This rather dramatic slowing down of the HH 34 “heads” is due to the fact that a precession of the outflow axis results in a direct interaction of the successive heads with undisturbed environmental material (^{Masciadri et al. 2002}).

As the knots along the HH 1 jet are very well aligned, we would not expect them to slow down due to frontal interaction with the surrounding, stationary environment (as occurs in the giant HH 34 jet, see above). One might still have “side entrainment” into the HH 1 jet, resulting in some amount of slowing down at increasing distances from the source. This effect has recently been evaluated by ^{Raga (2016)}, in terms of a somewhat uncertain “*α* prescription” for the entrainment velocity), who finds that in order to obtain a substantial slowing down one needs a surrounding environment (in contact with the jet beam) ≈ 10 to 100 times denser than the jet. This is unlikely to be the case in the optically visible HH 1 jet, which has already emerged from the dense core surrounding the outflow source. Also (as discussed by ^{Raga 2016}), the effect of buoyancy (which includes the gravity and the environmental pressure gradient) is negligible for the high velocities of HH outflows.

We therefore interpret the decreasing proper motion velocities (with increasing distances from the outflow source) along the HH 1 jet as ballistic motions resulting from an increasing ejection velocity as a function of time. We then take the positions of the HH 1 jet knots in the 2007 frame, and calculate the dynamical ejection times

where *x* is the distance from the outflow source and *v* is the proper motion velocity. In Figure 6, we then plot *v* as a function of *t*
_{
dyn
} , which is the “ballistic knot” prediction of the past ejection time-variability history of the outflow source. The ejection velocity has a general trend of increasing velocities towards more recent times, which we fit with a straight line, giving:

where *u*
_{0}(*τ*) is the ejection velocity in km s^{−1} and *τ* is the ejection time in years (*τ* = 0 corresponding to 2007, since we have used the knot positions of this epoch). The linear least squares fit (equation 4) has been calculated with the method described in Appendix A.

The ejection velocity clearly should also have a short-term variability that produces the knots that we observe along the HH 1 jet, so that the trend of equation (4) would actually correspond to a longterm variability superimposed on the “knot producing” mode (see, e.g., ^{Raga et al. 2015c}).

An interesting question is whether the relatively low velocity of knot H (see Table 1 and Figure 6) is evidence that the more recently ejected, optically detected material is starting to show a decreasing ejection velocity *vs.* time trend. Given the large error of our H knot proper motion (see Figure 6), it is hard to conclude that this is indeed the case.

Also, we can obtain a second estimation of the motion of knot H as follows. We take the separation between knots H and F in the 1998.15, [Fe II] 1.64 *µ*m image of ^{Reipurth et al. (2000)} (in which knot H is already visible), and compare it with the separation between these two knots in our 2014.63 image. From this comparison, we find that knot H has an axial motion 17±2 km s^{−1}, faster than knot F. Combining this result with our knot F proper motion (see Table 1), we obtain a (276 ± 18) km s^{−1} velocity for knot H, which is consistent (within the errors) with the proper motion obtained from the optical images (see Table 1), but does not support the existence of a drop in ejection velocity associated with this knot.

In order to use the present day positions and proper motions to predict the future evolution of the HH 1 jet, in Figure 7 we plot the ballistic trajectories of the knots on a (*x,t*) plot (where *x* is the position of the knots as a function of time *t*). The *t* = 0 axis corresponds to the 2007 knot positions. For knot H (the trajectory with the smallest *x* at *t* = 0), we have used the 276 km s^{−1} velocity estimated from IR images (see above).

In Figure 7, we see that the knot trajectories have crossing points in the *x* = 0 → 80^{′′} distance range, at times smaller than ≈ 500 yr. Therefore, by the time the HH 1 jet knots have reached the presentday position of HH 1 (also shown in Figure 7), many knot-merging events will have occurred. This result is similar to the one found by ^{Raga et al. (2012a)} for the HH 34 jet.

In order to visualize the effect of the knotmerging events, we use the simple momentum conserving knot-merging model of ^{Raga et al. (2012b)}. We take the 2007 HH 1 jet knot positions and velocities, and assign equal masses to all knots. We then follow the knot trajectories, merging colliding knots using mass and momentum conservation conditions. The knots are not assigned sizes, so that knot collisions take place at the points of trajectory crossings (see Figure 7).

The result of this exercise is shown in Figure 8. In this figure we show the knot positions at 150 yr intervals (*t* = 0 corresponding to 2007). The knots are represented as circles (centred on the knot positions) with radii proportional to the mass of the knots. It is clear that by *t* ≈ 450 yr most of the knots have merged, and that at this time the position of the merged knots is slightly upstream of the present position of HH 1 (at *x* ≈ 4*.*5 × 10^{17} cm or 75^{′′} at a distance of 400 pc). Therefore, the material that is being ejected now (in the HH 1 jet) will eventually form a new “head” close to the present-day position of HH 1.

7. SUMMARY

We present a discussion of the methods used in the past for measuring proper motions of HH outflows (§ 2), and a description of a new method recently developed (§ 3). This method has been implemented in a Python pipeline.

We then use this pipeline to determine proper motions of the knots along the HH 1 jet in the four available epochs of [S II] HST images (obtained in 1994, 1997, 2007 and 2014). We find proper motions that are well aligned with the outflow axis, and with values ranging from ≈ 230 to 290 km s^{−1}, with the faster velocities mostly in the knots closer to the outflow source (see § 4 and Table 1).

For each knot, we calculate a dynamical ejection time, and then plot the outflow (proper motion) velocity as a function of ejection time (see Figure 6). This plot shows that (under the assumption of ballistic knot motions) the outflow velocity has increased towards more recent times.

This “acceleration” of the ejection (as a function of ejection time) implies that the knots presently observed along the HH 1 jet will merge into a large working surface. This is seen in the crossings of the ballistic knot trajectories (see Figure 7) and in the momentum/mass conserving, “merging knot” model shown in Figure 8. From this model, we see that in ≈ 450 yr most of the HH 1 jet knots will have merged, and that at this time the position of the merged knots will be slightly upstream of the present position of HH 1 (at *x* ≈ 4*.*5 × 10^{17} cm or 75^{′′} at a distance of 400 pc). This result is qualitatively consistent with the suggestion of ^{Gyulbudaghian (1984)} that the diverging proper motions of the condensations of HH 1 directly imply that it was formed not far upstream from its present-day position.

This kind of morphology (a large working surface at large distances, and a short chain of knots that will merge at the position of the present-day large working surface) is to be expected from models of two-mode ejection variabilities (see the analytic discussion of ^{Raga et al. 2015c}). Also, a qualitatively most similar situation has been previously found for the HH 34 outflow (see ^{Raga et al. 2012a}, ^{b}), in which the knots along the jet will merge when they reach the present-day position of HH 34S.

An important question is whether or not this kind of configuration (of knots along a jet predicted to merge at the present-day position of a large “head”) is found in other HH jets. For some HH outflows, it is possible that proper motion data of sufficient accuracy might be already available, and a detailed study of the available data might yield interesting results. In other cases, future observations might be necessary in order to resolve this question.