Ejection History of the IRAS 04166+2706 Molecular Jet

TL;DR

This paper investigates the velocity structure of the IRAS 04166+2706 molecular jet, revealing a saw-tooth pattern explained by a spherically expanding wind model and supported by numerical simulations, shedding light on its ejection history.

Contribution

It introduces a geometric model of a spherically expanding wind to explain the jet's velocity pattern, supported by numerical simulations, providing new insights into the jet's ejection history.

Findings

The velocity pattern can be explained by a spherically expanding wind model.

The inclination angle of 52° explains the observed velocity features.

Numerical simulations support the geometric interpretation.

Abstract

The high-velocity molecular jet driven by Class 0 protostar IRAS 04166+2706 exhibits a unique saw-tooth velocity pattern. It consists of a series of well-aligned symmetric knots with similar averaged speeds, whose speeds at peaks of emission decreases roughly linearly away from the origin. Recent ALMA observations of knots R6 and B6 reveal kinematic behavior with expansion velocity increasing linearly from the axis to the edge. This pattern can be formed by a spherically expanding wind with axial density concentration. In this picture, the diverging velocity profile naturally possesses an increasing expansion velocity away from the axis, resulting in a tooth-like feature on the position-velocity diagram through projection. Such geometric picture predicts a correspondence between the slopes of the teeth and the outflow inclination angles, and the same inclination angle of 52$^\circ$ of the IRAS 04166+2706 can generally explain the whole pattern. Aided by numerical simulations in the framework of unified wind model by Shang et al. (2006), the observed velocity pattern can indeed be generated. A proper geometrical distribution of the jet and wind material is essential to the reconstruction the ejection history of the system.

Equations

$$d_L=D_0\cos (\theta +\alpha ),$$

$$d_R=D_0\cos (\theta -\alpha ),$$

$$\eta =\frac{v_L-v_R}{d_L-d_R}=\left(\frac{V_0}{D_0}\right)\frac{\sin (\theta +\alpha )-\sin (\theta -\alpha )}{\cos (\theta +\alpha )-\cos (\theta -\alpha )}.$$

$$\eta =\left(\frac{V_0}{D_0}\right)\frac{-1}{\tan \theta }.$$

$$\bar{v}=$$

$$\bar{d}=$$

$$\frac{\bar{v}}{\bar{d}}=\left(\frac{V_0}{D_0}\right)\tan \theta .$$

$$\eta =\left(\frac{\bar{v}}{\bar{d}}\right)\frac{-1}{{\tan }^2\theta }.$$

$$\hat{\eta }\equiv \frac{\eta }{\bar{v}/\bar{d}}=\frac{\Delta v/\bar{v}}{\Delta d/\bar{d}}=\frac{(v_L-v_R)/(v_L+v_R)}{(d_L-d_R)/(d_L+d_R)}.$$

$$\hat{\eta }=\frac{-1}{{\tan }^2\theta }.$$

Full text

ISM: individual objects (IRAS 04166$+$2706) $–$ ISM: jets and outflows $–$ ISM: kinematics and dynamics $-$ stars: formation

1 Introduction

Mass outflow phenomena have been an integral part of the star formation picture since the first discovery of Herbig-Haro objects in 1950s (Herbig, 1951; Haro, 1952). During the gravitational collapse of dense molecular cloud, a rotating disk forms around the protostar along with the infalling gas. Theoretical considerations suggest that both mass accretion onto the protostar and the driving of a high-velocity jet from the system may occur through the inner disk under the influence of magnetic fields (e.g., Shu et al., 2000; Bally, 2016). Studies of outflows may therefore help shed light on the accretion activities of a young protostellar system.

Observations have revealed increasingly more details of outflow morphology and kinematics. They have been generally classified as parabolic shell-like or collimated jet-like according to their morphology, which also motivated two types of theoretical models. The group of outflows associated with the youngest Class 0 protostars (Andre et al., 1993) appear to be more collimated in general and possessing additional emission peaks at higher velocities than the classical outflows featuring broad wings in the spectra. Such component is referred to as the extremely high velocity (EHV) component to be distinguished from the classical, standard high-velocity gas (SHV) (Bachiller, 1996). High resolution interferometric mapping further revealed that the EHV gas forms highly collimated jet-like structures along the axis while the SHV gas forms wider, cavity-like structures surrounding the former. Examples of such Class 0 outflows showing both EHV and SHV components include L1448C (e.g., Hirano et al., 2010), HH 211 (e.g., Gueth & Guilloteau, 1999), and HH 212 (e.g., Lee et al., 2007). The simultaneous presence of a high-velocity jet and low-velocity shell in the young outflows motivate models that can explain the formation of both features like the unified wind model of Shang et al. (2006).

A knotty appearance is common to the EHV jet in the youngest Class 0 outflows, which are resolved into a series of well-aligned emission peaks under high resolution observations through the SMA (Ho et al., 2004) or ALMA. One of the plausible knot-forming mechanism is through a variable velocity history of the flow. Signatures of shocks have been detected in some sources, which suggests an initial variability in the ejections. Some proposed physical origins for the quasi-periodic variability include: stellar magnetic cycles or global magnetospheric relaxations of the star-disk system (of few to tens of years time scale), potential perturbations by binary companions, or EXOr-FUOr outbursts (see Frank et al., 2014, and the references therein).

The molecular outflow driven by the class 0 protostar IRAS 04166$+$2706 (hereafter I04166) was the first highly collimated, extremely high velocity bipolar outflow identified in the nearby Taurus molecular cloud (140 pc) (Tafalla et al., 2004). Its CO $J=$ 2 – 1 emission show an overall length-to-width ratio of $>10$, and its EHV peaks are clearly detected in the spectra at velocities of $>30$ km s*-1*, higher than the classical outflow wings of $<10$ km s*-1* relative to the system velocity. The two spectral features are loosely connected by rather weak intermediate emissions.

