text:
# Simulations of Evolving or Outbursting Molecular Protostellar Jets

## 1 Introduction

Jets are being increasingly associated with protostars and with bipolar molecular outflows. Prominent examples of jet-driven collimated outflows include HH 111 (Nagar et al., 1997; Hartigan et al., 2001), HH 211 (Gueth & Guilloteau, 1999; Chandler & Richer, 2001), HH 212 (Lee et al., 2000) and HH 288 (Gueth et al., 2001). For many other outflows, there is evidence that their momentum is also derived from the thrust of jets. For example, transverse shock structures behind bow shocks are interpreted as the sites of jet impact (Smith et al., 2003). Nevertheless, estimated jet thrusts appear insufficient to drive some outflows (Fridlund & Liseau, 1998) and ballistic bullets and wide winds have been invoked. Alternatively, as investigated here, the presently observed jet may not reflect the ancestor jet that has previously driven the outflow.

Our goal is to determine the observable characteristics of outflows driven by molecular jets with evolving mass flow rates. Jet evolution was considered as an interpretation of the systematic decrease in flow velocity with distance in the atomic HH 34 outflow (Cabrit & Raga, 2000). The deduced requirements appeared unlikely and an alternative precessing-fragmenting scenario was favoured. A systematically increasing jet speed was simulated by Lim et al. (2002). Their idea was to explain how molecules can be accelerated to high speed without destruction. Here, we choose to evolve the mass flow rate rather than the speed since it is clear that the mass flow rate must change by many orders of magnitude as a protostar evolves. The taxonomy of protostars is defined by the categories Class 0 through Class 3, which have been linked to an evolutionary sequence (Lada, 1987; Andre et al., 1993). One scenario connects this sequence to simultaneous evolutions in both the accretion rate and the outflow rate (Smith, 2000). Furthermore, certain protostars undergo outbursts (Narayanan & Walker, 1996), the consequences of which may also show up later in the outflow properties (Arce & Goodman, 2001).

We investigate a limited set of three-dimensional simulations that correspond to significant and permanent changes to the input mass flux. This complements studies of outflow development generated by non-evolving heavy molecular jets (Suttner et al., 1997; Davis et al., 1998; Völker et al., 1999) as well as a range of jet conditions (Rosen & Smith, 2003). We first consider a rapid decrease in the mass flux, which we suspected would produce a bullet-like outflow. A decrease in mass flux of some variety must accompany the transitional stages of a protostar into a star. In addition, we model a jet that increases density after some initial propagation, and this model correlates with the opposite end of a protostar’s evolution. During the initial collapse of the core into a protostellar disk, a very young jet must ignite and, therefore, be characterised by an increase in mass flux.

Certain features of the molecular line emission from outflows have been attributed to an evolving jet mass outflow rate. For example, it is believed that the energy radiated from shock waves provides a measure of the present driving power while the mechanical luminosity provides a record of the driving history. In this way, the ratio of emissions from vibrationally excited molecular hydrogen and rotationally excited carbon monoxide could determine the evolutionary phase of the protostellar exhausts (Davis et al., 1998; Yu et al., 2000). We can determine here the validity of this statement, although these simulations are of a particular type and scale and so only begin to illustrate the full range of possibilities. For example, we take here a hydrodynamic highly-collimated jet and a uniform ambient medium that are both initially fully molecular. In the next section, we discuss more fully the characteristics of the jet simulations.

## 2 ZEUS-3D with Molecular Cooling and Chemistry

### 2.1 The code

We have modified the ZEUS-3D code, which usually updates variables with an explicit method, to include a semi-implicit method that calculates the molecular and atomic hydrogen fractions, according to the prescription of Suttner et al. (1997). We have also added a limited equilibrium C and O chemistry to calculate the CO, OH and H2O abundances (Rosen & Smith, 2003). Equilibrium CO and H2O is a reasonable estimate at high density, consistent with the low numerical resolution of the cooling layers behind shock waves (Rosen & Smith, 2003).

The details of the many components of the cooling function and chemistry network are discussed in the appendices of Smith & Rosen (2003). As an overview, we include cooling through rotational and vibrational transitions of H2, CO, and H2O, H2 dissociative cooling and reformation heating, gas-grain cooling/heating, and a time-independent atomic cooling function that includes non-equilibrium effects. We include an additional 10% by number of helium atoms, so the mean particle mass is 2.32 $`\times `$ 10-24 g. The dust temperature is fixed at 20 K.

The ZEUS code suffers from numerical errors. Recently, Falle (2002) has explored the appearance of rarefaction shocks which, as the name suggests, are particularly resistant and their influence cannot be removed by increasing the resolution. They may occur in locations where the local Courant number is near 0.5 and where the velocity is negative. Inaccurate fluid patterns could result for adiabatic or MHD flows but isothermal flows are not significantly in error. The present jet simulations are thus immune to these errors, with strong forward motion and strong cooling, especially in the zones with relatively large Courant number (i.e. approaching 0.5). In any case, isothermal flows are an approximation in which regions undergoing rarefaction are being artificially heated (to counter expansion cooling). Cooling also generates structure through thermal instabilities on many scales which are not adequately featured in the numerical approximation. Nevertheless, we persist with the ZEUS code because of its high versatility, allowing the introduction of new physics, the discovery and removal of specific errors over many years, and its speed of execution.

### 2.2 The fixed conditions

Two models for the time evolution of the jet are studied. We designate these as Simulations D and I, with a decreasing and increasing jet density, respectively.

In both 3D simulations, the computational grid spans 2 $`\times `$ 1017 cm along the jet axis ($`x`$) and 2 $`\times `$ 1016 cm in both transverse dimensions, with 1000 $`\times `$ 100 $`\times `$ 100 uniform zones in each direction. We initialise the jet with a nozzle radius, Rj, of 1.7 $`\times `$ 1015 cm or 7.5 zones. The initial internal Mach number of the jets is 140 with a temperature for the fully molecular jets of 100 K; this combination gives a speed of 100 km s-1. The ambient medium is molecular with a hydrogen nucleus density of 104 cm-3.

In addition, both jets are pulsed with an amplitude of $`\pm `$ 30% and a period of 60 yr, with a $`\mathrm{sin}(\omega t)`$ dependence. A small precession with a 1 half-precession angle and a period of 50 years is included. A radial shear that reduces the speed to 70% of the core velocity at the jet edge is also included. The sinusoidal form of the pulsation and the parabolic form of the shear follow Eqns. 1 and 2 of Völker et al. (1999).

### 2.3 The evolving parameters

Two models for the time evolution of the jet are studied. We designate these as Simulations D and I, with a decreasing and increasing jet density, respectively.

In Simulation D, we decrease the jet density abruptly, setting it to zero after a short period of time. We resorted to this dramatic variation of the flow because earlier attempts with more moderate decreases failed to show a significant deceleration. The end of the run may thus correspond to a highly evolved jet-driven outflow. Specifically, we allow the “jet” to propagate for a mere 20 years and then sever its ties to the life-giving momentum source at the inlet. The jet is initialised with a hydrogenic nucleon density of 105 cm-3 ($`\mathrm{log}n_j`$ = 5). In addition we have simulated other aborted jets, with flows terminated after 20 years but the flow was not initially precessed, and also a case where the flow was terminated after only 5 years and also not initially precessed.

In Simulation I, we wish to show the consequences of increasing the mass density in a light jet flow that had already propagated some distance. In this case, the density evolution follows a less abrupt change, determined by the following:

$$\rho _j(t)=\rho _0+\rho _1\mathrm{tanh}\left(\frac{tt_1}{t_0}\right)$$
(1)
where we set $`\rho _1`$ = 1.1484 $`\times `$ 10-19 g cm-3, $`\rho _0`$ = 1.1716 $`\times `$ 10-19 g cm-3, $`t_1`$ = 800 years, and $`t_0`$ = 50 years. This results in the logarithm of the jet hydrogenic number density rising from $`\mathrm{log}n_j`$ = 3 to $`\mathrm{log}n_j`$ = 5 over a 200 year period centered at a program time 800 years after the simulation has started. Given the axial grid length of 2 $`\times `$ 1017 cm, this leads to the bow shock from the initially light jet being roughly half way across the grid when the mass flux is increased.

## 3 Internal Structure

Slices of jet density down the midplane are presented in Figure 1. Both simulations feature a refocused leading region and wider trailing “shoulders”. We remark that the shoulder regions in both simulations shown in Fig. 1 have progressed to only 7 $`\times `$ 1016 cm, a small fraction of the jet length. In our study of non-evolving jets (Rosen & Smith, 2003), we were able to attribute the refocusing to the combination of both shear and strong molecular cooling, since simulations with atomic jets or molecular jets with no shear do not refocus. The ratio of the distance of the shoulders from the nozzle relative to that of the entire jet was found to be proportional to the jet-to-ambient density ratio. Hence, without a continuous strong power source in Simulation D, the bloated cocoon is stunted.

Furthermore, the average speed of the leading bow shock, which almost traverses the grid in 1300–1400 years, is $``$ 45 km s-1. This is similar to a “light” jet simulation (jet-to-ambient density ratio of 0.1 in Rosen & Smith 2003). This average speed is still, however, nearly 90% faster than expected from ram pressure arguments due to the concentration of the jet thrust caused by the refocusing.