The rather distinct morphology of the two components was revealed by the observations of Santiago-García et al. (2009). At low velocities, a symmetric pair of red and blue-shifted V-shaped cones are identified with their apex overlapping near the position of the central driving source. On the other hand, a series of well aligned emission peaks are found in the high-velocity channels, forming a collimated jet-like structure. Even considering the gradually increasing size of the high-velocity peaks, its $\simeq 10$° widening is still much smaller than the $\simeq 32$° opening angle of the low-velocity cones. Because of its relatively clean appearance, I04166 may be regarded as an archetypical example of the youngest molecular outflows.

More interestingly, a peculiar kinematic pattern is revealed in the position–velocity (PV) diagram cut along the axis of the I04166 outflow (Santiago-García et al., 2009). We plot the same PV diagram in Figure 1, which clearly shows that the EHV emission consists of a series of peaks like a sawtooth. While having a similar average speed of $\sim 40$ km s*-1*, each peak possesses an internal gradient where the observed velocity drops about linearly with distance from the driving source. The “slow head”–“fast tail” behavior found in each knot entity is not easily understood since one would naively expect to see higher velocity gas near the head simply because it is faster.

The discrete emission peaks of the I04166 jet suggests variable or episodic ejection activities in the past. Inferring the velocity history of the system is, however, not trivial. This is because an interpretation of the observed line-of-sight velocity will depend on a correct understanding of the flow geometry. Clues of the I04166 jet geometry may come from the two extended knots B6 and R6, whose kinematics have been mapped in details using the Atacama Large Millimeter/submillimeter Array (ALMA) by Tafalla et al. (2017). Their modeling showed that the observed velocity field is consistent with a slightly parabolic disk with expansion velocity increasing linearly from the center to the edge. They are consistent with the same interpretation of Wilson (1984), Raga et al. (1990), and Stone & Norman (1993), in which the expansion of a knot is driven by high pressure, shocked gas in an internal working surface. Numerical simulations of pulsed jets can produce a sawtooth-like velocity pattern similar to that observed in the I04166 outflow (see Figure 16 of Stone & Norman, 1993), with a higher jet velocity (235 km s*-1*) and high gas temperature ($10^{4}$ K).

In this work, we propose an alternative framework of a spherical wind with mass concentration around the outflow axis. Not only can it naturally explain the tooth-like velocity pattern of individual knots, the systemic variation of tooth slopes with distance across the entire jet is also automatically reproduced. Although the notion that a molecular jet presents a spherical wind-like velocity field may seem contradicting, this is actually the expected property in magnetocentrifugal wind theories like the X-wind model (Shu et al., 2000). Besides formulating an analytic description of the sawtooth velocity pattern, we also carry out numerical simulations in framework of the unified wind model of Shang et al. (2006). The model consists of a purely radially directed wind with cylindrically stratified density profile that is highly concentrated toward the axis. This is the asymptotic behavior predicted by the X-wind model and turns out to be ideal for explaining the I04166 jet.

Several observations toward the I04166 outflow are critical leading to current work. Santiago-García et al. (2009) observed the I04166 outflow with the IRAM Plateau de Bure interferometer (PdBI) in its CD configuration during 2004–05 (blue lobe) and 2005–06 (red lobe). CO $J=$ 2 – 1 (230.5 GHz) and SiO $J=$ 2 – 1 (86.8 GHz) emissions were mapped at angular resolutions of between 2″  and 4″. Single-dish data were added to the visibility to improve the quality of the CO $J=$ 2 – 1 images, and the beam-size of the combined data is $3\arcsec$. Tafalla et al. (2017) observed the two knots B6 and R6 of the I04166 outflow at 230.5 GHz using ALMA in December 2014. The array was in its C$32-2$ configuration and no observations with the compact Morita Array were made because of the low declination of the source. CO $J=$ 2 – 1, SiO $J=$ 5 – 4 and SO $J_{N}=$ $6_{5}$–$5_{4}$ transitions were observed simultaneously, and the beam-size of the CO data is $1\farcs 5\times 1\farcs 1$. The CO $J=$ 3 – 2 data of Wang et al. (2014) taken in 2010 by the SMA covered the inner four pairs of knots and the synthesized beam-size is $\sim 1\farcs 0\times 0\farcs 8$.

We organize the paper as follows. In Section 2, we explain how the proposed model explains the observed velocity pattern. An analytic description of the pattern slope is developed, which relates the slope of the tooth pattern to the underlying outflow inclination angle. In Section 3, we present numerical simulations based on the framework of unified wind model and realistically reproduce the observed pattern. Implications on the outflow velocity history are discussed in Section 4. Section 5 summarizes the results.

2 Origin of the Velocity Pattern: An Analytic View

The nature of the sawtooth velocity pattern on the Position-Velocity diagram (PV Diagram) observed in the I04166 jet has posted an intriguing challenge to models. Each emission peaks of both the blue and red shifted lobes presents an internal velocity gradient of “slow head–fast tail” appearance on the PV diagram along the outflow axis. Moreover, the slopes of these tooth features (the apparent velocity gradient on the PV diagram) appear to vary systematically with distance from the driving source, where the inner ones are steeper than the outer ones (Figure 1). Such ordered behavior is intriguing and is in need of an explanation. In this section, we propose an analytic model in which each emission peak of the I04166 jet corresponds to a pulse of ejection of a spherically expanding wind with an axial density concentration. We show here how the velocity pattern can be naturally explained in this picture. A quantitative description for the slope will be developed, which provides some insight into the model.

A cartoon illustration of the model and its corresponding PV diagram is shown in Figure 2. The star symbol in the upper panel indicates the driving source, and the outflow axis is inclined from the plane of the sky by an angle $\theta$. The ejecta consist of a spherically expanding layer of gas whose velocity is purely in the $r$-direction, pointing away from the origin as indicated by the gray arrows. The mass is only contained within a small opening angle $\alpha$ from the outflow axis, so in three dimensions it will look more like a slightly curved circular plate than a spherical shell. The illustration only shows its cross section in two dimension. In radio observations, an observer (viewing from below) will see a red-shifted lobe in this case.

The lower panel of the same figure illustrates the expected velocity pattern of the model, where a tooth-like pattern is formed. The origin of this quasi-linear velocity gradient can be intuitively understood by considering the gas residing at the left and right edges of the ejecta which we label as $L$ and $R$. In this geometric model, it is apparent that the gas on the left edge will be found at a smaller projected distance from the driving source than the other on the right edge ($d_{L}<d_{R}$). The projected velocity, on the other hand, will be larger for gas at the left edge than that at the right ($v_{L}>v_{R}$). This is because the ejecta expand at the designed speed but the flow direction of the left edge is more inline with the line-of-sight direction. Together, this results in an apparent velocity pattern of “slow head–fast tail” behavior on the PV diagram. In this picture, the observed velocity gradient is a natural consequence of projection.

For a more quantitative understanding, we characterize a tooth pattern by its slope on the PV diagram by considering the two points $L$ and $R$. For an inclination angle $\theta$ from the plane of the sky, the observable quantities can be written as

[TABLE]

where $D_{0}$ is the intrinsic distance of the ejecta from the source, $V_{0}$ is its intrinsic speed, and $\alpha$ is the half opening angle from the axis. The slope $\eta$ of the corresponding velocity pattern is defined as difference in line-of-sight velocity $\Delta v$ over the difference in projected distance $\Delta d$ in the plane of the sky, namely, $\eta\equiv\Delta v/\Delta d$. Considering the two points $L$ and $R$, we have

[TABLE]

With some algebra, we arrive at a simple relation:

[TABLE]

Note in the current formulation, a positive $\theta$ value describes a red-shifted gas while a negative $\theta$ corresponds to a blue-shifted one. The apparent position offsets $d_{L}$ and $d_{R}$, on the other hand, are considered positive for both cases, which is different from the continuous coordinate in Figure 1 where one of the two lobes will be negative. Since $D_{0}$ and $V_{0}$ are positive, according to Equation 1 the slope $\eta$ will be negative for a red-shifted gas and positive otherwise. Note the opening angle $\alpha$ happens to cancel out in this approximation.

Aside from the outflow inclination angle, Equation 1 clearly shows that the slope of a tooth feature is determined by the intrinsic velocity of the ejecta and distance from the driving source. This relation is most interesting when applied to a series of knots in the same outflow. In case of a fixed inclination angle and similar intrinsic speed, the formula predicts that $\eta$ should varies inversely with the distance of knot from the origin. This means that the tooth pattern farther away from the center should present a smaller $|\eta|$ and appear shallower or less steep. This agrees with the trend seen in the I04166 jet.

Figure 3 intuitively explains why this trend is naturally produced in the current geometric model by comparing two pieces of ejecta $a$ and $b$ at different distances from the origin. For simplicity, we assume that the two ejecta have same intrinsic velocities and sustain same opening angles, which is a good approximation for the I04166 system. In this case, the difference between the line-of-sight velocity components at the left and right edges will be exactly the same for the two ejecta: $|v_{L}-v_{R}|$. However, as is shown in the illustration, the outer knot $b$ extends a larger projected distance than the inner one $a$, namely, $|d_{Lb}-d_{Rb}|>|d_{La}-d_{Ra}|$. Hence the velocity pattern of $a$ will turn out to be steeper than that of $b$ because the same velocity difference is distributed across a smaller distance difference in $a$ compared to $b$. Note that the apparent steepness of the pattern depends on the velocity change per unit length ($\Delta v/\Delta d$) and not the actual size of the ejecta. For example, increasing the opening angle of $a$ would result in a larger knot but the pattern on the PV diagram will simply extend and does not become shallower. It would still appear steeper than the outer tooth $b$ regardless of the projected size.

To more directly describe variation of pattern slope on a PV diagram, we can replace $V_{0}$ and $D_{0}$ in Equation 1 by observable quantities as follows. Let $\bar{v}$ and $\bar{d}$ denote the mean velocity and distance of the two representative points, namely, $\bar{v}\equiv(v_{L}+v_{R})/2$ and $\bar{d}\equiv(d_{L}+d_{R})/2$. We have

[TABLE]

And therefore,

[TABLE]

Substituting this into Equation 1 and eliminate $V_{0}/D_{0}$, we derive another expression for the slope:

[TABLE]

Equation 2 expresses $\eta$ in terms of the mean projected velocity and mean projected distance of the two points $L$ and $R$. This explicitly shows that the slopes of the ejecta should vary inversely with their position offset from the central source for a same inclination angle and similar mean observed velocity $\bar{v}$.

This prediction can be tested with the I04166 jet. Figure 4 shows a comparison between model and the observational data. At each knot position, we overlay three sets of velocity pattern computed with different inclination angles on top of the observed CO $J=$ 2 – 1 data. Given an inclination angle, we can compute the intrinsic distance $D_{0}$ and the intrinsic speed $V_{0}$ of each ejected blob from its observed position and velocity. And with the values of $D_{0},V_{0}$ and $\theta$ the velocity pattern is decided. In this demonstration we display the model pattern up to a half-opening angle of $\alpha=7.5$° from the rough center of each tooth.