From a plot of a time-averaged advance speed of the bow shock (Figure 3), we see that the advance speed of the bow shock is not constant throughout either simulation. In Run D, the advance speed of the jet falls below 60 km s-1 within the first 100 years, and follows a linear decrease in velocity with time until t $``$ 800 years, where it continues a linear decrease but at a faster rate. The change of the deceleration rate is related to a broadening of the tip of the bow shock, which is at the leading edge of a region of refocused and accelerated flow after $`t`$ 200 years. At the end of Simulation D the bow shock is progressing at a speed of $``$ 30 km s-1. This is an example of the ability of all of our aborted molecular jet simulations to maintain its momentum for the length of the computational grid (0.06 pc). Even in the case of the jet aborted after 5 years, we still found a speed of $``$ 20–25 km s-1 after the 2000 year interval in which it crossed the grid.

In contrast, the progress of the bow shock in Simulation I, as measured from the molecular emission maps presented below, increases early in the simulation. The speed initially drops to 20 km s-1, but after only 50 years begins to rise. There is a significant increase at $`t`$ 300 years, after which the speed of advance of the bow shock is $``$ 50 km s-1. The increase does not coincide with the onset of the higher jet density but, instead, coincides with the initial refocusing of the bow shock at its leading edge. Thus, the stripping of the low-speed jet sheath is an interesting means by which an acceleration in the proper motions of jet knots may result. Even after the jet density has increased, the average advance speed for the bow shock remains nearly constant until the simulation is terminated.

The dense material in Simulation D is confined to the bow shock and swept-up shell while dense material is also associated with the internal shocks of the pulsed jet in Simulation I. The cocoon surrounding the jet behind the shoulder region has the lowest density, with the still-connected jet in Simulation I having a lower minimum density by a couple of orders of magnitude than the aborted jet in Run D. These overdense and underdense regions within the jet are also apparent in the axial cross sections of density in Figure 2. With the ability of the jet in Run I to generate a very low density cocoon, the relatively small overdensity (still an order of magnitude) in the bow shock is difficult to see. Note twin concentric shell structures are formed in the cross section of Run I, visible in Fig. 2 at distances Rj = 54–91, where the powerful jet transmits invigorated shocks through the cocoon.

## 4 Simulated Molecular Line Emission

### 4.1 The Spatial Distributions

Images of line emission from molecules are shown for Simulation D in Fig. 4 and for Simulation I in Fig. 5. Integration is across the $`y`$-direction transverse to the jet axis, and the simulation time corresponds to the previous figures. The H2 line emission is based on a non-LTE approximation to the vibrational populations (Hollenbach & McKee, 1979), and the CO emission is calculated from the formulae of McKee et al. (1982). In addition, we set a lower velocity limit of 1 km s-1 on the integration of CO emission lines in order to subtract out the undisturbed cloud.

In Figure 4, we show the H2 1–0 S(1) line, the CO R(1) and R(5) lines for the CO ground vibrational state ($`\upsilon `$ = 0). In the H2 line, only the leading edge of the bow shock remains visible and, since Simulation D is relatively devoid of warm (T $`>`$ 5000 K) gas, the image in the H2 2–1 S(1) line is not shown. We display instead the CO R(5) line, which should highlight gas at intermediate temperatures (T $``$ 100 K) between the H2 line and the CO R(1) line (which emphasizes the coldest, T $``$ 5 K gas). The CO R(5) line image not only has an intermediate appearance but also shows some filamentary structure, which is brightest at the leading edge of both the bow shock and the shoulders. This relatively warm gas is generated by the accelerating refocused region. In both of the simulations presented here, the moving CO R(1) shows the outline of the bow shock.

Simulation I is sufficiently energetic to generate copious amounts of 10,000 K gas. Hence, in addition to the lines shown in Fig. 4, we also show the relative luminosity of the H2 2–1 line to the 1–0 line in Fig. 5. The emission in the H2 1–0 image is sparse and filamentary, showing several ‘minibows’ especially at the edges of the shoulder region near $`x=`$ 9 $`\times `$ 1016 cm. Only a single pulse near the inlet is seen in the H2 1–0 map. Earlier pulses, farther down the jet, are easier to see in the ratio of the H2 lines, and are noticeable in the CO R(1) line, but very prominent in the map for the CO R(5) line. While the CO R(1) emission map does show a wider shoulder region well behind the leading edge of the bow, this region is much less pronounced than the appearance in the midplane density slice.

### 4.2 Totals

The evolution of the totals of the H2 emission lines and the CO R(1) line are shown in Fig. 6. In Simulation D, a dramatic reduction in the H2 totals begins after $`t`$ = 100 yr. A flattening of the CO luminosity, which monotonically increases for non-evolving molecular jets, also occurs but not until $`t`$ = 700 yr. Thus, a significant change in the ratio of the H2/CO total emission occurs shortly after the flow has been turned off. It is only when the H2 1–0/CO R(1) ratio approaches unity that the total emission from each emission line tends to flatten out. Before this occurs the ratio is a good age indicator for the aborted jet flow.

The 80 year delay between the jet demise and the H2 weakening is clearly the time to elapse before the leading bow stops being driven by the jet impact. The 700 year delay between the termination of the jet flow and the flattening of CO R(1) emission indicates that new ambient material continues to be snowploughed, as the momentum of a small fraction of high-speed gas is redistributed into a large fraction of low-speed gas. This is a true ‘bullet’ or ballistic phase. To test this, the deceleration of the bow shock in Run D as shown in Fig. 3 can be compared to the analytic motion of a bullet of speed $`v_b`$, fixed mass M, cross-sectional area A, and initial density $`\rho _{bo}`$, moving into a medium of density $`\rho _a`$, with a drag coefficient Dc. The equation of motion is

$$\frac{dv_b}{dt}=\frac{v_b^2}{L}$$
(2)
where $`L=M/(\mathrm{D}_c\rho _aA)`$. Here L defines an effective column length of ambient gas required to slow down the bullet. First, we take ram pressure confinement of a bullet which maintains its shape as the pressure decreases. This yields

$$L=L_o\left[\frac{\rho _{bo}}{\rho _a}\right]^{2/3}\left[\frac{v_b}{c_s}\right]^{4/3},$$
(3)
where $`c_s`$ is the bullet sound speed and $`L_o`$ is the initial value of L. Substitution and integration gives

$$t=3\frac{L_o}{v_{bo}}\left[1\left(\frac{v_b}{v_{bo}}\right)^{1/3}\right]$$
(4)
and

$$x=\frac{3}{4}L_o\left[1\left(\frac{v_b}{v_{bo}}\right)^{4/3}\right].$$
(5)
Alternatively, if the bullet is either disk shaped or does not expand at all, L is a constant and we find

$$\frac{v_b}{v_{bo}}=\frac{1}{1+t(v_{bo}/L)}.$$
(6)
For the particular conditions of Run D, $`L_o=(6\times 10^{16}/D_c)`$ cm. To decelerate from 60 km s-1 to 30 km s-1 in 1200 yr, then requires Dc = 0.17–0.28, depending on the bullet model utilised. This drag coefficient is then consistent with the aerodynamic shape of the leading bow and thus shows that a large mass can be mildly disturbed by a high-speed jet.

The duration of this phase is 35 times longer than the jet-driving phase, suggesting an effective jet column density which exceeds the deflected ambient column density by 35. We can roughly attribute this to the high jet density (10) and the refocusing (3.5).

It should be remarked that decelerating bullets have been investigated in various astrophysical contexts including Herbig Haro objects (Norman & Silk, 1979), radio galaxies (Smith & Norman, 1980) and planetary nebula (Palmer et al., 1996). Depending on how the bullet forms, stability is maintained through inertial confinement for a sound-crossing time or a shock-wave crossing time. Therefore, the high Mach number of the bullet in Run D permits a bullet to travel over 1017 cm without great expansion.

In Run I, periodic brightening in H2 luminosity occurs as the pulses in the flow pass through the shoulder region and energise the bow shock envelope. The increase in jet density that becomes significant for $`t>`$ 750 yr shows up as an increase in the H2 output very soon thereafter ($`t`$ = 800 yr). No accompanying increase in the CO emission is noticeable. There is an unrelated early break in the rate of increase in the CO R(1) total luminosity, at t = 100 yr. The slightly increased rate is related to the shape of the envelope, which shows the effects of refocusing near the start of the simulation.

The final integrated fluxes of the two simulations have contrasting behaviours. In Run D, the absolute luminosity in the 1–0 S(1) line of H2 is continued at the low level of L(1–0 S(1)) $``$ 10-5 L despite the lack of a jet. This compares to the value of $`10^3`$ L for the equivalent simulation in which the jet power is maintained (Rosen & Smith, 2003). In Run I, L(1–0 S(1)) reaches 3 $`10^4`$ L and is set to increase. It is clear that a complex relationship exists between the immediate jet power and L(1–0 S(1)). Whereas we found this ratio to lie in the range 80–600 for non-evolving jets (Rosen & Smith, 2003), we here find values between zero (Run D, after 20 yr) and 1.2 $`\times `$ 105 (Run I, 800 yr).