The thick dark lines show the resulting pattern of the model that assumes an inclination angle of $38$° from the plane of the sky ($52$° from the line-of-sight). This is the average angle estimated by Tafalla et al. (2017) based on the elliptical shapes of the knots B6 and R6 in their ALMA observations. Figure 4 shows that such angle is not only suitable for predicting the pattern of these two particular knots, but also results in good fits for most others. On the other hand, the patterns computed with $\pm 10$° from this value are less compatible with the observation. The predicted slope is too steep in the case of $28$° and too shallow in the case of $48$° . The overall behavior of the whole sawtooth velocity pattern can therefore be generally explained with one same inclination angle of $\simeq 38$° from the plane of the sky. Since the model dictates that the pattern slope should vary inversely with knot distance from the center (Equation 2), the simultaneous agreement across the knots is actually not trivial. The overall agreement implies that the proposed geometric model does capture the feature of the I04166 jet kinematics.

Figure 4 also shows that the model velocity pattern are not exactly linear. They are slightly curved and this can be understood with Equation 1. Consider a single dark line segment computed for an ejected blob of fixed $D_{0}$ and $V_{0}$ with an inclination angle of $\theta=38$° from the plane of the sky. Although $D_{0}$ and $V_{0}$ are the same for different points along this stripe, the angle of $38$° does not apply to the entire piece. As we have assumed a half opening angle of $7.5$° , the upper and lower ends of the model tooth pattern on the PV diagrams are actually inclined by $45.5$° and $30.5$° , respectively. According to Equation 1, $|\eta|$ should vary inversely with $\tan\theta$ for a fixed $V_{0}/D_{0}$. Therefore at the lower end of the velocity pattern, $|\eta|$ is always $\sim 1.73$ times larger than that at the upper end in this example. Because the precise inclination angles of different parts of an ejected blob are different, the pattern appears curved.

Finally, we go one step further and absorb the term $\bar{v}/\bar{d}$ into $\eta$. A dimensionless slope $\hat{\eta}$ can then be defined as

[TABLE]

From Equation 2, we have

[TABLE]

The dimensionless slope is simply a function of the outflow inclination angle $\theta$ from the plane of the sky. Its value is negative for both the red and blue-shifted gas.

Since $\hat{\eta}$ is an observable quantity that can be calculated from the observed tooth-like velocity pattern, Equation 3 actually provides a way to estimate the outflow inclination angle. For any two chosen points on a tooth pattern, it gives you the average inclination angle of the two specific points (see Figure 2). To estimate the inclination angle of the outflow axis, we therefore need to have $L$ and $R$ symmetrically displaced from the outflow axis by a same half-opening angle. Another caveat is that this method only applies when the pattern forms by a spherical wind-like velocity field as assumed. It will not just apply to any linear velocity pattern on a PV diagram.

To see how this method works in practice, we apply it to the pattern of I04166. We manually pick out two points for each tooth feature as shown in the upper two panels of Figure 5 where the selected points are labeled with open triangles. The corresponding dimensionless slopes are computed for each tooth and plotted on the bottom left panel. The value of $\hat{\eta}$ falls in the range of $\sim-3$ to $-1$. The bottom right panel plots the corresponding inclination angles computed using Equation 3. The mean inclination angle values are $38.4$° and $37.6$° for the blue and red shifted lobes, respectively, which seems to be a reasonable estimate. Note that there is a trend where $\hat{\eta}$ (and $\theta$) appear to decrease with distance from the source which is not really expected.

Given the finite resolution of observations, we might not always be able to choose perfect representative points for a tooth feature that best reveal its inclination. To assess the uncertainty associated with this potential error, we displace the upper selected point in each pattern to the left and right by 140 AU ($\sim$1″  at 140 pc distance) and compute the associated $\hat{\eta}$ and $\theta$ again. The varied values are then presented as the upper and lower bounds of the error bars in the panels of $\hat{\eta}$ and $\theta$ in Figure 5 to illustrate the associated error range. The larger error bars for the inner data points indicate that they are more sensitive to the position measurement. This is because the inner knots generally present a steeper velocity pattern and a smaller $\Delta d$. Therefore a same positional error would cause a more significant change in the estimated $\theta$. Considering the size of the error bar, the inferred $\theta$ from most of the knots appear to be consistent with the mean value of $\sim 38$° . The decreasing inferred $\theta$ from inner to outer knots, however, is puzzling as it appears too systematic to be random measurement errors. Since it is unlikely that the inclination angle of the I04166 outflow has really changed by $15$° over the years, this probably suggests that there are some missing elements in the current simple model. One potentially possible cause of this behavior is mentioned in the next section.

3 Numerical Simulations of Variable Velocity Wind

The kinematics of the I04166 molecular jet is compatible with a spherical wind of axial density concentration. We will further demonstrate the scenario by reproducing such velocity patterns using numerical simulations. Our simulations build upon the unified wind model framework outlined in Shang et al. (2006) with detailed simulation setups following the two-temperature generalization described in Wang et al. (2015).

3.1 The Unified Wind Model for Young Protostellar Outflows

In the framework of the unified wind model, a toroidally magnetized wind is launched into an ambient envelop of toroid. The primary wind is constructed based on the asymptotic behavior of the $X$-wind (Shu et al., 2000, 1995), which is characterized by a radially directed velocity field and cylindrically stratified density profile that is highly peaked toward the axis ($\rho\propto\varpi^{-2}$, where $\varpi$ is the cylindrical distance from the outflow axis). Meanwhile, the ambient mass distribution is described by the singular isothermal toroid solutions of Li & Shu (1996). Its density is the highest near the equator and decreases to nearly zero in the polar regions. Combining the wind and toroids, the unified wind model can naturally explain the dual high-velocity jet and low-velocity shell components commonly observed in the youngest group of molecular outflows associated with the Class 0 protostars.

Using 2D axisymmetric simulations, Shang et al. (2006) shows that the dense axial region of the primary wind can form a collimated jet-like structure along the axis, and the diverging part of the wide-angle wind will gradually sweep up the ambient mass to form a cavity-like structure near systemic velocity. The follow-up work of (Wang et al., 2015) decouples the wind temperature from the ambient by using a wind tracer field. This two-temperature scheme allows the study of interplay between wind magnetization and thermal pressure.

To more realistically compare the simulation and observations, we construct synthetic maps and PV diagrams for the simulations. In estimating the CO $J=$ 2 – 1 emission lines, we calculate the level population of different rotational $J$ states assuming statistical equilibrium between processes of spontaneous emission and collision (de-)excitation with molecular hydrogen. For simplicity, effects of radiative excitation and de-excitation on the level populations were ignored so the line emissions can be calculated locally. Such effects may be more important for line ratio studies, but should have little impact on the gas kinematics we are interested in here. The main difference between a statistical equilibrium result and a thermal population is found at lower density regimes. When collisional processes are not frequent enough to thermalize the population, the lower rotational states will be relatively more occupied. To obtain image cubes, we perform ray tracing along line-of-sight direction and assume a line width of 1 km s*-1* when accounting for the gas velocity.

The primary wind in the unified wind model actually bears similarities with the picture illustrated in Figure 4. Its wind velocity field is also directed along spherical-$r$ direction and has an angle-independent magnitude. Its density profile is cylindrically stratified and is highly concentrated toward the axis, which effectively makes the wind only visible close to the axis, resembling the finite size (opening angle) of the ejecta. Because of the similar setup, it is quite straightforward to apply the simulations to our current problem. Our goal is to reproduce the observed kinematic feature of the I04166 jet, so we will only focus on the high-velocity feature from the axially concentrated part of the primary wind. A few modifications are adopted to better fit the observations of I04166.

The observed knot width of the I04166 jet is better reproduced when we include an inner cone for the wind. Within $3$° half-opening angle, the rising density and toroidal magnetic field strength is capped at constant values. The gas behaves and interacts differently according to the toroidal magnetic field strength. The tooth-like velocity pattern is better maintained for weaker magnetization. We adopt a large Alfvén Mach number of $M_{A}=150$ for purposes of better modeling. The gas kinematic temperature is assumed to be $100$ K.

3.2 Modeling the Ejection History

An ejection history needs to be designed to reproduce the observed velocity pattern. In the model framework, this corresponds to specifying a velocity and density history for the wind. We adopt an inclination angle of $52$° (or $38$° from the plane of the sky), and make a few assumptions to simplify this task. The first is that we fix the wind density throughout a simulation and only vary the wind velocity. And the second is that we will try to follow a periodicity of $\sim 93$ year, which is estimated from the roughly equal spacing of the first few knots of I04166 jet. The dynamical age of the fitted regions (up to knots B7 and R7) is about 1200 years.

In the simplest of possible scenarios, we can model each tooth pattern using a single burst. This also most closely resembles the picture discussed in the previous section. Periodic bursts in protostars may occur through episodic accretion, which could arise, for example, in a gravitationally unstable circumstellar disk disturbed by binary companions. Evidence of submillimeter variability has been reported for the Class I protostar EC 53 in Serpens Main (Hodapp et al., 2012; Yoo et al., 2017). Theoretically, material can also build up in a disk before it rapidly accretes and results in a burst. Details of the processes are less well understood especially in the protostellar phase (Hartmann et al., 2016). Regardless of the origins, as the ejections come from physically small regions close to the star, the instantaneous perturbation to the velocity can be considered constant from the view of modeling purposes.

Figure 6 shows the model velocity history in which $t=0$ represents the time when the outflow is observed. The thicker dark line segments mark the time and speed of the bursts. They are separated by $\sim 93$ years in time, each having a constant speed and lasts for $20$ years long. The velocities of the ejecta are designed so as to let them travel ballistically to the observed location by the time of $t=0$. Between the bursts, the velocity is zero so no material is ejected during the break. Note that the simulations are intended to explore the feasibility of the concentrated spherical wind geometry, so the velocity histories shown here fits to the history of I04166 outflow. For example, we have actually skipped one cycle in both the blue and red-shifted lobes to better match the observed pattern.

Figure 7 shows the resulting synthetic PV diagrams along the outflow axis. They clearly present a series of linear features similar to the analytic results in Figure 4. The slope of the pattern varies with offset of distance, where the inner ones are apparently steeper than the outer ones. The simulation therefore show that the velocity pattern of spherical wind-like ejecta is indeed like a sawtooth. For a realistic comparison with observation, we also convolve the synthetic image cube with a Gaussian beam of $3$″  full width at half maximum to mimic the resolution of the CO $J=$ 2 – 1 data (Santiago-García et al., 2009). The lower two panels compares the convolved data (dark contours) with observations (dotted color contours). Although not perfectly, our numerical simulations successfully reproduce the observed velocity pattern.

While things seem to work out of the box, some tweaks on the velocity history are actually needed to achieve this desired outcome. The main problem is that the emission of an ejecta will gradually become weaker as it travels away from the origin. This is because its structure is becoming more diffused due to both the expansion and the effect of gas pressure. To achieve a comparable emission to the observations, we have put multiple ejecta together to form the knots beyond $\sim 7000$ AU. We do this by arranging a $\sim 3-5$ km s*-1* difference in intrinsic velocity between successive periods so that they end up running into each others to form a more massive feature. This is, of course, not the only way of achieving this. Directly raising the wind density for the particular knots may also work. It is probably not easy to tell which approach is closer to reality with the current data.