The predicted values of L(1–0 S(1)) are comparable to those derived in the unbiased statistical study of (Stanke et al., 2002), where the majority of outflows possess L(1–0 S(1)) in the range 10-4–10-3 L. Nevertheless, the absolute fluxes are quite low in comparison to a number of observations of well-known outflows. This is partly due to the compactness of the simulated outflows, limited by our need to follow molecule dissociation. We can argue that scaling up the spatial dimensions by a factor of 10 would scale the luminosities by a factor of 100. It remains to be shown, however, that the flow patterns would be maintained.

## 5 Velocity Distribution of Mass and Molecular Line Emission

### 5.1 Mass “Spectra”

CO velocity data for observed outflows are now invariably analysed by plotting the differential mass in velocity bins, with a few simple assumptions made to compute mass from CO intensity. One may classify how this mass decreases with velocity by fitting a power law, with the slope defined by $`dm/dvv^\gamma `$. Two surveys containing many sources (Yu et al. 2000 and Ridge & Moore 2001) conclude that a broken power law is often necessary, with a shallow slope of 1.0–2.0 at low velocities and a steeper slope, in some cases $`\gamma `$ is close to 10, at larger velocities. The velocity where the break occurs is usually $``$ 10 km s-1. One interpretation of this break speed is that it is the projected component of the critical shock speed, $`v_d\mathrm{cos}\theta `$ with $`v_d`$ = 23 km s-2 (e.g., Smith 1994), above which all molecules are assumed to dissociate within the shock.

We have performed a similar analysis for the velocity distributions of mass, of CO R(1) luminosity and of H2 1–0 luminosity for each of the aborted jet simulations, including Simulation D, and for Simulation I. In the case of Simulation D, the analysis was performed at a program time of $`t`$ = 1300 yr, corresponding to other Simulation D illustrations presented here. In the increased-density jet of Simulation I, the tip of the jet has just barely reached the outer $`x`$-boundary at $`t`$ = 1400 yr. We, therefore, chose to inspect data at $`t`$ = 1300 yr. We list the results for a sample of viewing angles, defined as the angle of the jet axis out of the plane of the sky and toward the observer, in Table 1. We display the data for one viewing angle (15) in Figure 7.

Having only found low values of $`\gamma `$ in previous simulations involving molecular jets (although high values are possible for an atomic jet), we have chosen to investigate whether a reduced mass flux would give larger values for $`\gamma `$. In Simulation D and the other aborted jet simulations, however, we still find shallow dependences of the mass distribution on velocity i.e. small values for $`\gamma `$. The values for the aborted jet simulations are in general larger than their counterparts in the increasing jet-density simulation, with the values for both sets usually nearer 1.0 than 2.0. However, large slopes are computed for high velocities in the only cases where a break could be seen, which are the CO-derived values for the aborted jet simulation. These results are not significantly different than those found for non-evolving jets (Rosen & Smith, 2003).

In the previous analysis for non-evolving molecular jets, the $`\gamma `$ determined for CO was nearly uniformly larger than that for mass. For the evolving jet simulations here, the dependence of CO luminosity on velocity is also usually steeper than the one for mass, but in some cases it is only equal (e.g., Run D at a viewing angle of 30) or, for one aborted jet simulation, was lower. This occurred in a simulation of an aborted jet that propagated into an equal pressure ambient medium. The CO-derived $`\gamma `$ is lower than the mass-derived one at all viewing angles, and by as much as 0.6. However, in this simulation, the CO distribution can be better fitted with two power laws, with the higher velocity one much higher than the single linear fit to the mass distribution.

Very little dependence of $`\gamma `$ on viewing angle is found for both of these simulation. This is contrary to previous jet simulations (Smith et al. 1997, Lee et al. 2001, Downes & Cabrit 2003, and Rosen & Smith 2003), which have shown an inverse relationship between viewing angle and $`\gamma `$. However, the CO-derived slopes in the high velocity range are again an exception. This group disagrees more strongly with the previous results and, for small viewing angles, $`\gamma `$ varies in direct proportion to the viewing angle.

An excess at the highest velocities is found in Simulation I, as also found in previous studies of non-evolving molecular jets. This bump peaks at log $`v`$/km s-1 = 1.4 (for a viewing angle of 15) and is associated with the internal shocks from the pulsed jet. In the aborted jet simulations, both the internal shocks and the high velocity bump are absent.

Recent observations of H2 emission from protostellar outflows have been similarly analysed (Salas & Cruz-González, 2002). As with the CO data, the H2 flux distributions show a break in behaviour, with the break between 2–17 km s-1 for the selection of sources. A flat or slightly rising flux is found at low velocities and a fall-off in flux with $`\gamma `$ in the range between 1.8 and 2.6 at high velocities.

From Fig. 7, we see that Run I may best match these observational constraints. However, the H2 flux distribution in Run I rises one order of magnitude in brightness from the first H2 datum in the figure to where it turns over near $`v`$ = 10 km s-1. This is much larger than that found in the observations, although it is a smaller rise than found in the non-evolving molecular jets simulations. Additionally, the best fit line between $`\mathrm{log}v`$ = 1.0 and 1.3 has $`\gamma `$ = -1.6, which is in line with the observations. This could be an adequate explanation for sources viewed with their jet axis near the plane of the sky, such as HH 212. However, at larger viewing angles (Fig. 7 shows a viewing angle of 15 only), the low-velocity rise steepens. In addition, at intermediate velocities the comparison becomes less good as the viewing angle is increased. C-type shock physics may solve these problems since ambipolar diffusion would generate considerably more emission near zero radial velocity in the absence of discontinuous acceleration. The effect this has on the global properties of an outflow remain to be calculated.

### 5.2 Position-Velocity Maps

We display position-velocity maps in H2 and CO emission lines for Simulation D (Fig. 8) and Simulation I (Fig. 9). The H2 1–0 position-velocity map for Simulation D shows emission only at the leading edge of the jet. This emission is spread over 10 km s-1 and 1 $`\times `$ 1016 cm along the projected jet axis and is characterised by two peaks. Note that steady state bow shock models (see Fig. 14b of Smith et al. 2003) predict two peaks, which correspond to (1) the projected edge-on leading edge near zero radial velocity and (2) the blue-shifted emission from the warm material projected well behind the bow apex.

The pulsed jet in Simulation I shows intermittent H2 1–0 and 2–1 emission along its length, with many emission sites spreading out over a large range of velocities. The H2 2–1 map shows that the sites of bright H2 1–0 emission are also bright in the more highly excited line.

The emission in the CO lines shows a pair of “Hubble-law” like dependences for Simulation D, but disconnected from the nozzle. The outer envelope of the first increases relatively quickly with position, and is followed by a shallower increase for much of the remainder of the jet. While the emission in the CO R(1) line is nearly uniform for non-zero velocities, there is some structure in the CO R(5) position-velocity diagram. The bright spot near $`x^{}`$ = 7 $`\times `$ 1016 cm is from the advancing edge of the shoulder region seen in Fig. 4, and the peak at $`x^{}`$ = 5 $`\times `$ 1016 cm is from the trailing filamentary structure. Both of these bright features have radial speeds near zero.

For Simulation I, the CO R(5) diagram shows more structure than the diagram for CO R(1). Both show a prominent series of triangular shapes that distinguish the Hubble Law behaviour of a pulsed jet. There is evidence for a connective loop of emission between successive Hubble Law regions, until the triangles merge downstream, filling in all velocities, and the Hubble Law is then hidden. The maximum luminosity in the CO R(5) map falls in between those of the H2 1–0 and the CO R(1) maps. Naturally, all of the maxima are higher than the corresponding values in the aborted jet simulation.

### 5.3 Velocity Channel Maps

We display CO R(1) velocity channel maps for Simulation D in Fig. 10 and Simulation I in Fig. 11 with a resolution of 1 km s-1. While the CO emission from Simulation D does not span a large velocity range, the morphology does change as one moves from relatively high radial velocity (10 km s-1) to low radial velocity. At high radial velocities, only the leading edge of material is visible, and as the velocity is decreased more emission appears closer to the inlet. This is another aspect of the Hubble Law profile seen in the position-velocity map. Near zero velocity the entire envelope of the jet and cocoon of shock ambient material is represented, along with the shoulder region seen in the density cross sections and integrated CO maps. This defines a T-type morphology. Even the filamentary structure seen in the integrated emission map for CO R(5) is seen to some extent in these CO R(1) velocity channel maps.

The CO R(1) velocity channel maps for Simulation I reveal morphologies in the different velocity ranges more typical of our previous non-evolving jets (Rosen & Smith, 2003). At high radial velocities, only the internal shocks from the pulsed jet are seen and at low velocities the entire envelope of jet plus cocoon material is present. At low velocities, the gap between the nozzle and the onset of emission decreases and a V-type global morphology is identifiable. At intermediate velocities, a spinal structure is present.