While the single burst design can successfully produce the main velocity structure, the generated tooth patterns (Figure 7) appear to be too isolated from each other compared to the observations. This may suggest that the true velocity history is more continuous than just a series of distinct bursts. As knots are commonly formed through the catching up process in a variable velocity flow, we now try to model the same pattern with a more continuous velocity history that consists of many low–high velocity cycles similar to a square wave. The revised velocity history is shown with blue and red line segments in Figure 6. A $40$ year higher velocity episode now follows a $40$ year lower velocity episode preceding it within each $\sim 93$ year period. Nothing is launched during the remaining gaps. The mean velocity of each cycle is similar to the velocity in the single burst case, and the difference between high and low velocities corresponds to a wave height of $0.03$. In this design, the higher velocity material is expected to run into previous lower velocity material, generating a denser layer in between.

The PV diagrams from the low-high cycle approach is presented in Figure 8. The new velocity history produces sawtooth pattern at large distance just like the single burst model, but the pattern is more complex in the inner region. Although the materials ejected during the low and high velocity episodes of a cycle are initially distinct, they gradually merge together and move ballistically as one single packet. In the current case, they appear as a single feature after a few thousand AU projected distance from the origin.

After convolution, the detailed velocity history are no longer discernible in the PV diagrams, but the differences in velocity history do affect the appearance of the tooth pattern especially for the inner ones (the lower two panels of Figure 8). The first difference is that the tooth now appear slightly wider. This is because the outflow mass is now more continuously distributed along the outflow axis instead of being only confined in discrete ejecta right from the beginning. The second difference is that the apparent slope of the tooth pattern is modified. This is because we are now not just observing the denser layer of mass. Emissions from the leading low-velocity gas and the trailing high-velocity gas also contribute to the velocity pattern, and they modify the tooth in a way that makes it look less steep.

Comparing the results in Figure 7 and 8, we see that the more evolved outer knots are indeed well described by the model of thin spherical shell. But for the younger knots closer to the driving source, one needs to worry about the continuous velocity variation as the gas has not yet evolved or collided into a thin layer. Interestingly, this effect may potentially explain the systematically smaller dimensionless slope magnitude $|\hat{\eta}|$ of the inner knots compared to the more evolved outer ones (Figure 5).

Finally, we check how the synthetic channel maps compare to the observations for completeness (Figure 9). Four channels centering at $\pm 37.5$ and $\pm 42.5$ km s*-1* relative to the system velocity are presented to compare the knot morphology of the EHV gas. The simulated knot positions generally agree with the observation and the strength are also similar. The most apparent differences are found in the size of the knot especially in the blue shifted lobe. For the outermost blue knot (B7) in the top panel and the second outermost blue knot (B6) around $\sim 7900$ AU in the second panel, the observational data are clearly more extended than the simulation result. For the knot B7, we did not really attempt to reproduce its emission because its velocity pattern along the axis is quite fragmented in the first place. It is therefore not surprising that the synthetic channel maps do not match well there. For the knot B6, it turns out that its relatively short tooth feature on the PV diagram we have successfully reproduced is misleading about its lateral size. The channel map shows that it is actually quite extended in the lateral direction. The asymmetric knot is not well represented by the current axisymmetric model.

Numerical simulations in framework of the unified wind model can successfully reproduce the sawtooth-like velocity pattern observed in the I04166 molecular jet. This supports the proposed spherical wind geometry interpretation of the pattern.

4 Discussion

4.1 Wind or Jet

Historically, models for outflows fall into two main paradigms: of wide angle wind and collimated jet. The former can explain the wide opening of outflow cavities while the latter fits the behavior of more collimated systems. In Class 0 outflows, both jet and shell components are simultaneous present and the unified wind model provides a more natural explanation. In this picture, the outflow is a spherically directed wide angle wind and the central jet-like feature arise from the high density concentration near the axis. The I04166 molecular outflow has been found to present low and high-velocity emissions that form spatially distinct components. The much wider opening of the low-velocity shell suggests it is not driven by the collimated central jet component, and probably represents swept-up material from some underlying wide angle wind that is not directly visible (Santiago-García et al., 2009). The analysis of this work shows that the ordered sawtooth velocity pattern seen in the jet component can actually represent a spherical wind-like velocity field, which further reveal the wind nature of the source. The observations of the I04166 outflow is therefore in support of the picture outlined in the unified wind model of Shang et al. (2006), in which a magnetocetrifugal wind is launched into an ambient toroid.

While a sawtooth velocity pattern is well explained by this wind picture, an authentic jet paradigm can actually also generating a similar pattern. The key to have lateral velocity component from a parallel jet is to have post-shock gas ejected from an internal working surface arising from a variable velocity flow. Santiago-García et al. (2009) has pointed out that the numerical simulations of pulsed jets in Stone & Norman (1993) can result in a sawtooth-like velocity features similar to that observed in the I04166 jet. The kinematics of the two extended knots R6 and B6 (near $\sim 7900$ AU offsets from the source) mapped by ALMA in are also explained by the geometric model proposed in Tafalla et al. (2017). The authors found that, after subtracting out the mean motion of the moving knots, the residual velocity field can be well described by a slightly curved expanding disk with an expansion speed increasing linearly from [math] km s*-1* at the center to $13$ km s*-1* at the edge, and such expansion velocity may be potentially attributed to side-way ejection from an internal working surface.

It turns out that the kinematics of the model proposed in this work can appear quite similar to the model of Tafalla et al. (2017). This becomes apparent if we also switch our reference frame to one where the tip of the blob is at rest. Consider an arbitrary part of the ejectum that is at an angular offset of $\alpha$ from the outflow axis. Its distance from the center along the surface is proportional to $\alpha$, and its perpendicular velocity component is proportional to $\sin\alpha$. For a small angle, $\sin\alpha$ is approximately equal to $\alpha$, which means that the perpendicular (or “side-way”) velocity component will increase linearly with distance from the axis, similar to the picture of an expanding disk discussed in Tafalla et al. (2017). To explicitly show that the high resolution ALMA data of knots B6 and R6 can be explained by both the expansion disk scenario and our spherical wind scenario, we construct a simple 2D axisymmetric model for knot B6 and R6 with the meridianal velocity distributions specified by a spherical wind. The model density and velocity distributions are shown in Figure 10, and the corresponding synthetic images are presented in Figure 11.

The upper and lower rows of Figure 10 are models for knots B6 and R6, respectively. We use exactly the same velocity field and they differ only in the knot position from the origin. We mimic the setup of the numerical simulations presented in Section 3 so the wind is purely in the spherical-$r$ direction and of a constant magnitude. The density profile drops off with distance from the axis as $\rho\propto\varpi^{-2}$, but is constant inside an inner cone of $3$° half opening angle. To model the density distribution of a thin mass layer, a Gaussian profile of FWHM $490$ AU ($3.5$″) centering at a fixed distance from the driving source at the origin is applied along spherical-$r$ direction. The magnitude of the perpendicular velocity component is determined by the distance of the ejecta from the driving source (and the adopted inclination angle). Its magnitude increases roughly linearly from 0 to $\sim 15$ km s*-1* from the axis to the outer edge of the white contour as shown in the second columns of the figure.

Figure 11 compares the synthesized model channel maps and velocity maps for the knot R6 and B6 to the observations. The knot positions within each channel and the characteristic velocity gradient along the axis are well reproduced. The main difference is that the observed emissions appear to be more fragmented than the model, especially in the blue lobe. The observations also appear to be relatively more curved than the model emission, which may suggest that the ejecta are slightly more curved than those from spherical shell.

While both scenarios can similarly explain a tooth-like velocity pattern (because of their approximately equivalent apparent kinematics), the spherical wind model is more constraining especially when applied to a series of knots. It predicts that the dimensionless slope $\hat{\eta}$ of all tooth features should all be the same and only depend on the outflow inclination angle (Equation 3). As illustrated in Figure 4, when the projected velocities and offsets are both fixed in a given observation, the outflow inclination angle is the only free parameter that can be varied to fit the sawtooth velocity pattern. The reasonable agreement across the entire sawtooth pattern is therefore a compelling sign of an underlying wind-like geometry. On the other hand, in the paradigm of parallel jets with internal working surfaces, the slope of the tooth feature should depend critically on the strength of the side-way ejection velocity (besides the jet inclination angle) because this is the cause of the apparent velocity gradient on the PV diagram. As the shock conditions may simply differ for different knots, there is no guarantee that the series of teeth pattern should behave coherently, and the observed systematic variation of slope with distance is not automatically explained in this paradigm. In other words, the origin of lateral velocity components in a spherical wind scenario is purely geometrical, so it is not possible to vary the predicted slope of some of the teeth while keeping others unchanged. The systematic and coherent behavior of the teeth slopes in the I04166 jet is therefore a diagnostic of a spherical wind-like velocity field, although this does not rule out the possibility of a jet shock. We note that a thermal sound speed of $\sim 15$ km s*-1* corresponds to a thermal temperature of $>60000$ K. A jet-shock scenario would therefore predict very high temperature at least for the freshly-shocked inner knots before heavy cooling can take place.

4.2 Implications on Velocity History

A correct geometric picture is crucial for inferring the intrinsic outflow velocities from the observed line-of-sight components, which is fundamental for the study of outflow velocity history. For example, observations of the I04166 jet reveal an apparent velocity variation on order of $\sim 5$–$10$ km s*-1* within each emission peak. Without any prior knowledge, one might simply de-project the velocity using a fixed inclination angle and conclude that the flow velocity varies periodically with a wave height of $\sim 8$–$16$ km s*-1*. The catch of this seemingly natural conclusion is that it implicitly assumes a geometric picture of unidirectional flow with varying speed along the axis, which is probably not true. The modeling of Tafalla et al. (2017) showed that a tooth pattern likely results from the lateral velocity of the flow, and in this work we show that the entire pattern can be consistently explained by a spherical wind like geometry. This leads to a very different conclusion where the observed $\sim 5$–$10$ km s*-1* apparent velocity variation in each peak may have little to do with outflow velocity variation. As demonstrated in the simulation of the single-burst velocity history model (Figure 7), ejecta of constant speeds can reproduce the observed feature well.

A “slow head–fast tail” velocity pattern on a PV diagram could in principal form in three ways. The first is a diverging velocity field like the spherical wind discussed in this work. The second is the side-way ejected gas from internal working surfaces of shocks. The third is an intrinsic velocity variation along the flow direction, which is also the true velocity history of general interests. In the particular case of I04166 outflow, the series of clearly detected tooth-like pattern have allowed us to reveal its wind nature. In a more general case without multiple clearly-detected teeth spanning a range of distance from the origin, however, discriminating between the potential origins is difficult. The three mechanisms are also not mutually exclusive, and could be all at work in a same source. A long as the flow velocity field is diverging in nature, there will be a corresponding dispersion in the observed line-of-sight velocity component due to projection. As long as there is compressed, high density region, lateral velocity could be generated from the pressure gradient force like in the case of an internal working surface, and the real question is whether it is significant or negligible. And as long as the outflow launching velocity is time dependent, there will be velocity gradient along the flow direction, whether it is an unidirectional jet or a diverging wind. Disentangling contributions from these different mechanisms is therefore required for an accurate understanding.