The twin protuberances at the front of the jet are from the differing advance speeds associated with the extrema of the evolving jet density in Simulation I. As confirmed from a sequence of integrated H2 1–0 emission line maps (Fig. 12), a faster-moving bow shock leaves the path of the earlier slower-moving initial bow shock. This departure occurs at a fairly late time in the evolution of the jet ($`t`$ = 1300 yr, or 550 yr after the onset of the increase in jet density at the inlet) and well into the refocused leading part of the bow ($`x`$ = 1.6 $`\times `$ 1017 cm). The average speed of the new bow is therefore close to 100 km s-1. An average speed for the propagation of the new bow shock equal to that of the (time-averaged) jet speed is not unreasonable in this simulation with the very light initial jet initially evacuating the path and with an acceleration from the refocused region.

## 6 Application: an evolutionary diagnostic?

The structures and features detailed above can provide clues to distinguish the evolutionary phase of a specific outflow. Statistically, however, the integrated properties can provide the most employable measures. In particular, the ratio L(H2 1–0 S(1))/L(CO R(1)) can be determined from simulations and compared to the observationally derived L(H2)/Lmech. The total molecular hydrogen emission L(H2) can be estimated from L(H2 1–0 S(1)) (Smith, 1995), provided the near-infrared extinction can be determined. The mechanical luminosity for both lobes in a bipolar outflow is derived from formulae of the form Lmech = Mco $`v_{co}^3`$/(2 $`l`$) where Mco is derived from L(CO R(1)), assuming a CO abundance and temperature, $`v_{co}`$ represents the appropriately averaged speed of the mass that dominates the outflow energy (and must account for the source orientation when derived) and $`l`$ is the outflow size.