The velocity record is better preserved in regions closer to the driving source before gas traveling at different speeds collides with each other and merges. For example, the resulting synthetic PVs of the model of single burst history in Figure 7 and the model of more continuous cycle history in Figure 8 are notably different within a few thousand AU but appear almost identical beyond that. High angular resolution observations toward the inner region of the outflow is needed to assess the velocity history. For example, Wang et al. (2014) mapped the inner few thousand AU of the I04166 jet in CO $J=$ 3 – 2 using the Submillimeter Array at $\sim 1$″  resolution. The PV diagrams are reploted in Figure 12, which show finer structures especially in the blue-shifted lobe.

The velocity feature between $\sim 1500$–$2000$ AU in the blue-shifted lobe is particularly interesting. Unlike the main tooth pattern that are likely associated with a spherical wind-like geometry, the emission between them appears to show a reverse pattern. It may be possible that it represents the true velocity variation over time. In this case, the increase of velocity with distance may be interpreted as a decrease of outflow launching velocity with time over the particular era. Unfortunately, the emission is weak and fragmented and would need to be confirmed by higher sensitivity observation before further modeling could help understand its nature. The innermost red-shifted emission is also interesting. There seems to be multiple overlapping teeth pattern within 1000 AU from the source and it is possible that we are looking at a series of ejecta launched with increasingly higher speed. Still, higher sensitivity data is needed to confirm its behavior.

The fact that the outer knots do not preserve clues of the detailed velocity variation is actually a benefit when testing whether the pattern is consistent with a spherical wind scenario or not. After gas of different velocities merges into a single packet and move more-or-less at a constant speed, the observed sawtooth patterns will most likely be associated with lateral velocity components and we can worry less about the velocity variation along the flow direction. The model predictions are hence more reliable for the other knots. Establishing the wind nature of the velocity field can also help interpret the detailed velocity structures found close to the source. Without such knowledge, it would be difficult to tell whether the tooth-like patterns in the inner knots in Figure 12 are signs of a wind-like velocity field or an axial velocity variation. In this sense, the sawtooth velocity pattern of the I04166 jet poses an unique clue for the study of its ejection history.

Similar quasi-linear velocity patterns have also been seen in some other Class 0 outflows. SMA SiO $J=$ 8 – 7 and CO $J=$ 3 – 2 observations toward the L1448C outflow revealed tooth-like velocity patterns in at least two pairs of its EHV peaks (see Figure 7 of Hirano et al., 2010). These “fast tail–slow head” features might have the same geometric origin as the I04166 outflow, but it would be difficult to confirm. Unlike I04166, we do not have multiple knots distributed across a large distance. Several quasi-linear features have also been identified in the HH 211 molecular jet (see Figure 7 of Lee et al., 2009). In particular, the SiO $J=$ 8 – 7 emission of the blue-shifted knots BK2 and BK3 present a “fast tail–slow head” velocity structure that were interpreted as the intrinsic velocity variation along the jet axis during a catching-up process (Jhan & Lee, 2016). It would be interesting to see whether, or how much, perpendicular velocity component (either due to geometry or from sideway ejection of an internal working surface) contributes to the observed velocity pattern. The interpretation may be further complicated by the wiggling of HH 211 jet (Moraghan et al., 2016). In any case, disentangling contributions from different potential mechanisms is critical for a thorough understanding, and would require high quality observations toward many different outflows.

5 CONCLUSION

We propose that the observed sawtooth-like velocity pattern of I04166 molecular jet can be associated with an expanding spherical wind geometry of density concentrating near the axis, similar to that predicted by the X-wind model. We derive an analytical formula for the pattern slope, which describes its dependence on the outflow inclination angle, the ejecta velocity, and its distance from the driving center. The formula can reasonably explain the systematic variation of slopes across the whole pattern using one same underlying outflow inclination angle, which is not automatically explained in the paradigm of jet and internal working surface. The reasonable agreement is therefore a quite compelling signpost of an underlying spherical wind-like velocity field in the I04166 jet. The dimensionless slope of the tooth-like velocity pattern can also be used to estimate the outflow inclination in this picture.

Using numerical simulations of variable velocity wind in framework of the unified wind model of Shang et al. (2006), we are able to generate synthetic CO $J=$ 2 – 1 PV diagrams that reproduce the observed pattern well. This realistically demonstrates the feasibility of the proposed wind scenario. While a spherical ejecta picture can nicely describe the kinematics of the outer knots, clues of the velocity history of general interest are better preserved in the inner knots where gas of different speeds has not yet completely merge together. Different potential mechanisms may underlie the observed velocity pattern, including the intrinsic velocity variation along the flow direction, the potential side-way velocity from internal working surfaces of shocks, and the diverging (wind-like) geometry of the flow. A correct model of flow geometry is therefore critical for an accurate understanding of the flow history. The unique sawtooth pattern of the I04166 jet has allowed us to associate it with a wind interpretation, which could help disentangle the velocity pattern in the innermost region. High resolution high sensitivity observations could help confirm the detailed flow structure and shed light on its velocity history.

We thank M. Tafalla and J. Santiago-García for providing the CO $J=$ 2 – 1 data of the molecular outflows IRAS 04166$+$2706. We also thank an anonymous referee whose comments have helped improve the clarity and presentation of the manuscript. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00304.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work also uses data from The Submillimeter Array. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. This work utilized tools developed by the CHARMS group and computing resources and cluster in TIARA at Institute of Astronomy and Astrophysics in Academia Sinica (ASIAA). We also acknowledge grant support from Ministry of Science and Technology in Taiwan by MOST 102-2119-M-001-008-MY3 and MOST 105-2119-M-001-037-.