From the ratio of mass to Lco in the log $`v`$ bins of Figure 7, we estimate from the simulations that MCO for one side of a bipolar outflow $``$ 300 L(CO R(1)) in solar units. Furthermore, given the consistently low values of $`\gamma `$, we take $`v_{co}`$ = 20 km s-1. Finally, we fix L(H2 1–0 S(1)/L(H2) = 0.05. It must be remarked that a rather high number of assumptions precede a comparison of simulation and observation. Nevertheless, it is clear that the ratio L(H2 1–0 S(1))/L(CO R(1)) $``$ 3–8 for Run D after 1000 yr, while L(H2 1–0 S(1))/L(CO R(1)) $``$ 30–40 for Run I after 1000 yr. Given the above factors these convert to: L(H2)/Lmech $``$ 0.02–0.05 for Run D and $``$ 0.2–0.3 for Run I.

We confirm that strong H2 flows are thus associated with Run I. For example, we find L(H2)/Lmech $``$ 0.3 for Cepheus E, one of the dynamically youngest outflows (Smith et al., 2003). Moreover, the H2 emission is widespread and the CO position-velocity diagram displays structure cutting through all radial velocities at all positions, rather than the triangular Hubble-law features (Ladd & Hodapp, 1997), consistent with the Run I predictions. Finally, we note that the CO $`\gamma `$ values for Cepheus E are 0.99 and 1.8, for the two lobes, although these values are quite sensitive to the chosen local standard of rest (Smith et al., 1997). Hence, observations of Cepheus E are consistent with our results for a newly-forming outflow.

To date, CO emission has not been fully resolved in many outflows. Velocity channel maps are thus rare. A highly collimated CO outflow emanates from IRAS 21391+5802, with a V-shaped morphology on velocity channel images (Beltrán et al., 2002). The V-shape suggests that it is also an immature outflow in which the central source is increasing in power. We would thus expect to find a strong H2 outflow, as is, indeed, the case: the total shocked emission is estimated to be 4 L whereas the mechanical power estimates lie in the range 0.15–1.2 L (Nisini et al., 2001).

In contrast, there are many examples of outflows with very weak total H2 emission. These include the HH 366/B5 IRS1 outflow in the Perseus cloud complex, which possesses high collimation. Hubble-law behaviour, consistent with a single outburst in the past is also present (Yu et al., 1999). The H outflow in OMC-2 is another example, where the presently-observed H2 jet does not appear to be the driver for the CO outflow (Yu et al., 2000).

IRAS 18148-0440 in the L483 core is thought to be a protostar undergoing the transition from Class 0 to Class 1. Its CO outflow may possess a T-shaped morphology, although background cloud subtraction makes this uncertain (Tafalla et al., 2000). It belongs to a class of objects characterised by low velocity CO and no CO bullets, as also predicted here for a transitional outflow (Tafalla et al., 2000).

## 7 Summary

We have reported on two jet simulations that model drastic jet density evolution over tens or hundreds of years. The basic physical structure is similar to that previously found for non-evolving sheared molecular jets (Rosen & Smith, 2003), with an accelerated leading bow shock followed by a wider, slower shoulder-like region. More specifically, both of the simulations on display here resemble a light molecular jet simulation. The primary reason for this is that the shoulder region makes only slow progress relative to the leading edge of the bow shock.

We have explored predicted features relevant to high resolution near-infrared (H2 rovibrational at 2.12$`\mu `$m emission), far-infrared (CO R(5) rotational emission at 435$`\mu `$m) and millimetre (CO R(1) rotational emission at 1.3mm) wavelength observations. As expected, the H2 lines reveal patchy details of the structure and the CO lines show the general outline of material encompassed by the bow. CO line emission from mid-level rotational states show both the general outline and the internal shock structure. Owing to the small amount of higher temperature gas in the aborted jet simulation, the H2 emission is limited to the very tip of the accelerated bow shock. The predicted values of the 1–0 S(1) 2.12$`\mu `$m luminosity are comparable to those derived for the majority of outflows in the unbiased statistical study of Stanke et al. (2002).

The evolution in the integrated luminosity in these molecular lines is evident in both simulations. The changes in the integrated values of H2 occur at program times that are within 100 yr after the modification of the inflow jet density. Of particular note is the ratio of the integrated H2 to the integrated CO R(1) line in the aborted jet case, which decreases dramatically with time. Thus, in principle this ratio can be used to estimate the age of an outburst or aborted molecular flow from a protostellar core.

We have also computed velocity distributions in mass, CO luminosity, and H2 luminosity. Here the main results are:

1. The slopes of these distributions are small, closer to one than two in most cases. Steeper slopes ($`\gamma `$ = 2–3) are seen at high velocities in CO in the aborted jet simulation.

2. The slopes for the aborted jet case are larger than those for the increasing jet density case, although not by a large amount for the simulations listed in Table 1.

3. Similar to non-evolving jets (Rosen & Smith, 2003), $`\gamma _{\mathrm{CO}}`$ is usually greater than $`\gamma _{\mathrm{mass}}`$, but there are significant exceptions.

4. In the aborted jet simulation, there is no excess in the amount of mass or CO at the highest velocities; we have seen this excess that is associated with the internal shocks of a pulsed jet.

The position-velocity maps for the aborted jet case are relatively featureless. In CO emission, a triangular Hubble-Law outline is filled in by the emission but does not include the low velocity emission near the jet inflow boundary. The CO R(5) position-velocity map shows some fine structure, with brighter regions associated with the warmer parts of the flow: the leading part of the bow shock and some internal structure near the front of the shoulder region.

In the increasing jet density simulation, the position-velocity CO maps show the series of Hubble Law regions associated with successive pulses in the jet. Additionally, the CO R(5) map has more structure than the CO R(1) map. In the H2 position-velocity map, there are numerous features that are narrow along the jet axis but cover a large range of velocities (the vertical lines in Fig. 9).

The CO R(1) velocity channel maps for the aborted jet simulation have emission confined to a small range of velocities. The sequence of channel images shows a T-morphology at low velocities and the characteristics typical of a Hubble Law at high velocities. The channel maps for the increasing jet density simulation possess a V-morphology, with the inclusion of an additional protuberance from the new denser flow that carves a fresh path near the previous flow.

These features, combined with the levels of integrated emissions, suggest that jets in certain outflows (e.g. IRAS 21391+5802 and Cepheus E, see Sect. 6) are increasing in power, while other jets (e.g. that from IRAS 18148-0440 in L483) are on the decline.

Not many detailed velocity channel maps have as yet reached the literature. This will change shortly, and with the advent of new telescopes such as SIRTF, FIRST and ALMA, imaging and spectroscopy will reveal the properties on fine scales in many molecular lines of intermediate and low excitation, such as in the CO R(5) and R(1) lines. Here we have determined how these properties can help determine the jet evolutionary stage.

## Acknowledgments

The numerical calculations were run on the local SGI Origin 2000 computer (FORGE), acquired through the PPARC JREI initiative with SGI participation. AR is funded by PPARC.


format: text

text:
# 1 Introduction

## 1 Introduction

Light-front wavefunctions are the amplitudes which interpolate between hadrons and their quark and gluon degrees of freedom in QCD. For example, the eigensolution of a meson, projected on the eigenstates $`\{|n\}`$ of the free Hamiltonian $`H_{LC}^{QCD}(g=0)`$ at fixed light-front time $`\tau =t+z/c`$ with the same global quantum numbers, has the expansion:

$`|\mathrm{\Psi }_M;P^+,\stackrel{}{P}_{},\lambda `$ $`=`$ $`{\displaystyle \underset{n2,\lambda _i}{}}{\displaystyle \mathrm{\Pi }_{i=1}^n\frac{d^2k_idx_i}{\sqrt{x_i}16\pi ^3}}`$
$`\times 16\pi ^3\delta \left(1{\displaystyle \underset{j}{\overset{n}{}}}x_j\right)\delta ^{(2)}\left({\displaystyle \underset{\mathrm{}}{\overset{n}{}}}\stackrel{}{k}_{\mathrm{}}\right)`$
$`\times |n;x_iP^+,x_i\stackrel{}{P}_{}+\stackrel{}{k}_i,\lambda _i\psi _{n/M}(x_i,\stackrel{}{k}_i,\lambda _i).`$
The set of light-front Fock state wavefunctions $`\{\psi _{n/M}\}`$ represents the ensemble of quark and gluon states possible when the meson is intercepted at the light-front. The light-front momentum fractions $`x_i=k_i^+/P_\pi ^+=(k^0+k_i^z)/(P^0+P^z)`$ with $`_{i=1}^nx_i=1`$ and $`\stackrel{}{k}_i`$ with $`_{i=1}^n\stackrel{}{k}_i=\stackrel{}{0}_{}`$ represent the relative momentum coordinates of the QCD constituents; the scalar light-front wavefunctions $`\psi _{n/p}(x_i,\stackrel{}{k}_i,\lambda _i)`$ are independent of the proton’s momentum $`P^+=P^0+P^z`$, and $`P_{}`$. The physical transverse momenta are $`\stackrel{}{p}_i=x_i\stackrel{}{P}_{}+\stackrel{}{k}_i.`$ The $`\lambda _i`$ label the light-front spin $`S^z`$ projections of the quarks and gluons along the quantization $`z`$ direction. The spinors of the light-front formalism automatically incorporate the Melosh-Wigner rotation. Light-cone gauge $`A^+=0`$ is used to eliminate unphysical gauge degrees of freedom. The gluon polarization vectors $`ϵ^\mu (k,\lambda =\pm 1)`$ are specified in light-cone gauge by the conditions $`kϵ=0,\eta ϵ=ϵ^+=0.`$ The quark and gluon degrees of freedom are all physical; there are effectively no ghost or negative metric states.

An important feature of the light-front formalism is that the projection $`J_z`$ of the total angular momentum is kinematical and conserved. Each light-front Fock state component satisfies the angular momentum sum rule: $`J^z=_{i=1}^nS_i^z+_{j=1}^{n1}l_j^z.`$ The summation over $`S_i^z`$ represents the contribution of the intrinsic spins of the $`n`$ Fock state constituents. The summation over orbital angular momenta

$$l_j^z=\mathrm{i}\left(k_j^x\frac{}{k_j^y}k_j^y\frac{}{k_j^x}\right)$$
(2)
derives from the $`n1`$ relative momenta. This excludes the contribution to the orbital angular momentum due to the motion of the center of mass, which is not an intrinsic property of the hadron. The light-front eigensolution corresponds to a spin $`J`$ particle in the hadron rest frame $`P^+=P^{}=M,\stackrel{}{P}_{}=\stackrel{}{0}`$, not the constituent rest frame $`_i\stackrel{}{k}_i=\stackrel{}{0}`$ since $`_ik_i^zP^z.`$ The numerator structure of the light-front wavefunctions is in large part determined by the angular momentum constraints. Thus wavefunctions generated by perturbation theory can provide a template for the numerator structure of nonperturbative light-front wavefunctions.

Hadronic amplitudes can be computed by inserting a sum over a complete sets of free Fock states for each external hadron, thus representing the dynamics of each hadron as a convolution of its light-front wavefunctions with the corresponding $`n`$-particle irreducible quark-gluon matrix elements, summed over $`n.`$ For example, in the case of spacelike form factors, the matrix elements of local currents are given by a simple overlap of light-front wavefunctions. If one chooses the frame with $`q^+=0`$, then matrix elements of currents such as $`j^+`$ in electroweak theory and have only diagonal matrix elements $`n^{}=n.`$ Thus once one has solved for the light-front wavefunctions, one can compute hadron matrix elements of currents between hadronic states of arbitrary momentum. Remarkably, quantum fluctuations of the vacuum are absent if one uses light-front time to quantize the system, so that matrix elements such as the electromagnetic form factors only depend on the currents of the constituents described by the light-cone wavefunctions. As I discuss below, the degrees of freedom associated with vacuum phenomena such as spontaneous symmetry breaking in the Higgs model have their counterpart in light-front $`k^+=0`$ zero modes of the fields.

Matrix elements of spacelike currents such as spacelike electromagnetic form factors thus have an exact representation in terms of simple overlaps of the light-front wavefunctions in momentum space with the same $`x_i`$ and unchanged parton number $`n`$. The Pauli form factor and anomalous moment are spin-flip matrix elements of $`j^+`$ and thus connect states with $`\mathrm{\Delta }L_z=1.`$ Thus, these quantities are nonzero only if there is nonzero orbital angular momentum of the quarks in the proton. The Dirac form factor is diagonal in $`L_z`$ and is typically dominated at high $`Q^2`$ by highest states with the highest orbital angular momentum. In the case of nuclear form factors, Fock states with “hidden color” play an important role, particularly at large momentum transfer. The formulae for electroweak current matrix elements of $`j^+`$ can be easily extended to the $`T^{++}`$ coupling of gravitons. In, fact, one can show that the anomalous gravito-magnetic moment $`B(0)`$, analogous to $`F_2(0)`$ in electromagnetic current interactions, vanishes identically for any system, composite or elementary. This important feature, which follows in general from the equivalence principle , is obeyed explicitly in the light-front formalism.

The light-front Fock representation is especially advantageous in the study of exclusive $`B`$ decays. For example, we can write down an exact frame-independent representation of decay matrix elements such as $`BD\mathrm{}\overline{\nu }`$ from the overlap of $`n^{}=n`$ parton conserving wavefunctions and the overlap of $`n^{}=n2`$ from the annihilation of a quark-antiquark pair in the initial wavefunction. The off-diagonal $`n+1n1`$ contributions give a new perspective for the physics of $`B`$-decays. A semileptonic decay involves not only matrix elements where a quark changes flavor, but also a contribution where the leptonic pair is created from the annihilation of a $`q\overline{q^{}}`$ pair within the Fock states of the initial $`B`$ wavefunction. The semileptonic decay thus can occur from the annihilation of a nonvalence quark-antiquark pair in the initial hadron. Intrinsic charm $`|b\overline{u}c\overline{c}`$ states of the $`B`$ meson, although small in probability, can play an important role in its weak decays because they facilitate CKM-favored weak decays. The “handbag” contribution to the leading-twist off-forward parton distributions measured in deeply virtual Compton scattering has a similar light-front wavefunction representation as overlap integrals of light-front wavefunctions.

In the case of hadronic amplitudes involving a hard momentum transfer $`Q`$, it is often possible to expand the quark-gluon scattering amplitude as a function of $`k_{}^2/Q^2`$. The leading-twist contribution then can be computed from a hard-scattering amplitude $`T_H`$ where the external quarks and gluons emanating from each hadron can be taken as collinear. The convolution with the light-front wavefunction and integration $`\mathrm{\Pi }_id^2k_i`$ over the relative transverse momentum projects out only the $`L_z=0`$ component of the light-front wavefunctions. This leads to hadron spin selection rules such as hadron helicity conservation. Furthermore, only the minimum number of quark and gluon quanta contribute at leading order in $`1/Q^2.`$ The nominal scaling of hard hadron scattering amplitudes at leading twist then obeys dimensional counting rules. Recently these rules have been derived to all orders in the gauge coupling in conformal QCD and large $`N_C`$ using gauge/string duality. There is also evidence from hadronic $`\tau `$ decays that the QCD coupling approaches an infrared fixed-point at low scales. This may explain the empirical success of conformal approximations to QCD. The distribution amplitudes $`\varphi (x_i,Q)`$ which appear in factorization formulae for hard exclusive processes are the valence LF Fock wavefunctions integrated over the relative transverse momenta up to the resolution scale $`Q`$. These quantities specify how a hadron shares its longitudinal momentum among its valence quarks; they control virtually all exclusive processes involving a hard scale $`Q`$, including form factors, Compton scattering, semi-exclusive processes, and photoproduction at large momentum transfer, as well as the decay of a heavy hadron into specific final states.

The quark and gluon probability distributions $`q_i(x,Q)`$ and $`g(x,Q)`$ of a hadron can be computed from the absolute squares of the light-front wavefunctions, integrated over the transverse momentum. All helicity distributions are thus encoded in terms of the light-front wavefunctions. The DGLAP evolution of the structure functions can be derived from the high $`k_{}`$ properties of the light-front wavefunctions. Thus given the light-front wavefunctions, one can compute all of the leading twist helicity and transversity distributions measured in polarized deep inelastic lepton scattering. Similarly, the transversity distributions and off-diagonal helicity convolutions are defined as a density matrix of the light-front wavefunctions.

However, it is not true that the leading-twist structure functions $`F_i(x,Q^2)`$ measured in deep inelastic lepton scattering are identical to the quark and gluon distributions. It is usually assumed, following the parton model, that the $`F_2`$ structure function measured in neutral current deep inelastic lepton scattering is at leading order in $`1/Q^2`$ simply $`F_2(x,Q^2)=_qe_q^2xq(x,Q^2)`$, where $`x=x_{bj}=Q^2/2pq`$ and $`q(x,Q)`$ can be computed from the absolute square of the proton’s light-front wavefunction. Recent work by Hoyer, Marchal, Peigne, Sannino, and myself shows that this standard identification is wrong. In fact, one cannot neglect the Wilson line integral between currents in the current correlator even in light-cone gauge. In the case of light-cone gauge, the Wilson line involves the transverse gluon field $`A_{}`$ not $`A^+.`$ Gluon exchange between the fast, outgoing partons and the target spectators affects the leading-twist structure functions in a profound way. The final-state interactions lead to the Bjorken-scaling diffractive component $`\gamma ^{}ppX`$ of deep inelastic scattering. The diffractive scattering of the fast outgoing quarks on spectators in the target in turn causes shadowing in the DIS cross section. Thus the depletion of the nuclear structure functions is not intrinsic to the wave function of the nucleus, but is a coherent effect arising from the destructive interference of diffractive channels induced by final-state interactions. Similarly, the effective Pomeron distribution of a hadron is not derived from its light-front wavefunction and thus is not a universal property. Many properties involving parton transverse momentum are also affected by the Wilson line.

Measurements from the HERMES and SMC collaborations show a remarkably large single-spin asymmetry in semi-inclusive pion leptoproduction $`\gamma ^{}(q)p\pi X`$ when the proton is polarized normal to the photon-to-pion production plane. Hwang, Schmidt, and I  have shown that final-state interactions from gluon exchange between the outgoing quark and the target spectator system lead to single-spin asymmetries in deep inelastic lepton-proton scattering at leading twist in perturbative QCD; i.e., the rescattering corrections are not power-law suppressed at large photon virtuality $`Q^2`$ at fixed $`x_{bj}`$. The existence of such single-spin asymmetries requires a phase difference between two amplitudes coupling the proton target with $`J_p^z=\pm \frac{1}{2}`$ to the same final-state, the same amplitudes which are necessary to produce a nonzero proton anomalous magnetic moment. The single-spin asymmetry which arises from such final-state interactions does not factorize into a product of distribution function and fragmentation function, and it is not related to the transversity distribution $`\delta q(x,Q)`$ which correlates transversely polarized quarks with the spin of the transversely polarized target nucleon. In general all measures of quark and gluon transverse momentum require consideration of final-state interactions as incorporated in the Wilson line.

These effects highlight the unexpected importance of final- and initial-state interactions in QCD observables—they lead to leading-twist single-spin asymmetries, diffraction, and nuclear shadowing, phenomena not included in the light-front wavefunctions of the target. Alternatively, as discussed by Belitsky, Ji, and Yuan, one can augment the light-front wavefunctions by including the phases induced by initial and final state interactions. Such wavefunctions correspond to solving the light-front bound state equation in an external field.

## 2 The Light-Front Quantization of QCD

In Dirac’s “Front Form”, the generator of light-front time translations is $`P^{}=i\frac{}{\tau }.`$ Boundary conditions are set on the transverse plane labelled by $`x_{}`$ and $`x^{}=zct`$. Given the Lagrangian of a quantum field theory, $`P^{}`$ can be constructed as an operator on the Fock basis, the eigenstates of the free theory. Since each particle in the Fock basis is on its mass shell, $`k^{}k^0k^3=\frac{k_{}^2+m^2}{k^+},`$ and its energy $`k^0=\frac{1}{2}(k^++k^{})`$ is positive, only particles with positive momenta $`k^+k^0+k^30`$ can occur in the Fock basis. Since the total plus momentum $`P^+=_nk_n^+`$ is conserved, the light-cone vacuum cannot have any particle content.

The Heisenberg equation on the light-front is

$$H_{LC}|\mathrm{\Psi }=M^2|\mathrm{\Psi }.$$
(3)
The operator $`H_{LC}=P^+P^{}P_{}^2,`$ the “light-cone Hamiltonian”, is frame-independent. This can in principle be solved by diagonalizing the matrix $`n|H_{LC}|m`$ on the free Fock basis:

$$\underset{m}{}n|H_{LC}|mm|\psi =M^2n|\mathrm{\Psi }.$$
(4)
The eigenvalues $`\{M^2\}`$ of $`H_{LC}=H_{LC}^0+V_{LC}`$ give the squared invariant masses of the bound and continuum spectrum of the theory. The light-front Fock space is the eigenstates of the free light-front Hamiltonian; i.e., it is a Hilbert space of non-interacting quarks and gluons, each of which satisfy $`k^2=m^2`$ and $`k^{}=\frac{m^2+k_{}^2}{k^+}0.`$ The projections $`\{n|\mathrm{\Psi }\}`$ of the eigensolution on the $`n`$-particle Fock states provide the light-front wavefunctions. Thus solving a quantum field theory is equivalent to solving a coupled many-body quantum mechanical problem:

$$\left[M^2\underset{i=1}{\overset{n}{}}\frac{m^2+k_{}^2}{x_i}\right]\psi _n=\underset{n^{}}{}n|V_{LC}|n^{}\psi _n^{}$$
(5)
where the convolution and sum is understood over the Fock number, transverse momenta, plus momenta, and helicity of the intermediate states. Light-front wavefunctions are also related to momentum-space Bethe-Salpeter wavefunctions by integrating over the relative momenta $`k^{}=k^0k^z`$ since this projects out the dynamics at $`x^+=0.`$

A review of the development of light-front quantization of QCD and other quantum field theories is given in the references. The light-front quantization of gauge theory can be most conveniently carried out in the light-cone gauge $`A^+=A^0+A^z=0`$. In this gauge the $`A^{}`$ field becomes a dependent degree of freedom, and it can be eliminated from the Hamiltonian in favor of a set of specific instantaneous light-front time interactions. In fact in $`QCD(1+1)`$ theory, this instantaneous interaction provides the confining linear $`x^{}`$ interaction between quarks. In $`3+1`$ dimensions, the transverse field $`A^{}`$ propagates massless spin-one gluon quanta with polarization vectors which satisfy both the gauge condition $`ϵ_\lambda ^+=0`$ and the Lorentz condition $`kϵ=0`$.

Prem Srivastava and I have presented a new systematic study of light-front-quantized gauge theory in light-cone gauge using a Dyson-Wick S-matrix expansion based on light-front-time-ordered products. The Dirac bracket method is used to identify the independent field degrees of freedom. In our analysis one imposes the light-cone gauge condition as a linear constraint using a Lagrange multiplier, rather than a quadratic form. We then find that the LF-quantized free gauge theory simultaneously satisfies the covariant gauge condition $`A=0`$ as an operator condition as well as the LC gauge condition. The gluon propagator has the form

$$0|T(A_{}^{a}{}_{\mu }{}^{}(x)A_{}^{b}{}_{\nu }{}^{}(0))|0=\frac{i\delta ^{ab}}{(2\pi )^4}d^4ke^{ikx}\frac{D_{\mu \nu }(k)}{k^2+iϵ}$$
(6)
where we have defined

$$D_{\mu \nu }(k)=D_{\nu \mu }(k)=g_{\mu \nu }+\frac{n_\mu k_\nu +n_\nu k_\mu }{(nk)}\frac{k^2}{(nk)^2}n_\mu n_\nu .$$
(7)
Here $`n_\mu `$ is a null four-vector, gauge direction, whose components are chosen to be $`n_\mu =\delta _{\mu }^{}{}_{}{}^{+}`$, $`n^\mu =\delta _{}^{\mu }{}_{}{}^{}`$. Note also

$`D_{\mu \lambda }(k)D_{}^{\lambda }{}_{\nu }{}^{}(k)=D_\mu (k)D_{}^{}{}_{\nu }{}^{}(k)`$ $`=`$ $`D_{\mu \nu }(k),`$ (8)
$`k^\mu D_{\mu \nu }(k)=0,n^\mu D_{\mu \nu }(k)`$ $``$ $`D_\nu (k)=0,`$
$`D_{\lambda \mu }(q)D^{\mu \nu }(k)D_{\nu \rho }(q^{})`$ $`=`$ $`D_{\lambda \mu }(q)D^{\mu \rho }(q^{}).`$
The gauge field propagator $`iD_{\mu \nu }(k)/(k^2+iϵ)`$ is transverse not only to the gauge direction $`n_\mu `$ but also to $`k_\mu `$, i.e., it is doubly-transverse. Thus $`D`$ represents the polarization sum over physical propagating modes. The last term proportional to $`n_\mu n_\nu `$ in the gauge propagator does not appear in the usual formulations of light-cone gauge. However, in tree graph calculations it cancels against instantaneous gluon exchange contributions.

The remarkable properties of (the projector) $`D_{\nu \mu }`$ greatly simplifies the computations of loop amplitudes. For example, the coupling of gluons to propagators carrying high momenta is automatic. In the case of tree graphs, the term proportional to $`n_\mu n_\nu `$ cancels against the instantaneous gluon exchange term. However, in the case of loop diagrams, the separation needs to be maintained so that one can identify the correct one-particle-irreducible contributions. The absence of collinear divergences in irreducible diagrams in the light-cone gauge greatly simplifies the leading-twist factorization of soft and hard gluonic corrections in high momentum transfer inclusive and exclusive reactions since the numerators associated with the gluon coupling only have transverse components.

The interaction Hamiltonian of QCD in light-cone gauge can be derived by systematically applying the Dirac bracket method to identify the independent fields. It contains the usual Dirac interactions between the quarks and gluons, the three-point and four-point gluon non-Abelian interactions plus instantaneous light-front-time gluon exchange and quark exchange contributions

$`_{int}`$ $`=`$ $`g\overline{\psi }^i\gamma ^\mu A_{\mu }^{}{}_{}{}^{ij}\psi ^j`$ (9)
$`+{\displaystyle \frac{g}{2}}f^{abc}(_\mu A_{}^{a}{}_{\nu }{}^{}_\nu A_{}^{a}{}_{\mu }{}^{})A^{b\mu }A^{c\nu }`$
$`+{\displaystyle \frac{g^2}{4}}f^{abc}f^{ade}A_{b\mu }A^{d\mu }A_{c\nu }A^{e\nu }`$
$`{\displaystyle \frac{g^2}{2}}\overline{\psi }^i\gamma ^+(\gamma ^{^{}}A_{^{}})^{ij}{\displaystyle \frac{1}{i_{}}}(\gamma ^{}A_{})^{jk}\psi ^k`$
$`{\displaystyle \frac{g^2}{2}}j_{}^{+}{}_{a}{}^{}{\displaystyle \frac{1}{(_{})^2}}j_{}^{+}{}_{a}{}^{}`$
where

$$j_{}^{+}{}_{a}{}^{}=\overline{\psi }^i\gamma ^+(t_a)^{ij}\psi ^j+f_{abc}(_{}A_{b\mu })A^{c\mu }.$$
(10)
The renormalization constants in the non-Abelian theory have been shown  to satisfy the identity $`Z_1=Z_3`$ at one-loop order, as expected in a theory with only physical gauge degrees of freedom. The renormalization factors in the light-cone gauge are independent of the reference direction $`n^\mu `$. The QCD $`\beta `$ function computed in the noncovariant LC gauge agrees with the conventional theory result. Dimensional regularization and the Mandelstam-Leibbrandt prescription for LC gauge were used to define the Feynman loop integrations. There are no Faddeev-Popov or Gupta-Bleuler ghost terms.

The running coupling constant and the QCD $`\beta `$ function have also been computed at one loop in the doubly-transverse light-cone gauge. It is also possible to effectively quantize QCD using light-front methods in covariant Feynman gauge. It is well-known that the light-cone gauge itself is not completely defined until one specifies a prescription for the poles of the gauge propagator at $`nk=0.`$ The Mandelstam-Liebbrandt prescription has the advantage of preserving causality and analyticity, as well as leading to proofs of the renormalizability and unitarity of Yang-Mills theories. The ghosts which appear in association with the Mandelstam-Liebbrandt prescription from the single poles have vanishing residue in absorptive parts, and thus do not disturb the unitarity of the theory.

A remarkable advantage of light-front quantization is that the vacuum state $`|\mathrm{\hspace{0.17em}0}`$ of the full QCD Hamiltonian evidently coincides with the free vacuum. The light-front vacuum is effectively trivial if the interaction Hamiltonian applied to the perturbative vacuum is zero. Note that all particles in the Hilbert space have positive energy $`k^0=\frac{1}{2}(k^++k^{})`$, and thus positive light-front $`k^\pm `$. Since the plus momenta $`k_i^+`$ is conserved by the interactions, the perturbative vacuum can only couple to states with particles in which all $`k_i^+`$ = 0; i.e., so called zero-mode states. Bassetto and collaborators have shown that the computation of the spectrum of $`QCD(1+1)`$ in equal time quantization requires constructing the full spectrum of non perturbative contributions (instantons). In contrast, in the light-front quantization of gauge theory, where the $`k^+=0`$ singularity of the instantaneous interaction is defined by a simple infrared regularization, one obtains the correct spectrum of $`QCD(1+1)`$ without any need for vacuum-related contributions. Zero modes of auxiliary fields are necessary to distinguish the theta-vacua of massless QED(1+1) , or to represent a theory in the presence of static external boundary conditions or other constraints. Zero-modes provide the light-front representation of spontaneous symmetry breaking in scalar theories.

## 3 Light-Front Quantization of the Standard Model

Prem Srivastava and I have also shown how light-front quantization can be applied to the Glashow, Weinberg and Salam (GWS) model of electroweak interactions based on the nonabelian gauge group $`SU(2)_W\times U(1)_Y`$. This theory contains a nonabelian Higgs sector which triggers spontaneous symmetry breaking (SSB). A convenient way of implementing SSB and the (tree level) Higgs mechanism in the front form theory was developed earlier by Srivastava. One separates the quantum fluctuation fields from the corresponding dynamical bosonic condensate  (or zero-longitudinal-momentum-mode) variables, before applying the Dirac procedure in order to construct the Hamiltonian formulation. The canonical quantization of LC gauge GWS electroweak theory in the front form can be derived by using the Dirac procedure to construct a self-consistent LF Hamiltonian theory. This leads to an attractive new formulation of the Standard Model of the strong and electroweak interactions which does not break the physical vacuum and has well-controlled ultraviolet behavior. The only ghosts which appear in the formalism are the $`nk=0`$ modes of the gauge field associated with regulating the light-cone gauge prescription. The massive gauge field propagator has good asymptotic behavior in accordance with a renormalizable theory, and the massive would-be Goldstone fields can be taken as physical degrees of freedom.

For example, consider the Abelian Higgs model. The interaction Lagrangian is

$$=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+|D_\mu \varphi |^2V(\varphi ^{}\varphi )$$
(11)
where

$$D_\mu =_\mu +ieA_\mu ,$$
(12)
and

$$V(\varphi )=\mu ^2\varphi ^{}\varphi +\lambda (\varphi ^{}\varphi )^2,$$
(13)
with $`\mu ^2<0,\lambda >0.`$ The complex scalar field $`\varphi `$ is decomposed as

$$\varphi (x)=\frac{1}{\sqrt{2}}v+\phi (x)=\frac{1}{\sqrt{2}}[v+h(x)+i\eta (x)]$$
(14)
where $`v`$ is the $`k^+=0`$ zero mode determined by the minimum of the potential: $`v^2=\frac{\mu ^2}{\lambda }`$, $`h(x)`$ is the dynamical Higgs field, and $`\eta (x)`$ is the Nambu-Goldstone field. The quantization procedure determines $`A=MR`$, the ’t Hooft condition. One can now eliminate the zero mode component of the Higgs field $`v`$ which gives masses for the fundamental quantized fields. The $`A_{}`$ field then has mass $`M=ev`$ and the Higgs field acquires mass $`m_h^2=2\lambda v^2=2\mu ^2.`$

A new aspect of LF quantization, is that the third polarization of the quantized massive vector field $`A^\mu `$ with four momentum $`k^\mu `$ has the form $`E_\mu ^{(3)}=n_\mu M/nk`$. Since $`n^2=0`$, this non-transverse polarization vector has zero norm. However, when one includes the constrained interactions of the Goldstone particle, the effective longitudinal polarization vector of a produced vector particle is $`E_{\mathrm{eff}\mu }^{(3)}=E_\mu ^{(3)}k_\mu kE^{(3)}/k^2`$ which is identical to the usual polarization vector of a massive vector with norm $`E_{\mathrm{eff}}^{(3)}E_{\mathrm{eff}}^{(3)}=1`$. Thus, unlike the conventional quantization of the Standard Model, the Goldstone particle only provides part of the physical longitudinal mode of the electroweak particles.

In the LC gauge LF framework, the free massive gauge fields in the electroweak theory satisfy simultaneously the ’t Hooft conditions as an operator equation. The sum over the three physical polarizations is given by $`K_{\mu \nu }`$

$`K_{\mu \nu }(k)`$ $`=`$ $`{\displaystyle \underset{(\alpha )}{}}E_\mu ^{(\alpha )}E_\nu ^{(\alpha )}=D_{\mu \nu }(k)+{\displaystyle \frac{M^2}{(k^+)^2}}n_\mu n_\nu `$
$`=`$ $`g_{\mu \nu }+{\displaystyle \frac{n_\mu k_\nu +n_\nu k_\mu }{(nk)}}{\displaystyle \frac{(k^2M^2)}{(nk)^2}}n_\mu n_\nu `$
which satisfies: $`k^\mu K_{\mu \nu }(k)=(M^2/k^+)n_\nu `$ and $`k^\mu k^\nu K_{\mu \nu }(k)=M^2`$. The free propagator of the massive gauge field $`A_\mu `$ is

$`0|T\left(A_\mu (x)A_\nu (y)\right)|0=`$ (16)
$`{\displaystyle \frac{i}{(2\pi )^4}}{\displaystyle d^4k\frac{K_{\mu \nu }(k)}{(k^2M^2+iϵ)}e^{ik(xy)}}.`$
It does not have the bad high energy behavior found in the (Proca) propagator in the unitary gauge formulation, where the would-be Nambu-Goldstone boson is gauged away.

In the limit of vanishing mass of the vector boson, the gauge field propagator goes over to the doubly transverse gauge, ($`n^\mu D_{\mu \nu }(k)=k^\mu D_{\mu \nu }(k)=0`$), the propagator found  in QCD. The numerator of the gauge propagator $`K_{\mu \nu }(k)`$ also has important simplifying properties, similar to the ones associated with the projector $`D_{\mu \nu }(k)`$. The transverse polarization vectors for massive or massless vector boson may be taken to be $`E_{()}^\mu (k)D_{}^\mu (k),`$ whereas the non-transverse third one in the massive case is found to be parallel to the LC gauge direction $`E_\mu ^{(3)}(k)=(M/k^+)n_\mu `$. Its projection along the direction transverse to $`k_\mu `$ shares the spacelike vector property carried by $`E_{()}^\mu (k)`$. The Goldstone boson or electroweak equivalence theorem becomes transparent in the LF formulation.

The interaction Hamiltonian of the Standard Model can be written in a compact form by retaining the dependent components $`A^{}`$ and $`\psi _{}`$ in the formulation. Its form closely resembles the interaction Hamiltonian of covariant theory, except for the presence of additional instantaneous four-point interactions. The resulting Dyson-Wick perturbation theory expansion based on equal-LF-time ordering has also been constructed, allowing one to perform higher-order computations in a straightforward fashion. The singularities in the noncovariant pieces of the field propagators may be defined using the causal ML prescription for $`1/k^+`$ when we employ dimensional regularization, as was shown in our earlier work on QCD. The power-counting rules in LC gauge then become similar to those found in covariant gauge theory.

Spontaneous symmetry breaking is thus implemented in a novel way when one quantizes the Standard Model at fixed light-front time $`\tau =x^+.`$ In the general case, the Higgs field $`\varphi _i(x)`$ can be separated into two components:

$$\varphi _i(\tau ,x^{},\stackrel{}{x}_{})=\omega _i(\tau ,\stackrel{}{x}_{})+\phi (\tau ,x^{},\stackrel{}{x}_{}),$$
(17)
where $`\omega _i`$ is a classical $`k^+=0`$ zero-mode field and $`\phi `$ is the dynamical quantized field. Here $`i`$ is the weak-isospin index. The zero-mode component is determined by solving the Euler-Lagrange tree-level condition:

$$V_i^{}(\omega )_{}_{}\omega _i=0.$$
(18)
A nonzero value for $`\omega _i`$ corresponds to spontaneous symmetry breaking. The nonzero $`\omega _i`$ couples to the gauge boson and Fermi fields through the Yukawa interactions of the Standard Model. It can then be eliminated from the theory in favor of mass terms for the fundamental matter fields in the effective theory. The resulting masses are identical to those of the usual Higgs implementation of spontaneous symmetry breaking in the Standard Model.

The generators of isospin rotations are defined from the dynamical Higgs fields:

$$G_a=i𝑑x^{}𝑑x^{}(_{}\phi )_i(t_a)_{ij}\phi _j.$$
(19)
Note that the weak-isospin charges and the currents corresponding to $`G_a`$ are not conserved if the zero mode $`\omega _i`$ is nonzero since the cross terms in $`\phi ,`$ and $`\omega `$ are missing. Thus $`[H_{LF},G_a]0.`$ Nevertheless, the charges annihilate the vacuum: $`G_a|\mathrm{\hspace{0.17em}0}_{LF}=0,`$ since the dynamical fields $`\phi _i`$ have no support on the LF vacuum, and all quanta have positive $`k^+.`$ Thus the LF vacuum remains equal to the perturbative vacuum; it is unaffected by the occurrence of spontaneous symmetry breaking.

In effect one can interpret the $`k^+=0`$ zero mode field $`\omega _i`$ as an $`x^{}`$-independent external field, analogous to an applied constant electric or magnetic field in atomic physics. In this interpretation, the zero mode is a remnant of a Higgs field which persists from early cosmology; the LF vacuum however remains unchanged and unbroken.

## 4 Non-Perturbative Methods

As noted in section 2., solving a quantum field theory at fixed light-front time $`\tau `$ can be formulated as a relativistic extension of Heisenberg’s matrix mechanics. If one imposes periodic boundary conditions in $`x^{}=t+z/c`$, then the $`+`$ momenta become discrete: $`k_i^+=\frac{2\pi }{L}n_i,P^+=\frac{2\pi }{L}K`$, where $`_in_i=K`$. For a given “harmonic resolution” $`K`$, there are only a finite number of ways a set of positive integers $`n_i`$ can sum to a positive integer $`K`$. Thus at a given $`K`$, the dimension of the resulting light-front Fock state representation of the bound state is rendered finite without violating Lorentz invariance. The eigensolutions of a quantum field theory, both the bound states and continuum solutions, can then be found by numerically diagonalizing a frame-independent light-front Hamiltonian $`H_{LC}`$ on a finite and discrete momentum-space Fock basis. The continuum limit is reached for $`K\mathrm{}.`$ This formulation of the non-perturbative light-front quantization problem is called “discretized light-cone quantization” (DLCQ). The method preserves the frame-independence of the Front form.

The DLCQ method has been used extensively for solving one-space and one-time theories, including applications to supersymmetric quantum field theories and specific tests of the Maldacena conjecture. There has been progress in systematically developing the computation and renormalization methods needed to make DLCQ viable for QCD in physical spacetime. For example, John Hiller, Gary McCartor, and I  have shown how DLCQ can be used to solve 3+1 theories despite the large numbers of degrees of freedom needed to enumerate the Fock basis. A key feature of our work is the introduction of Pauli Villars fields to regulate the UV divergences and perform renormalization while preserving the frame-independence of the theory. A recent application of DLCQ to a 3+1 quantum field theory with Yukawa interactions is given in the references. One can also define a truncated theory by eliminating the higher Fock states in favor of an effective potential. Spontaneous symmetry breaking and other nonperturbative effects associated with the instant-time vacuum are hidden in dynamical or constrained zero modes on the light-front. An introduction is given by McCartor and Yamawaki

The pion distribution amplitude has been computed using a combination of the discretized DLCQ method for the $`x^{}`$ and $`x^+`$ light-front coordinates with a spatial lattice  in the transverse directions. A finite lattice spacing $`a`$ can be used by choosing the parameters of the effective theory in a region of renormalization group stability to respect the required gauge, Poincaré, chiral, and continuum symmetries. Dyson-Schwinger models  can also be used to predict light-front wavefunctions and hadron distribution amplitudes by integrating over the relative $`k^{}`$ momentum of the Bethe-Salpeter wavefunctions. Explicit nonperturbative light-front wavefunctions have been found in this way for the Wick-Cutkosky model, including spin-two states. One can also implement variational methods, using the structure of perturbative solutions as a template for the numerator of the light-front wavefunctions.

## 5 A Light-Front Event Amplitude Generator

The light-front formalism can be used as an “event amplitude generator” for high energy physics reactions where each particle’s final state is completely labelled in momentum, helicity, and phase. The application of the light-front time evolution operator $`P^{}`$ to an initial state systematically generates the tree and virtual loop graphs of the $`T`$-matrix in light-front time-ordered perturbation theory in light-cone gauge. Given the interactions of the light-front interaction Hamiltonian, any amplitude in QCD and the electroweak theory can be computed. For example, this method can be used to automatically compute the hard-scattering amplitudes $`T_H`$ for the deuteron form factor or $`pp`$ elastic scattering.

At higher orders, loop integrals only involve integrations over the momenta of physical quanta and physical phase space $`d^2k_idk_i^+`$. Renormalized amplitudes can be explicitly constructed by subtracting from the divergent loops amplitudes with nearly identical integrands corresponding to the contribution of the relevant mass and coupling counter terms – the “alternating denominator method”. The natural renormalization scheme to use for defining the coupling in the event amplitude generator is a physical effective charge such as the pinch scheme. The argument of the coupling is then unambiguous. The DLCQ boundary conditions can be used to discretize the phase space and limit the number of contributing intermediate states without violating Lorentz invariance. Since one avoids dimensional regularization and nonphysical ghost degrees of freedom, this method of generating events at the amplitude level could provide a simple but powerful tool for simulating events both in QCD and the Standard Model. Alternatively, one can construct the $`T`$matrix for scattering in QCD using light-front quantization and the event amplitude generator; one can then probe its spectrum by finding zeros of the resolvant. It would be particularly interesting to apply this method to finding the gluonium spectrum of QCD.

## Acknowledgments

Work supported by the Department of Energy under contract number DE-AC03-76SF00515. This talk is in large part based on collaborations with the late Professor Prem Srivastava. I am grateful to Professor Koichi Yamawaki and the other organizers of this meeting for their outstanding hospitality in Nagoya.