text:
# Molecular Hydrogen Outflows in W51

## 1 Introduction

The high-mass star-forming region W51 was discovered as a region of radio continuum emission by Westerhout (1958) and as a massive molecular cloud by Wilson, Jefferts, & Penzias (1971). It is one of the most luminous star-forming regions in our Galaxy and may be the closest analog to the extremely luminous star forming regions found in many other galaxies (e.g. 30 Doradus in the LMC). W51 is far more luminous than the well-studied relatively nearby Orion star forming region. The mass of the W51 complex has been estimated to be about 106 $`M_{}`$ (Carpenter and Sanders, 1998), based both on virial mass estimates and 12CO intensities. Mass estimates of molecular clouds in a complex environment obviously depend on where the boundaries are set, and this estimate refers to a set of emission features of approximately spherical shape. W51 is in the top 1% of all galactic molecular clouds by size, and the top 5 - 10% by mass.

W51 lies in the plane of the Galaxy at the substantial distance of 7.0 $`\pm `$ 1.5 Kpc measured by Genzel et al. (1981) through proper motion studies of the W51 Main H2O maser. W51 is therefore heavily obscured by interstellar extinction (AV = 24 $`\pm `$ 3 mag, as measured by Goldader & Wynn-Williams (1994) on the basis of the color of the extended emission). W51 is about 14 times farther away than Orion and the distance modulus is almost 6 magnitudes larger than for Orion, in addition to an extra 2 magnitudes of extinction in the K band, so that studies of W51 are necessarily limited to very bright objects and rather large extended features.

Numerous infrared sources were discovered in the radio continuum region W51. One of these, W51 IRS 2 (Wynn-Williams, Becklin, & Neugebauer, 1974), corresponding to the radio-continuum region W51e (Martin, 1972) is now known to be a young embedded star cluster larger than the Orion-Trapezium cluster, as first suggested by Genzel et al. (1982). A detailed study of massive stars in W51 was done by Okumura et al. (2000). They identified 4 subgroups within W51, roughly aligned parallel to the Galactic plane, with average ages in those subgroups ranging from less than 1 Myr to 2.3 Myr. Past studies of W51 have largely focussed on its global properties, e. g., the total luminosity and mass, and have traced the effects of star formation integrated over the lifetime of O-stars, by studying radio continuum or Br$`\gamma `$ emission .

This paper, by contrast, seeks to obtain a snapshot of star formation activity in the most recent past. Stars in their main accretion phase still have a substantial fraction of their final (ZAMS) mass residing in a circumstellar disk. They are therefore characterized by spectral energy distributions (SEDs) peaking at sub-mm wavelengths (André, Ward-Thompson, and Barsony, 1993) and are not directly detectable at near-infrared wavelengths. Fortunately, these class 0 sources usually also have strong outflows that can be traced by the shock-excited emission of H2 resulting from the interaction of the outflow with the ambient molecular cloud. Examples of successful searches for outflows using narrow-band imaging in the H2 1–0 S(1) line are, among many others, Hodapp & Ladd (1995), one of the earlier S(1) imaging surveys, and Stanke, McCaughrean, & Zinnecker (1998) who conducted a survey of large areas in the Orion GMC. Very strong outflow activity and shock-excited molecular hydrogen emission persists only through the class 0 and early class I phases, i.e. for less than 105y. In contrast to the study of HII regions, H2 outflows therefore trace contemporary lower mass star formation.

We will discuss the technical aspects of the observations in chapter 2. The results will be discussed in chapter 3, separated into discussions of the overall distribution of outflow sources in chapter 3.1 and of the morphology of individual outflows in chapter 3.2. Velocity resolved spectroscopy of two outflows will be discussed in detail in chapters 3.3 and 3.4, and some results on the embedded IRS 2 cluster will be presented in chapter 3.5.

## 2 Observations and Data Reduction

The data for this study of outflows in W51 were obtained in three observing runs. The initial wide-field imaging data in H, H2 1-0 S(1) at 2.122 $`\mu `$m, and K were obtained in 1997 and formed the basis for identifying previously unknown outflows. Some of the brightest outflows discovered in this survey were subsequently imaged at higher spatial resolutions in S(1) and \[FeII\] in June 2000. High resolution spectroscopy of two outflows was obtained at Keck in May 2000 using NIRSPEC (McLean et al., 1998).

### 2.1 Wide-Field H, S(1), and K-band Images

Our search for emission line objects in W51 is based on a wide-field mosaic obtained in the broad band H and K filters and in a narrow-band filter centered on the H2 1–0 S(1) line at 2.122 $`\mu `$m. The data were obtained with the infrared camera QUIRC (Hodapp et al., 1996) on 1997 August 16 – 19 (UT) at the UH 2.2 m telescope f/10 focus. The pixel scale was 0$`\stackrel{}{\mathrm{.}}`$188 per pixel, resulting in a field of view of 193$`\mathrm{}\times `$193$`\mathrm{}`$. To cover the large area of the W51 molecular cloud in a reasonable time, we used 10 s integration time in the H and K filters and 20 s in the molecular hydrogen S(1) narrow-band filter. The imaging data were obtained in a 10$`\times `$9 mosaicing pattern with 90$`\mathrm{}`$ spacing, so that every position on the sky is covered by at least 4 individual exposures. The total integration time for each point in the image is therefore 40 s in the broad-band filters and 80 s in the S(1) narrow-band filter.

The individual frames were flatfielded using domeflats and known bad pixels were masked off. The individual reduced frames were assembled into the large mosaic by aligning stars in the overlap region between of two frames. The full K-band mosaic was assembled first, and then used as a template to position the individual frames obtained in the other filters. This method ensures that all three colors are well registered relative to each other, so that the three mosaic images can be combined into a false-color image.

The three separate images in the three filters were then combined into a false-color image, with H represented in blue, S(1) in green, and K in red. Regions of S(1) line emission stand out as greenish features in this color scheme. A compressed version of this false-color image can be downloaded from http://www.ifa.hawaii.edu/$$hodapp/W51-sm-rgb.jpg. Visually identified potential shock front features were checked on the individual H, S(1), and K mosaics and false detections due to residual images, ghost images, or glare, as well as bad pixels or cosmic ray hits were discarded.

The extended emission permeating much of the survey area has been shown by Okumura et al. (2001) to be dominated by Br$`\gamma `$ and HeI emission, but to also contain many emission lines of H2, probably fluorescently excited by the UV field from the ionizing stars in the HII regions. In identifying shock-excited H2 emission in the survey images, the distinction between shock-excited features and features excited by fluorescence in the stellar UV field can only be based on the strength of the H2 emission, i.e. its contrast relative to the broad-band images, and its morphology. Objects identified as shock fronts in our survey have a strong contrast to the broad-band filters, indicating that their emission in the K band is dominated by H2 emission. They are either marginally resolved small knots or they clearly show the mophological features of bow shocks. In contrast, extended, rather diffuse features showing some H2 emission, but not being dominated by it, that are also mostly located on the perimeter of the HII regions outlined in broad-band K and the Br$`\gamma `$ image of Okumura et al. (2000), were not identified as potential stellar outflow shock fronts. These latter features are most likely excited by fluorescence in the UV radiation field generated by the ionizing stars in the HII regions, as first suggested by Okumura et al. (2001).

The extent of the survey area and the location of detected S(1) emission features fitting our search criteria is indicated in Fig. 1, the K-band image produced in our survey. The individual emission knots, or systems of knots in cases of clearly related (bipolar) outflow features, are shown in Fig. 2. For each object, the name follows the IAU conventions and includes the J2000 coordinates. Usually, the coordinates refer to the center of the extended emission region. In highly structured, irregular emission knots, the coordinates refer to the most prominent feature. For systems of emission features that morphologically appear to be part of one outflow, the coordinates refer to one of the features, and additional features are referred to by additional, self-explanatory designations (N, S, W, E, or C).

### 2.2 High Resolution Emission Line Imaging

To obtain images at higher spatial resolution of some of the likely outflow shock fronts discovered in the wide-field survey, tip-tilt corrected images were obtained on June 16 - 19, 2000 (UT). The tip-tilt corrected f/31 focus (Jim et al., 2000) of the UH 2.2m telescope was used where the QUIRC infrared camera has a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$062 per pixel and a field of view of 63$`\mathrm{}\times `$63$`\mathrm{}`$. Images were obtained in narrow-band \[FeII\] 1.644$`\mu `$m and H2 1–0 S(1) 2.122$`\mu `$m filters and the data were reduced in the same way as described above for the wide-field images. Under good seeing conditions, the tip-tilt corrected images at 2.122$`\mu `$m come close to the diffraction limit of the telescope and under more typical conditions, 0.3$`\mathrm{}`$ FWHM is achieved.

### 2.3 Infrared Spectroscopy

High resolution echelle spectroscopy of two of the outflow sources in W51 was obtained in the night of May 28, 2000 (UT) at the W.M. Keck Observatory, using the NIRSPEC (McLean et al., 1998) spectrograph. The weather in this night was non-photometric and unsuitable for flux calibration of the spectra.

The spectra were obtained with two grating and filter settings: The first setting was for observations in the K band, intended to cover several of the H2 shock-excited emission lines but concentrating on the 1-0 S(1) line at $`\lambda _{vac}`$=2.1218334$`\mu `$m (Bragg, Smith, & Brault, 1982). Observations in Br$`\gamma `$ were not part of the original observing plan, since this line is rarely excited in H2 outflows. Fortunately, the Br$`\gamma `$ line at $`\lambda _{vac}=2.166167\mu `$m was also recorded in our K-band spectra, albeit near the edge of the detector array so that wavelength calibration was difficult. The other grating setting was for observations in the H band, primarily to observe the \[FeII\] line at $`\lambda _{vac}`$=1.64400$`\mu `$m (Johannson, 1978) .

Since the objects are extended, sky spectra were taken separately at positions away from the visible shock fronts, but still within the extended emission from the HII regions permeating much of the W51 star forming region. The raw data frames were flatfielded with frames taken with continuum illumination of the spectrograph slit. We did not subtract dark current frames from the individual exposures, but relied on the sky frames to subtract dark current.

Wavelength calibration was done using OH airglow lines recorded during the on-object integrations. The line positions were measured on the flatfielded frames before sky subtraction. In the H band, we used the OH airglow lines (wavelengths from the UKIRT website) at 1.63885$`\mu `$m, 1.64147$`\mu `$m, 1.64476$`\mu `$m, and 1.65023$`\mu `$m for wavelength calibration. In the K band, OH lines at 2.11766$`\mu `$m, 2.12325$`\mu `$m, and 2.12497$`\mu `$m were used for wavelength calibration near the H2 1–0 S(1) line. As mentioned above, Br$`\gamma `$ was not part of the original observing plan, but was nevertheless recorded near the edge of the detector array. Wavelength calibration for Br$`\gamma `$ is based on extrapolation from the two OH lines at 2.17111$`\mu `$m and 2.18022$`\mu `$m. Unfortunately, we did not record any cataloged OH lines at wavelengths shorter than Br$`\gamma `$ in the same echelle order, limiting the accuracy of the wavelength calibration. Only the \[FeII\] line, H2 1–0 S(1) and Br$`\gamma `$ were evaluated. Other transitions of H2 were detected in the spectrum, but not with sufficient signal to noise to allow detailed analysis. The slit width corresponds to a velocity resolution of 22 km/s.

In W51, Okumura et al. (2001) have discussed near-infrared spectroscopy of the two compact HII regions W51 IRS 2 East and West. They conclude that the faint S(1) line emission in or around the HII region is excited by fluorescence rather than shocks. This diffuse molecular hydrogen emission is distinct from the shock-excited emission discussed here for the outflow objects. However, the presence of UV-excited fluorescent S(1) line emission gives a background level of line emission that is clearly detected in our spectra and not always fully removed by the sky-subtraction process due to its spatial variation.

## 3 Results and Discussion

### 3.1 Distribution of outflow sources

The locations of the shock-excited H2 emission features are indicated on the K-band image of the wide-field survey data (Fig. 1). One group of outflows is located in dense filaments near the main HII regions associated with IRS 1 and IRS 2. The others are distributed in almost a circular pattern at larger distances of about 6$`\mathrm{}`$ away from this central group, on the perimeter of the HII region complex. This distribution of outflow sources may be partly explained by an observational selection effect, since our detection method is clearly not as sensitive in regions superposed on or lying behind high levels of extended flux. More importantly, however, shock excited H2 emission cannot occur in fully ionized regions. In the likely case that the distribution of outflows is not entirely dominated by selection effects, it suggests that the class 0 and class I sources responsible for these outflows may have formed as a result of secondary star formation triggered by winds from the central HII regions and the cluster around IRS2.

### 3.2 Individual Molecular Outflows and Emission Knots

The individual knots of H2 emission found in this survey are named following the IAU convention: The name indicates their association with the W51 complex and their discovery by emission of the H2 molecule and gives their J2000 coordinates with 1$`\mathrm{}`$ precision. In the following, we will discuss the morphology of each candidate emission knot in detail.

#### 3.2.1 W51H2 J192318.1+142738

The H2 emission W51H2 J192318.1+142738 is a single, slightly elongated patch of S(1) emission, about 1.5$`\mathrm{}`$ in NS extent (Fig. 2a).

#### 3.2.2 W51H2 J192318.5+142656 and W51H2 J192319.6+142654

The emission knots W51H2 J192318.5+142656 and W51H2 J192319.6+142654 are only separated by about $``$ 10$`\mathrm{}`$. We show both of them in the high-resolution S(1) image Fig. 3 and discuss their relationship below.

Many H2 outflows in nearby star-forming regions are associated with class 0 sources that are not directly visible at near-infrared wavelengths. Since mass accretion and the associated outflow power peak during the class 0 phase, these extremely young sources often produce the strongest S(1) shock fronts. It would therefore be entirely plausible that the two shock fronts W51H2 J192318.5+142656 and W51H2 J192319.6+142654 are part of the same bipolar outflow, even though there is no positive detection of a central star that would drive the outflow.

The arguments against this scenario are based on the morphology of other extended emission features in the immediate vicinity of these two brightest emission knots. The star $``$ 9$`\mathrm{}`$ north of W51H2 J192318.5+142656 is the reddest object within about $`1\mathrm{}`$ and shows some reflection nebulosity (based the neutral color of the extended emission in the H, H2 S(1) and K image) extending in the general direction of the W51H2 J192318.5+142656 emission knot. It is therefore a plausible alternative scenario that W51H2 J192318.5+142656 originates in an outflow from that red star (probably a class I source). The emission knot W51H2 J192319.6+142654 shows extended emission north and bending towards east from the main emission knot. This emission is faintly indicated in Fig. 2a, showing that it is mainly S(1) emission, and is clearly shown in the high-resolution image in Fig. 3. Based on these two morphological arguments, we list the two emission knots as separate objects, probably belonging to different outflows.

#### 3.2.3 W51H2 J192322.0+143333

The S(1) line emission in W51H2 J192322.0+143333 appears marginally resolved in the high-resolution image in Fig. 4. The only other faint extended feature in this image is the elongated nebulosity 1$`\mathrm{}`$ east and 3.5$`\mathrm{}`$ south of W51H2 J192322.0+143333 that points roughly to the latter. On morphological grounds, it could be nebulosity associated with the driving source of W51H2 J192322.0+143333, but is too faint to be detected in the H, H2 S(1), and K survey images, so its color is unknown.

#### 3.2.4 W51H2 J192323.6+142515

The pair of emission regions W51H2 J192323.6+142515 S and N are clearly part of the same bipolar outflow. The high-resolution S(1) image in Fig. 5 reveals a feature exactly in the middle of the bright N and S regions that consists of an unresolved source and some nebulosity, separated by about $`1\mathrm{}`$ along the axis between W51H2 J192323.6+142515 S and N, at coordinates 19:23:23.6 +14:25:15 (J2000). We label these features as W51H2 J192323.6+142515 C (center). Note that in this case the two shock fronts of this outflow are listed as the north and south features of an outflow listed under the coordinates of its central source. These central features are very red and are not visible at all in the H band wide-field image. They are weakly indicated in the K-band image. This is exactly the morphology seen in many closer bipolar outflow sources, where two shock fronts are seen equidistant to some scattered light near the central source. An example of this is the outflow in L 1634 (Hodapp & Ladd, 1995). For comparison, the extent of the W51H2 J192323.6+142515 outflow is $`30\mathrm{}`$ between the two shock fronts, while the L 1634 outflow in Orion at a distance of 500 pc is $`6.2\mathrm{}`$ long, nearly the same linear projected extent of $``$1 pc.

#### 3.2.5 W51H2 J192325.9+143703 and W51H2 J192327.9+143701

Two faint emission knots in close vicinity were found at the northern edge of our survey field. The two panels in Fig. 2b showing W51H2 J192325.9+143703 and W51H2 J192327.9+143701 actually overlap. There is no morphological indication that the two emission knots are part of the same outflow, so, conservatively, we list them separately, as was discussed above.

#### 3.2.6 W51H2 J192335.0+143028 and W51H2 J192336.6+143014

The emission knot W51H2 J192335.0+143028 is a bright, compact knot of H2 1–0 S(1) line emission. Close to it is a second emission knot W51H2 J192336.6+143014. It is noteworthy that these two emission knots are collinear with one of the reddest stars detected in our whole W51 survey area, labeled RS in Fig. 6, at 19:23:32.7 +14:30:48 J2000. The red star is undetected in our H-band image, and appears much brighter in K (2.2$`\mu `$m) than in S(1) (2.12$`\mu `$m), compared to other stars in the field. Fig. 6 is a larger subframe of the survey images and illustrates this collinear relationship between these features. This alignment suggests that the red star may be the driving source of a well-collimated jet of $``$ 2.2 pc projected length that produces the two emission knots. The degree of collimation of this jet would be quite remarkable, but not unprecedented.

#### 3.2.7 W51H2 J192338.3+143047

The faint emission knot W51H2 J192338.3+143047 has no morphological features (Fig. 2c) that would help in identifying its driving source.

#### 3.2.8 W51H2 J192339.7+143131

The largest shock-front detected in S(1) is W51H2 J192339.7+143131. It appears morphologically to be associated with a driving source in the general direction of the IRS 2 young cluster to its south. A plausible candidate for the driving source is the extremely red star at 19:23:39.8 +14:31:21 that lies close to the symmetry axis of the bow shock fronts. More detailed images in the emission lines of \[FeII\] and S(1) are shown in Fig. 7. The S(1) emission shows the typical shape of multiple bow shocks, while \[FeII\] is generally more concentrated on the axis of the outflow, i.e. the areas of highest shock velocity and excitation. It is noteworthy, however, that the S(1) emission exhibits pronounced peaks near the apexes of two of the bow shocks, different from what is found in most other S(1) bow shocks, e.g., (Tedds, Brand, & Burton, 1999). This object will be discussed in more detail below (3.4) in the context of our spectroscopic data.

#### 3.2.9 W51H2 J192345.5+143537

The two emission knots found in W51H2 J192345.5+143537 are treated as part of one system, since their close proximity makes a chance superposition of two unrelated shock fronts very unlikely. The high-resolution image in Fig. 8 shows one knot clearly resolved, the other marginally resolved, but no further conclusions about the location of the driving source can be derived from these images.

#### 3.2.10 W51H2 J192347.2+142944

The star labeled as W51H2 J192347.2+142944 stands out due to its strong flux in the S(1) filter. While the star appears slightly more extended than other stars on the S(1) frame and more extended than its K-band image, the distribution of S(1) emission around the star cannot be mapped out. We conclude that the molecular emission in this region originates in the immediate vicinity of the star.

#### 3.2.11 W51H2 J192403.3+143255

The emission knot W51H2 J192403.3+143255 is extended, but small. No obvious association with other emission knots or red stars was found that might help in identifying the driving source.

### 3.3 Velocity-resolved spectroscopy of W51H2 J192323.6+142515 S

A high spatial resolution image of the two H2 emission features W51H2 J192323.6+142515 N and S is shown in Fig. 5 and was discussed in section 3.2.4. A high-resolution H2 spectrum was obtained towards only the southern component at the slit position indicated in Fig. 5 and is shown in Fig. 9.

In knot S-A, farthest from the likely driving source, broad, double-peaked H2 profiles peaking at $`V_{\mathrm{LSR}}60`$ km s-1 (just blue of the cloud systemic velocity of 70 km s-1) and $`V_{\mathrm{LSR}}0`$ km s-1 (blue-shifted by $`70`$ km s-1) are observed. This velocity-split emission corresponds to two spatially separate emission knots in Fig. 5. Note however that the separation of the two velocity components is much larger than the spectral resolution of 22 km/s, so the velocity split is real and not just a projection of two spatially separate knots in the slit. Behind (north of) the shock front S-A (offset 1″ in Fig. 9) the H2 velocities converge to an intermediate, blueshifted velocity $`V_{\mathrm{LSR}}20`$ km s-1, spatially coinciding with a single emission knot centered on the slit.

The fainter knot S-B (not covered well by the slit) also shows a double-peaked H2 profile with peaks at $`V_{\mathrm{LSR}}40`$ km s-1 and $`V_{\mathrm{LSR}}20`$ km s-1, i.e. about 20 km s-1 more blueshifted than knot S-C. Further back towards the driving source, shock front S-C is much more diffuse than S-A and S-B, and fainter. Fig. 9 shows a broad line centered at $`V_{\mathrm{LSR}}30`$ km s-1, without a clear indication of a separation into two peaks.

The observed double-peaked H2 profiles, with peak-to-peak separations of $`60`$ km s-1, are predicted by numerical models (Völker et al., 1999) and are implied by the analytical (Hartigan, Raymond, & Hartmann, 1987) J-type bow shock models of atomic line emission. The most extended, double-peaked profiles will be expected near the front of the bow shock; narrow, low-velocity peaks will instead be associated with the oblique wings. Indeed, the range of H2 velocities observed in W51H2 J192323.6+142515 S can be explained on purely geometrical grounds, if the flow is inclined towards the observer at an angle (with respect to the line of sight) of $`\varphi =40^{}70^{}`$. By comparison, a bow shock moving in the plane of the sky ($`\varphi =90^{}`$) would produce two (blended) peaks blue- and red-shifted by the same amount, while a bow viewed head-on would produce only one, blue-shifted component (see for example Plots IV and I in the appendix in Davis, Hodapp, & Desroches (2001).

Overall, the spectral features seen in W51H2 J192323.6+142515 S are quite similar to those found in HH 111 by Davis, Hodapp, & Desroches (2001). They can be explained by the straightforward geometric fact that the H2 emission arises in the oblique shocks in the wings of the bow shocks, where shock velocities are low despite a high velocity of the jet relative to the ambient cloud medium. The shock velocities seen in W51H2 J192323.6+142515 S are just above the dissociation limit for pure J shocks, which suggests a certain degree of magnetic cushioning to avoid rapid dissociation of the H2 molecule (Smith, M. D., 1994) .

### 3.4 Velocity-resolved spectroscopy of W51H2 J192339.7+143131

#### 3.4.1 Morphology and slit orientation

The bow shocks in W51H2 J192339.7+143131 indicate a much faster and powerful outflow than the outflow discussed above. High-spatial resolution images of W51H2 J192339.7+143131 in \[FeII\] and H2 1–0 S(1) are shown in Fig. 7. This object is situated 23$`\mathrm{}`$ north of the IRS 2 cluster.

In H2, W51H2 J192339.7+143131 resembles a sequence of at least three “nested” bow shocks separated by 2″–3″, labeled A, B and C in Fig.7. A very similar arrangement of nested bow shocks, with the larger ones being found in the wakes (i.e. towards the driving source) of smaller bow shocks, is observed in the low-mass L 1634 outflow (Hodapp & Ladd, 1995). The bow shock nature of W51H2 J192339.7+143131 is confirmed by high-resolution spectroscopy in the H2 1–0 S(1) and \[FeII\] lines (Fig. 10).

The \[FeII\] emission in W51H2 J192339.7+143131 (Fig. 7, left panel) appears to be generally more closely confined to the north-south outflow axis than the H2 emission. This is not unexpected since the \[FeII\] traces the higher-excitation bow shock caps while the H2 is excited in the oblique, lower-excitation bow-shock wings (Davis, Smith, & Eislöffel, 2000; Lorenzetti et al., 2001). However, contrary to this generalized statement, we find strong H2 emission from near the apexes of the bow shocks, which we believe to result from fluorescence, as will be discussed below.

We present \[FeII\], H2, and Br$`\gamma `$ (Fig. 10) spectra observed through the center of W51H2 J192339.7+143131. The slit was aligned to cover the brightest H2 1–0 S(1) emission knots (B and C) and is probably closely aligned with the outflow axis, if our identification of the driving source or at least the location of the driving source in the general area of IRS 2 is indeed correct. Millimeter-wave mapping of the cloud structure in W51A suggests that the LSR systemic velocity of the IRS 2 region is $`61`$ km s-1 (Carpenter and Sanders, 1998), although variations in the systemic velocities of cloud cores across the region are evident in the molecular cloud maps. Spatially extended S(1) emission around IRS 2, probably UV-excited in the radiation field of the 0 stars in IRS 2, and largely subtracted out from our sky-subtracted spectral images, is centered around 70 km/s. All three spectra were individually wavelength calibrated using OH airglow lines recorded on the object frames. The wavelength calibration of the Br$`\gamma `$ line is relatively poor, since Br$`\gamma `$ was not part of the original observing plan and was only evaluated when interesting structure was unexpectedly found. As a consequence, Br$`\gamma `$ was recorded near the edge of the detector array where optical distortions could be significant, and the wavelength calibration could only be done by extrapolation from OH-lines, not by interpolation, since there were no OH lines recorded shortward of 2.166 $`\mu `$m in the echelle order used. The small differences in the velocity of the narrow components of S(1) and Br$`\gamma `$ emission are probably due to these calibration problems and we do not believe they are significant. We adopt a systemic cloud velocity of $`70`$ km s-1 for IRS2 and W51H2 J192339.7+143131.

The position-velocity images in Figs. 10 show very distinct features in the three emission lines studied here. The H2 1–0 S(1) emission is largely confined to a narrow velocity range around the systemic velocity. The velocity structure seen in H2 1–0 S(1) is, however, correlated to the high-velocity features seen in \[FeII\] and Br$`\gamma `$, in the sense that knot B shows a slight blueshift and knot C a slight redshift. Only in the faint knot A, north of the main shock fronts, does the S(1) line split, producing a weak blueshifted component, as expected for the spectrum of a non-dissociative bow shock.

The \[FeII\] line is a more robust tracer of shocked gas than the rather fragile H2 molecule and more closely traces the full velocity field in the shock fronts, which in this case extends over a range of 300 kms-1. The detection of velocity features related to the shock fronts in Br$`\gamma `$ was a surprise. We had, of course, expected copious Br$`\gamma `$ emission from the nearby IRS2 HII regions, either projected as foreground or background emission or scattered into the line of sight. However, the Br$`\gamma `$ line in W51H2 J192339.7+143131 shows features clearly correlated to those seen in the other two lines studied here. These features are clearly related to the outflow shock fronts and make W51H2 J192339.7+143131 a rare case of outflow shock fronts with detectable Br$`\gamma `$ emission. We will now discuss the individual emission lines in order of wavelength.

#### 3.4.2 The \[FeII\] Line

The \[FeII\] emission is concentrated in the region between the H2 knots B and C and is less spatially extended perpendicular to the jet axis than the S(1) emission (Fig. 7). This is consistent with the fact that \[FeII\] emission should be concentrated in the high excitation regions directly at the head of the bow shock, a region where S(1) emission is suppressed by dissociation of the H2 molecule.

At the position of knot C and up to 1.5″north of it, the \[FeII\] emission is strongly redshifted, and very broad, extending from $`V_{\mathrm{LSR}}20`$ km s-1 to $`V_{\mathrm{LSR}}200`$ km s-1. Further north, up to the position of knot B, broad emission is observed with velocities ranging from blueshifted $`V_{\mathrm{LSR}}60`$ km s-1 to redshifted $`V_{\mathrm{LSR}}120`$ km s-1. Overplotted on the spectral image in Fig. 10 is the normalized spectrum integrated over the relevant parts of the slit, for comparison with the model presented by Hartigan, Raymond, & Hartmann (1987). The broad, only slightly asymmetric profile measured in W51H2 J192339.7+143131 closely matches the theoretical profile calculated for H$`\alpha `$ emission of their 200 km/s jet model inclined by 60 against the line of sight towards the observer (their Fig. 3c). In our case, the total velocity spread is larger, about 300 km/s and the redshifted component is somewhat stronger and more extended than the model would predict, but the overall agreement to the idealized model is remarkably good.

The \[FeII\] spectrum shows a much wider spread of velocities than the S(1) spectra at the same slit position. Our S(1) and \[FeII\] images suggest a spatial anticorrelation between the two emission lines. S(1) is primarily emitted in front and behind the main shock outlined in \[FeII\], since the high temperatures in this 300 km/s shock dissociate the H2 molecule.

In addition to the good agreement of the integrated \[FeII\] spectrum with the models by Hartigan, Raymond, & Hartmann (1987), the details of the \[FeII\] spectrum, i.e., the wide velocity range observed and the prominent redshifted emission can be explained by a fully developed hydrodynamical model.

The hydrodynamical models of Völker et al. (1999) discuss the gas motions near the apex of a jet in detail. Their Fig. 13 clearly shows that, in the reference frame of the jet, gas is being pushed sideways at the apex of the jet and streams backward in the wings of the bow shock. The same models also show that the velocity field in the bow of a shock can be complex and that knots of peculiar velocities exist. In the case of W51H2 J192339.7+143131 where the predominantly blueshifted and redshifted components of \[FeII\] emission are spatially separated, we have to postulate that the emission arises from distinct emission knots near the apex of the bow, the redshifted emission being on the far side of the jet and therefore receding relative to the jet, while the blueshifted emission is from the near side of the bow, moving towards the observer in excess of the jet velocity.

#### 3.4.3 The H2 1–0 S(1) Line

The H2 1–0 S(1) line in W51H2 J192339.7+143131 is strikingly narrow in an object that shows such substantial line width in \[FeII\]. Only the northernmost and faintest of the emission knots identified in Fig. 7 (knot A) shows the split line profile expected from a bow shock seen from the side. The presence of \[FeII\] emission supports J-type shock excitation, i.e., no significant cushioning by magnetic fields. In such shocks, H2 will be dissociated at relatively low shock velocities, of the order of 20-25 kms-1. The centers of the two velocity components in knot A are separated by about 50 kms-1, which is roughly double the dissociation shock velocity limit for H2. The other two knots (B and C) do not show a splitting of the lines, but show slightly broadened spectral features spreading over a maximum of 50 kms-1, i.e. within the dissociation limit. The velocity structure is related to the velocity field seen in \[FeII\] emission in the sense that knot B exhibits a slight blueshift while knot C is slightly redshifted. H2 will survive and be collisionally excited into emission only in the oblique wings of the bow shocks, exhibiting narrow line profiles. An example of this are the bullets in Orion (Tedds, Brand, & Burton, 1999) where the tips of the bows are traced in \[FeII\] and the wings in H2 S(1). However, in our images the H2 emission appears brightest on-axis, in the two compact knots B and C near the apex of the bow shocks, where molecules will be collisionally dissociated. We therefore postulate that the strong, spatially concentrated H2 emission near the apex of the bow with narrow line width represents fluorescent emission from quiescent molecular material in front of the bow shock that is excited by the UV radiation field generated by the high temperature in the shock front. The presence of Br$`\gamma `$ emission in the shock front is proof that temperatures high enough to dissociate H2 and ionize H exist locally in the shock front. Fluorescent emission from H2 in the UV has been found in IUE spectra of several HH objects by Böhm, Scott, & Solf (1991). In an object similar to W51H2 J192339.7+143131, the well studied HH 7, the higher excitation lines of H2 have been explained by fluorescence of H2 in the presence of a strong UV field generated in the shock fronts (Fernandes & Brand, 1995).

#### 3.4.4 The Br$`\gamma `$ Line

Emission in the Br$`\gamma `$ line is very rarely seen in the shock fronts of stellar outflows, even though H$`\alpha `$ is commonly observed in low-extinction Herbig-Haro jets and shock fronts (see Reipurth & Bally (2001) for a recent review). Br$`\alpha `$ has recently been observed (Fuller, Zijlstra, & Williams, 2001) in an outflow from a high-mass star.

Emission of Br$`\gamma `$ is, of course, expected near the HII region associated with IRS2 and could arise in the foreground or background of the outflow shock fronts of W51H2 J192339.7+143131, or be scattered into that line of sight. Remarkably, however, the Br$`\gamma `$ emission in W51H2 J192339.7+143131 shows features similar to those seen in \[FeII\] and H2, demonstrating that a significant fraction of the Br$`\gamma `$ flux originates in the outflow shock fronts themselves.

The brightest peak in the Br$`\gamma `$ position-velocity image is just north of knot B. The narrow-line Br$`\gamma `$ intensity south of knot B drops to 55% of its value (per pixel) just north of this knot. This indicates that $``$45% of this narrow-line flux must originate close to the shock front, while $``$55% may be foreground or background flux from elsewhere along the line of sight. The Br$`\gamma `$ emission in knot C is more spatially extended than the H2 emission, but they are clearly related.

At very low but clearly significant flux levels, spectrally very broad Br$`\gamma `$ emission essentially duplicates the features seen in \[FeII\]: strong blueshift in knot B and redshift in knot C, with about the same total velocity extent. This component of the Br$`\gamma `$ emission is clearly related to fast moving gas in the bow shock front. The broad high-velocity Br$`\gamma `$ emission is probably collisionally excited in the rapid, high-excitation shock fronts, similar to the \[FeII\] emission.

As explained above, absolute velocity calibration of the Br$`\gamma `$ is rather uncertain due observational limitations. We therefore do not ascribe any significance to the absolute velocity difference between the H2 1–0 S(1) and Br$`\gamma `$ lines in Fig. 10.

#### 3.4.5 Synopsis of W51H2 J192339.7+143131

Of the three emission lines discussed here, \[FeII\] is a pure and robust tracer of shock-excited gas. It shows a spread of velocities over a range of 300 kms-1 and its integrated line profile matches simple models of bow shock emission very well. The Br$`\gamma `$ line traces both the high-velocity gas, either by direct shock excitation or UV excitation from the shock, but also traces the low velocity components, by UV excitation from the shock in addition to fore and background emission. The H2 emission traces only low-velocity gas, since H2 dissociates in high velocity shocks. The strong H2 emission from the tips of the bow shocks can best be explained by fluorescent excitation of H2 molecules in the radiation field of the shock, just prior to their being dissociated when hit by the shock directly.

An outflow system appearing superficially similar to W51H2 J192339.7+143131 are the “bullets” observed to emerge from a source in the Orion-Trapezium cluster. A detailed study of this system, including high-resolution spectroscopy in the \[FeII\] and H2 S(1) lines has recently been presented by Tedds, Brand, & Burton (1999). The shock fronts in W51H2 J192339.7+143131 appear to have a more organized large scale shape than the system of small ”bullets” in Orion, but some small-scale clumpy structure is clearly present in W51H2 J192339.7+143131. The Orion ”bullets” do not show the very broad \[FeII\] profiles seen in W51H2 J192339.7+143131 and they do not exhibit the strong H2 emission at the apex of the bow shock. Rather, they show H2 only in the wings of the bow shocks, as expected from pure shock excitation.

Observations in \[FeII\] of the jets emanating from L1551 IRS5 have recently been reported by Pyo et al. (2002). Their spectra along the jet axis show a similarly broad velocity profile of the \[FeII\] line to the one reported here. In contrast to our observations, theirs concentrated on the jet close to its source of origin and did not include well developed bow shocks.

### 3.5 The young stellar cluster around IRS2

The narrow-band H2 1–0 S(1) line image was also used to study the young cluster surrounding IRS 2 in more detail. This image reaches about the same limiting magnitude and has better spatial resolution than the K-band image obtained of the wider W51 region, and is therefore best suited for a rough count of the stars contained in this cluster. The S(1) line image was photometrically calibrated using the standard star FS 27 (Hawarden et al., 2001) so that approximate K-band magnitudes could be obtained. Down to a limiting magnitude of K=15.3, we count 60 stars within a projected circle of 1 pc diameter. Okumura et al. (2000) have already pointed out that one foreground star is visible in front of the W51 IRS 2 cluster, but clearly most of the stars counted around IRS 2 are physically close to this embedded source. The image (Fig. 11) also shows areas very near IRS 2 that appear as dark patches against the extended flux of the HII region and with virtually no stars, indicating opaque extinction. Roughly half of the 1 pc diameter circle is covered by such extinction.

The distance of 7Kpc to W51, compared to 0.5 Kpc to the Trapezium cluster, makes stars appear 5.73 magnitudes fainter. Okumura et al. (2000) found a strong peak at AV=25 mag in their extinction histogram of ”region 3” that contains the young embedded cluster associated with IRS2. We take AV=25 here as an estimate of the extinction to those stars in the cluster that are not completely obscured, obviously an extreme simplification of the actual situation. Similar AV values were found for many of the stars in the immediate vicinity of IRS2 by Goldader & Wynn-Williams (1994). While infrared excess may be partly responsible for the red colors leading to these extinction estimates, we assume here that stars in the W51 IRS2 cluster suffer about AK=2 mag more extinction than the low-extinction stars seen in the Orion Trapezium (Herbig & Terndrup, 1986). This statement obviously excludes the high extinction objects in the BN/KL region. Combining the effects of distance modulus and extinction, stars in the W51 IRS2 cluster appear about 7.7 mag fainter than they would in the Orion Trapezium cluster. Our limiting K-band magnitude of 15.3 in the high-resolution image in Fig. 12 corresponds to a magnitude of 7.6 in Orion. On the 2MASS K-band frame we count 15 stars brighter than K=7.6 mag in a circle of projected diameter of 1 pc, using photometry in the Trapezium cluster from McCaughrean & Stauffer (1994).

In comparing the young stellar clusters associated with the Orion Trapezium and W51 IRS 2, we can reach a few tentative conclusions, despite the uncertainties in the estimates for extinction and the incompleteness of the star count due to opaque areas obscuring about half the cluster, if one assumes a roughly spherical intrinsic distribution of the stars. The star count in the W51 IRS2 cluster appears between 4 to 8 times higher than in the Trapezium cluster. This is roughly consistent with the count of O stars by Goldader & Wynn-Williams (1994) (4 in Trapezium vs. probably 9 in W 51 IRS 2).

## 4 Conclusions

Our H, H2 1–0 S(1), and K-band survey of the W51 GMC resulted in the discovery of 14 H2 shock fronts associated with stellar outflows. The outflows are found in dense molecular filaments near the central HII regions and the young embedded cluster associated with IRS 2, and near the perimeter of the molecular cloud complex, away from the HII regions. We speculate that the outflows near the perimeter represent a secondary, triggered phase of star formation in W51.

Detailed, high-resolution images were used to identify plausible candidates for the driving sources of some of the outflows.

For two of the outflow shock fronts, high resolution spectroscopy was obtained. Of particular interest are the shock fronts W51H2 J192339.7+14313. The very broad \[FeII\] line seen in this shock front is in good agreement with models developed by Hartigan, Raymond, & Hartmann (1987) for other atomic lines in shock fronts. Refinements of these models based on the hydrodynamic models of Völker et al. (1999) lead to a satisfactory agreement with the observations, even though the large observed line width and the large redshifted component remain quite remarkable. We also report the very rare detection of high-velocity Br$`\gamma `$ emission from this outflow, with a high-velocity components essentially matching the \[FeII\] profiles. The strong, narrow H2 1–0 S(1) emission found near the apex of the shock fronts strongly suggests fluorescence as the excitation mechanism for H2 at this location. The outflow W51H2 J192339.7+14313 is driven by a source associated with the IRS 2 young embedded cluster. Based on star counts and rough extinction estimates, we conclude that this cluster is more massive, richer in stars and richer in high-mass stars than the Orion-Trapezium cluster, in agreement with previous results by others.

The imaging data presented here were obtained at the University of Hawaii 2.2m telescope. The high-resolution spectroscopy presented here is based on data obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to extend special thanks to those of Hawaiian ancestry on whose sacred mountain we are privileged to be guests.

Figure Captions


format: text

text:
# Effect of an impulsive force on vortices in a rotating Bose-Einstein condensate

## I Introduction

Since the successful detection 1  of Bose-Einstein condensates (BEC) in dilute trapped bosonic atoms at ultra-low temperature, one problem of extreme interest is the formation of vortices in a rotating condensate in an axially symmetric trap. The experimental observation exp ; exp1 ; expm ; jila ; hal ; hal1 ; mit ; mit1  of vortices in the condensate has intensified theoretical investigations on vortices using the mean-field Gross-Pitaevskii (GP) equation gp . For large frequencies of rotation, a very large number of vortices arranged in a triangular lattice, frequently referred to as a vortex lattice, have been observed. The vortex lattices are highly excited collective states of BECs with an angular momentum of up to $`60\mathrm{}`$ per particle mit . There have been many theoretical studies on different aspects of BEC 2 ; 3 ; 11  and specially, on the vortices in a rotating BEC in axially symmetric traps 2a ; 2c ; 2d ; 2e ; sin ; 2f ; 2ff ; 2f1 ; 2g ; 2h ; fc ; ad ; ueda ; dis .

When an ordinary fluid is placed in a rotating container the fluid rotates with the container because of viscous force. For superfluid 4He(II) in a rotating container no motion of the fluid is observed below a critical rotational frequency. Above the critical frequency, lines of singularity appear in its velocity field. These lines of singularity correspond to vortices and appear in the form of a lattice specially for large rotational frequencies. Superfluidity exhibits many peculiar features: absence of viscosity, reduction of moment of inertia, occurrence of persistent current, appearance of quantized vortices among others 11 . The formation of vortex lattice in a rotating 4He(II) is intimately related to its superfluidity and is a manifestation of quantized vortices exp ; exp1 ; expm ; jila ; mit ; mit1 . However, because of the strong interaction between 4He atoms, the theoretical description of this system is not an easy task.

Similar vortices can be generated in theoretical mean-field models of trapped BEC based on the GP equations 2a ; 2c ; 2d ; 2e ; sin ; 2f ; 2ff ; 2f1 ; 2g ; 2h ; fc ; ad ; ueda  and have been observed experimentally exp ; exp1 ; expm ; jila ; hal ; hal1 ; mit ; mit1  in BEC with repulsive interaction. In contrast to liquid 4He(II), a trapped BEC is a very dilute and weakly interacting system, which makes a mean-field analysis appropriate. This analogy of a trapped BEC with liquid 4He(II) suggests the presence of superfluidity in the BEC. Many different ways for generating quantized vortices in a BEC have been suggested 2d  , e.g., via spontaneous vortex formation in evaporative cooling 2f , via controlled excitation to an excited state of atoms 2ff , by stirring a BEC using an external laser with an angular frequency above a critical value for generating a quantized vortex with a single quantum 2c , or by the rotation of an axially symmetric trap with an angular frequency above a similar critical value 2g .

Experimentally, quantized vortices have been generated by different groups in a rotating BEC using a variety of techniques. Vortices have been detected in a 87Rb condensate in a cylindrical trap by Madison et al at ENS exp , where angular momentum is generated by a stirring laser beam. At MIT Raman et al expm  also studied the nucleation of vortices in a Bose-Einstein condensate stirred by a laser beam. They have been observed by Matthews et al jila  at JILA in coupled BEC’s comprised of two spin states of 87Rb in a spherical trap, where angular momentum is generated by a controlled excitation of the atoms between the two states. The tilting of the vortex axis in a nonspherical trapping potential has been studied at JILA by Haljan et al hal . Vortices have also been generated by evaporatively spinning up a normal gas of 87Rb and then cooling below quantum degeneracy at JILA by Haljan et al hal1 . A large number of vortices in the form of a lattice have been observed at MIT by Abo-Shaeer et al expm ; mit  as well as at ENS by Madison et al exp  for sufficiently large rotational frequency. The interesting dynamics of the formation and decay of vortices has also been studied mit1 . Madison et al exp  also confirmed that before a vortex pattern is formed the condensate passes through dynamically unstable distorted configurations. The dynamical instability transforms the distorted condensate to an appropriately symmetric condensate with vortices sin ; ueda .

There have been theoretical studies of how vortices are generated in a BEC via these dynamically unstable configurations using the mean-field GP equations in two ueda  and three fc  dimensions. The theoretical study of the dynamical instability using mean-field approach in the generation of vortex is of interest as one can then compare the results of the models with experiments. This will provide a more stringent test of the mean-field theories than the comparison of stable configurations of BEC.

However, there is another type of instability which has been studied in relation to a nonrotating BEC under the application of an impulsive force. One can now suddenly change the trapping potential on a stable BEC. This can be done in experiment by varying the current in the coil responsible for confinement. One can also change the strength of interatomic interaction suddenly by changing the magnetic field near a Feshbach resonance yy1 . These set in dynamical instabilities in the BEC, which can be studied both experimentally ex  and theoretically th .

Similar impulsive force can be applied on a rotating BEC with vortices and the resultant oscillation studied. One can study how the number of vortices change after the application of the impulsive force on a rotating BEC. In this paper we perform such mean-field analyses of the dynamics of a rotating BEC subject to three types of impulsive forces. We consider a sudden change in (i) the strength of interaction near a Feshbach resonance yy1  induced by a variation of an external magnetic field, (ii) the trapping potential, and (iii) the rotation of the BEC. Specifically, we study the dynamics when the strength of interaction or the trapping potential is suddenly increased or reduced. We also study the decay of vortices once the rotation of the condensate is suddenly stopped. There have been certain experimental measurements on such decay of vortices exp ; mit1 . Finally, the trapping potential can be removed suddenly in a rotating BEC with vortices and the resultant expansion studied.

We base our study on the numerical solution 3 ; 2a ; 2h  of the nonlinear time-dependent mean-field GP equation gp , which should provide a faithful description of the formation and evolution of vortices. The experiments on the vortices used a highly anisotropic trap which may be simulated by a quasi two-dimensional model. Hence as in the recent studies on the subject in Refs. ad ; ueda  we employ the GP equation in two dimensions. We use a split-step method where in the time evolution of the GP equation the kinetic energy and angular momentum (derivative) terms in the two directions are dealt with in two independent steps by the Crank-Nicholson rule complimented by the known boundary conditions 2h ; koo . All the harmonic oscillator and the nonlinear terms are treated in a third separate step. We find that this time-dependent split-step approach leads to good convergence.

In Sec. II we describe briefly the time-dependent GP equation for a rotating BEC and a numerical method for its solution. In Sec. III we report the numerical results of the present investigation and finally, in Sec. IV we give a summary of our study.

## II Nonlinear Gross-Pitaevskii Equation

At zero temperature, the time-dependent Bose-Einstein condensate wave function $`\mathrm{\Psi }(𝐫;\tau )`$ in two dimensions at position $`𝐫`$ and time $`\tau `$ may be described by the following mean-field nonlinear GP equation gp ; 11

$`[{\displaystyle \frac{\mathrm{}^2}{2m}}^2`$ $`+`$ $`V(𝐫)+gN|\mathrm{\Psi }(𝐫;\tau )|^2\widehat{\mathrm{\Omega }}L_{\widehat{z}}`$ (1)
$``$ $`i\mathrm{}{\displaystyle \frac{}{\tau }}]\mathrm{\Psi }(𝐫;\tau )=0.`$
Here $`m`$ is the mass and $`N`$ the number of atoms in the condensate, $`g`$ the strength of interatomic interaction. The trapping potential in two dimensions may be written as $`V(𝐫)=\frac{1}{2}m\omega ^2(\widehat{x}^2+\widehat{y}^2)`$ where $`\omega `$ is the angular frequency. The rotational term $`\widehat{\mathrm{\Omega }}L_{\widehat{z}}=i\mathrm{}\widehat{\mathrm{\Omega }}(\widehat{x}_{\widehat{y}}\widehat{y}_{\widehat{x}})`$ where $`\widehat{\mathrm{\Omega }}`$ is the rotational frequency along the $`\widehat{z}`$ direction with $`L_{\widehat{z}}`$ the $`\widehat{z}`$ component of angular momentum. The normalization condition of the wave function is $`𝑑𝐫|\mathrm{\Psi }(𝐫;\tau )|^2=1.`$ The GP equation (1) can accommodate quantized vortex states with rotational motion of the condensate around the $`\widehat{z}`$ axis.

It is convenient to use the dimensionless variables defined by $`x=\widehat{x}/l`$, $`y=\widehat{y}/l`$, $`t=\tau \omega /2,`$ $`l\sqrt{\mathrm{}/(m\omega )}`$, $`\mathrm{\Omega }=2\widehat{\mathrm{\Omega }}/\omega `$, and $`\phi (x,y;t)=l\mathrm{\Psi }(𝐫;t)`$. In these units Eq. (1) becomes

$`[{\displaystyle \frac{^2}{x^2}}`$ $``$ $`{\displaystyle \frac{^2}{y^2}}+\left(x^2+y^2\right)+i\mathrm{\Omega }\left(x{\displaystyle \frac{}{y}}y{\displaystyle \frac{}{x}}\right)`$ (2)
$`+`$ $`𝒩|\phi (x,y;t)|^2i{\displaystyle \frac{}{t}}]\phi (x,y;t)=0,`$
where $`𝒩=2gNm/\mathrm{}^2.`$ The normalization condition of the wave function is

$$_{\mathrm{}}^{\mathrm{}}𝑑x_{\mathrm{}}^{\mathrm{}}𝑑y|\phi (x,y;t)|^2=1.$$
(3)
The explicit inclusion of the rotational energy in the GP hamiltonian above simulates the time evolution of a rotational BEC which permits multiple-vortex formation. The initial input solution to the GP equation is a circularly symmetric state obtained with $`\mathrm{\Omega }=0`$. The method of solution with $`\mathrm{\Omega }0`$ is time iteration using an appropriate algorithm. If one iterates the GP equation above in time using the full angular momentum term the circular symmetry is maintained always and a solution with multiple vortices is not generated. It is essential that the circular symmetry of the solution is broken at some stage before multiple vortices are generated during time iteration.

The circular symmetry of the solution can be broken in the process of time evolution of the GP equation above by introducing fc  a phase in the initial wave function for $`\mathrm{\Omega }=0`$ of the form

$$\delta =\underset{i}{}\mathrm{arctan}[(yy_i)/(xx_i)],$$
(4)
where $`x_i`$ and $`y_i`$ are a set of parameters such that the positions $`(x_i,y_i)`$ fall inside the boundary of the BEC. This initial phase preserves the norm of the wave function but accelerates the formation of vortex nucleation after time evolution of the GP equation. However, there is no direct relation between the parameters $`(x_i,y_i)`$ and the positions or the numbers of the vortices in the final BEC. Another way of facilitating vortex formation is to include a small absorptive term in the GP equation as has been advocated in Ref. 2f  and employed in Ref. ueda .

We solve the GP equation (2) with $`\mathrm{\Omega }=0`$ using a time-iteration method elaborated in Refs. 2h ; koo  starting from the analytic two-dimensional harmonic-oscillator solution for $`𝒩=0`$. The full GP Hamiltonian is conveniently broken into three parts $`H_x`$, $`H_y`$ and $`H_{\text{res}}`$ $``$ the first containing the $`x`$-dependent derivative terms, the second containing the $`y`$-dependent derivative terms with the harmonic-oscillator and nonlinear interaction terms contributing to the third residual part. The GP equations along the $`x`$ and $`y`$ directions are discretized on a two-dimensional set of grid points $`N_x\times N_y`$ using the Crank-Nicholson discretization method koo . The second (first) derivatives are discretized using a three-point (two-point) finite difference rule. Effectively, each time iteration of the GP equation is broken up into three parts $``$ using $`H_x`$, $`H_y`$ and $`H_{\text{res}}`$. The first two parts lead to two sets of tridiagonal equations along $`x`$ and $`y`$ directions which are solved alternately by the Gaussian elimination method koo . The third part does not contain any derivative and is solved essentially exactly. For a small time step $`\mathrm{\Delta }`$ the error involved in this break-up procedure along $`x`$ and $`y`$ directions is quadratic in $`\mathrm{\Delta }`$ and hence can be neglected. Using the above procedure the solution of the GP equation corresponding to $`\mathrm{\Omega }=0`$ and $`𝒩0`$ is obtained. This solution serves as the initial input for the simulation of multiple vortex generation using the time iteration of the GP equation.

The initial solution with $`\mathrm{\Omega }=0`$ is circularly symmetric. For the generation of multiple vortex the circular symmetry is broken by introducing the phase $`\delta `$ of Eq. (4) to the solution with $`\mathrm{\Omega }=0`$ and iterating the GP equation in time with full $`\mathrm{\Omega }`$. We also included a small dissipative term in the GP equation. As the formation and decay of vortices is intrinsically dissipative mit1 , it is more realistic to include a small dissipative term ueda ; dis . However, from a mathematical point of view, vortices and vortex lattice can be generated in numerical simulation without including any dissipative term fc . The time evolution of the GP equation is performed as described above using three independent steps. The system passes through dynamically unstable deformed configurations and during time iteration the multiple vortex centers appear for $`\mathrm{\Omega }`$ greater than a critical value.

## III Numerical Results

In the present numerical study on vortices we use an equilibrated condensate with $`𝒩=500`$ trapped in the stationary potential as the initial state. Figure 1 shows the dynamics of the condensate as it starts to rotate suddenly with an angular frequency $`\mathrm{\Omega }=1.5`$. In the simulation we use as input the wave function with phase $`\delta `$ given by Eq. (4). In the beginning the system passes through deformed configurations with unstable oscillating boundaries. The vortices are generated in the peripheral region via the oscillation of the boundary. More and more vortices are generated which move towards the center and which arrange themselves in a regular form in the condensate. For a sufficiently large number of vortices they are arranged as in a lattice and maintain that pattern during rotation. When such a regular equilibrated final pattern is formed the condensate attains the symmetry of the problem $``$ circular in the present case.

In actual experiment the size of the condensate reduces with time due to inelastic atomic collisions mit . However, throughout the present study we renormalize the number of particles during time iteration so that the size is conserved. This should not affect the generation of vortices.

For $`\mathrm{\Omega }=1.5`$ and $`𝒩=500`$ one has a hexagonal lattice structure with seven vortices at large times. In recent theoretical studies fc ; ad ; ueda  some general features of vortex formation have been established which are confirmed in the present study. For a fixed large $`𝒩`$ there is a critical $`\mathrm{\Omega }`$ for the appearance of a certain number of vortices exp ; mit . The critical value of $`\mathrm{\Omega }`$ increases as the number of vortices increases. For a fixed $`\mathrm{\Omega }`$ the number of vortices increases as $`𝒩`$ increases. It is interesting to recall that for $`𝒩=0`$ only states with a single vortex at the center can be formed.

An interesting set of experiments can be performed where on a stable rotating BEC with vortices the strength of atomic interaction or the harmonic oscillator trapping potential is suddenly changed. Once either of these parameters are changed suddenly on a vortex lattice with hexagonal structure, the system enters again a dynamically unstable asymmetric configuration with highly oscillating deformed boundary. Through the oscillating deformed boundary new vortices can enter the condensate or some of the existing vortices can get out. Eventually, a new equilibrated BEC with vortices is formed in the condensate with circular symmetry. However, as the condensate oscillates after the application of the impulsive force it is more difficult to obtain perfect symmetry numerically in the final state. Small imperfections remain near the boundary and in the arrangement of the vortices in the final state.

In the simulation presented in Fig. 2 the nonlinear parameter $`𝒩`$ on the hexagonal vortex lattice of Fig. 1 is suddenly changed to 1200 from 500 at $`t=0`$ and the evolution of the resultant condensate studied. This corresponds to an increase of the strength of interatomic interaction by a factor 2.4. The size of the condensate increases and it accommodates more vortices. The system seems to keep on oscillating with $`(12\pm 1)`$ vortices. This oscillation seems to be a natural consequence of the sudden change of the strength of interaction. Next we consider in Fig. 3 the sudden change in the nonlinear parameter $`𝒩`$ from 500 to 100 at $`t=0`$ on the final configuration of Fig. 1. With the reduction of the strength of interaction by a factor 0.2, the condensate shrinks in size by passing through dynamically unstable and asymmetric configurations and starts to oscillate. The condensate tries to reach a equilibrated stage by getting rid of vortices. However, at large times the system continues to oscillate with ($`6\pm 1`$) vortices.

Next we consider a sudden change in the harmonic oscillator potential. First, the harmonic oscillator potential is increased by a factor of 1.5 on the rotating BEC of Fig. 1 at $`t=0`$. The time evolution of the resultant condensate is shown in Fig. 4. With the sudden increase in the trapping potential the condensate starts to shrink and oscillate. It passes through dynamically unstable asymmetric configurations and reaches an equilibrium stage by getting rid of a couple of vortices. Eventually, a condensate emerges with four vortices. The system executes small oscillations at large times with a fixed number of vortices. In the next simulation the harmonic oscillator potential is reduced by a factor of 0.75 on the configuration of Fig. 1 at $`t=0`$. Because of the sudden decrease in the trapping potential the condensate starts to swell and oscillate. More vortices are generated near the boundary which move inwards. Eventually, one has a larger condensate with 11 vortices at large times.

We study the decay of the vortices of Fig. 1 once the rotation of the BEC is suddenly stopped at $`t=0`$. The dynamics of the decay of vortices has been studied experimentally in Refs. exp ; mit1 . The numerical simulation of the decay of vortices is shown in Fig. 6. The vortices do not decay immediately after stopping the rotation. Rather they stay for a long time after the rotation stops. As time passes vortices get out of the condensate until one with a single vortex results. This BEC with a single vortex survives for a long time. If we consider $`\omega =2\pi \times 219`$ s-1 as in the experiment of Madison et al exp , one unit of present dimensionless time is 0.000726 s. The system gets rid of six vortices in about 600 units of time or in 0.4 s. However, the state with one vortex stays up to $`t3000`$, which gives a lifetime of about 2400 units of time or about 1.6 s. Madison et al exp  measured the half life of a single vortex state and found it to be about 1 s, which is of the same order of magnitude as in the present simulation. However, a quantitative comparison between the two is not to the point because of approximate nature of the present simulation compared to experiment. The present calculation used a quasi two-dimensional model opposed to a full three-dimensional one. Also, the nonlinearity of the present model is not the same as in the experiment. However, in a rotating BEC with many vortices, most of the vortices decay within one second whereas a small number of vortices survive for a long time. This was also experimentally confirmed by Abo-Shaeer et al mit1  who used $`\omega 2\pi \times 85`$ s-1.

Finally, we consider the free expansion of the vortices in Fig. 1 when the harmonic oscillator trap is suddenly removed. The dynamics of free expansion is shown in Fig. 7. It is interesting that the vortex pattern is not destroyed immediately after the removal of the trap. On the other hand, upon expansion the radius of the vortex core is found to increase more rapidly than the radius of the condensate. Due to the very small radius of the vortex core it is difficult to observe the vortices and count them properly in an experiment mit1 . However, because of the increase in the radius of the vortex core after a free expansion the detection and counting of vortices will be easier from a picture of the condensate after expansion mit1 . In Fig. 7 the vortices are more prominently visible at $`t=1.6`$ than at $`t=0.2`$. Using the radial trapping frequency $`\omega =2\pi \times 85`$ s-1 as in the experiment of Abo-Shaeer et al mit1 , $`t=1.6`$ of Fig. 7 corresponds to about 3 ms. In the actual experiment mit1  after 42.5 ms of expansion the vortex cores were magnified 20 times their size in the trap. The present simulation has again led to the right order of magnitude.

## IV Summary

We have investigated the effect of an impulsive force on a rotating BEC with vortices using the mean-field GP equation. We have considered the effect of suddenly increasing or decreasing the strength of atomic interaction and the harmonic oscillator trapping potential. Upon application of such an impulsive force the condensate passes through highly unstable distorted configurations with oscillating boundary. New vortices can enter the condensate via the oscillating boundary. Also the condensate can get rid of some of the vortices in a similar fashion. Eventually, a (more) symmetric condensate with a different number of vortices result. When the atomic interaction is suddenly increased (decreased) the final condensate grows (shrinks) in size and accommodates more (less) vortices. When the trapping potential is suddenly increased (decreased) the final condensate shrinks (grows) in size and accommodates less (more) vortices. Experiments can be performed when similar impulsive force in applied on a rotating BEC with vortices and compared with the mean-field prediction.

We also study the decay of vortices in a rotating BEC when the rotation is suddenly stopped. Most of the vortices decay quickly (in a fraction of a second), but a few vortices survive during a relatively long time before they decay. In a typical experimental situation exp ; mit1  the single vortex at the end survives for a unusually long interval of time (a couple of seconds).

In another simulation of a free expansion of a rotating BEC we find that the vortices are not immediately lost during expansion. Rather, upon expansion the radius of the vortex core increases at a faster rate than that of the condensate. Hence the vortices are more prominently visible in a picture of the condensate during the free expansion of the condensate. Experimentally, it is difficult to observe and count the vortices correctly as the vortex cores have usually small radii mit1 . The present simulation suggests that the free expansion could aid in the observation and counting of vortices in a rotating BEC, specially, when there is a large number of vortices. This was noted in the experiment of Abo-Shaeer et al mit1 . In a typical experimental situation, on free expansion the vortices may be observed during few tens of milliseconds.

###### Acknowledgements.

The work is supported in part by the Conselho Nacional de Desenvolvimento Científico e Tecnológico and Fundação de Amparo à Pesquisa do Estado de São Paulo of Brazil.


text:
# Quark and Nucleon Self-Energy in Dense Matter

## I INTRODUCTION

In recent years there has been a great deal of interest in understanding the properties of quark matter at relatively low temperature and high density \[1-5\]. Since it is difficult to study QCD at finite chemical potential using lattice simulations, the model of choice for such studies has been the Nambu–Jona-Lasinio model \[6-8\]. Of particular interest is the prediction of diquark condensates and color superconductivity at high densities. (Such studies may be relevant to the properties of neutron stars.) Calculations made for dense matter using the NJL model are carried out in Minkowski space. These calculations are limited in that they do not include a model of confinement and are, therefore, unable to provide a comprehensive description of the hadronic phase present at low density and temperature. However, it is generally believed that a good deal of information may be gained by studying quark matter, with a proper study of the hadronic phase deferred until some future time.

As is well known, the NJL model provides a microscopic dynamical description of chiral symmetry breaking with the generation of associated quark vacuum condensates and constituent masses. In the standard versions of the NJL model \[6-8\], the constituent quark mass that is generated in the model is a constant. However, it is known from lattice simulations of QCD that the constituent mass goes over to the current quark mass when the quark momentum, $`p^2`$, is less than about -2 GeV2 . It is our belief that, if we are to use the NJL model to study dense matter, it is desirable to make the model as realistic as possible. To that end, we have introduced a nonlocal version of the NJL model  that is able to reproduce the Euclidean-space behavior of the quark mass seen in lattice simulations of QCD . To carry out that program we have introduced a momentum-dependent $`q\overline{q}`$ interaction in the calculation of the quark self-energy and have separated the regularization of the model from the specification of that interaction. (This procedure requires the introduction of additional parameters into the model.)

Recently, we have seen an attempt to obtain the quark self-energy in Minkowski space by analytic continuation of a Euclidean-space form based upon gluon exchange enhanced at small momentum transfers to simulate confinement . The resulting Minkowski-space values exhibit resonant-like structures and very large values of the constituent mass. Such results do not appear to have any natural physical interpretation.

We note that the quark self-energy may be written as

$`\mathrm{\Sigma }(k^2)=A(k^2)+B(k^2)\text{ / }k`$ (1.1)
in vacuum. In matter there is another four-vector, $`\eta ^\mu `$, that describes the motion of the matter rest frame. It is useful to put $`\eta ^2=1`$ and to note that, if we work in the matter rest frame, we can put $`\eta ^\mu =[1,0,0,0]`$. In matter, we have

$`\mathrm{\Sigma }(k^2,\eta k)=A(k^2,\eta k)+B(k^2,\eta k)\overline{\text{ / }k}+C(k^2,\eta k)\text{ }\eta \text{ }/,`$ (1.2)
where we have found it useful to introduce the four-vector

$`\overline{k^\mu }=k^\mu (k\eta )\eta ^\mu .`$ (1.3)
Note that $`\overline{k^0}=0`$ in the matter rest frame. In that frame, we may write

$`\mathrm{\Sigma }(k^0,\stackrel{}{k})=A(k^0,\stackrel{}{k})B(k^0,\stackrel{}{k})\stackrel{}{\gamma }\stackrel{}{k}+\gamma ^0C(k^0,\stackrel{}{k}).`$ (1.4)
As we will see, A and B satisfy coupled nonlinear equations, while C may be calculated independently. (Note that A and C have the same dimension, while B is dimensionless.)

We first discuss our results for the case in which we neglect the dependence of A, B, and C on $`k^0`$. We anticipate that the dependence on $`k^0`$ will be relatively weak and justify that assumption later in this work. As a further simplification, we will at first neglect B and study the behavior of $`A(\stackrel{}{k},\rho )`$, where $`\rho `$ is the density of quark matter which is taken to contain equal numbers of up and down quarks. In this case, the maximum value of $`A(\stackrel{}{k},\rho )`$ is found at $`\stackrel{}{k}=0`$, leading to a simpler presentation of our results. We then go on to the consideration of the coupled equations for $`A(\stackrel{}{k},\rho )`$ and $`B(\stackrel{}{k},\rho )`$. As usual, we may introduce a density and momentum-dependent mass defined by

$`M(\stackrel{}{k},\rho )={\displaystyle \frac{A(\stackrel{}{k},\rho )}{1B(\stackrel{}{k},\rho )}}.`$ (1.5)
We provide values of $`A(\stackrel{}{k},\rho )`$, $`M(\stackrel{}{k},\rho )`$ and $`C(\stackrel{}{k},\rho )`$ in the following sections.

We note that the inclusion of $`C(\stackrel{}{k},\rho )`$ precludes the passage to Euclidean space. The analogous problem arises when one introduces a finite chemical potential. As we will see, $`C(\stackrel{}{k},\rho )`$ is quite large at finite density and can only be neglected at very small densities. In vacuum, Lorentz invariance leads to a relation between $`C(k^2)`$ and $`B(k^2)`$. However, our formalism does not maintain Lorentz invariance. (For example, our regulator depends only upon $`|\stackrel{}{k}|`$.) Therefore, in the following we will only calculate the contribution to $`C(\stackrel{}{k},\rho )`$ from the matter, which takes the form of two Fermi seas of up and down positive-energy quarks with Fermi momentum $`k_F`$.

The Lagrangian of our model is

$`=`$ $`\overline{q}(i\text{ / }m^0)q+{\displaystyle \frac{G_S}{2}}{\displaystyle \underset{i=0}{\overset{8}{}}}[(\overline{q}\lambda ^iq)^2+(\overline{q}i\gamma _5\lambda ^iq)^2]`$ (1.6)
$`{\displaystyle \frac{G_V}{2}}{\displaystyle \underset{i=0}{\overset{8}{}}}[(\overline{q}\lambda ^i\gamma _\mu q)^2+(\overline{q}\lambda ^i\gamma _5\gamma _\mu q)^2]`$
$`+{\displaystyle \frac{G_D}{2}}\{det[\overline{q}(1+\gamma _5)q]+det[\overline{q}(1\gamma _5)q]\}`$
$`+_{\text{conf}}.`$
This Lagrangian was used in Ref.  with $`G_V=0`$. There, both the condensate values and self-energies were given for the up (down) and strange quarks. For the present work we neglect the ’t Hooft interaction $`(G_D=0)`$ and our model of confinement. We also drop the strange quark from consideration, so that the quark current mass matrix is $`m^0=\text{ diag}(m_u^0,m_d^0)`$ with $`m_u^0=m_d^0`$. Since the ’t Hooft interaction contributes to the self-energy , its neglect leads us to use a larger value of $`G_S`$ than that we would use if the ’t Hooft interaction were included in our analysis.

The organization of our work is as follows. In Section II we present the coupled nonlinear equations that determine $`A(\stackrel{}{k},\rho )`$ and $`B(\stackrel{}{k},\rho )`$ and also provide an expression for $`C(\stackrel{}{k},\rho )`$, and the density-dependent condensate $`\overline{u}u_\rho =\overline{d}d_\rho `$. We present values of $`A(\stackrel{}{k},\rho )`$ and the condensate, when we neglect $`B(\stackrel{}{k},\rho )`$. In Section III we consider finite values of $`B(\stackrel{}{k},\rho )`$ and provide results for $`A(\stackrel{}{k},\rho )`$, $`M(\stackrel{}{k},\rho )`$, $`C(\stackrel{}{k},\rho )`$ and the condensate. In Section IV, we consider the dependence of A, B and C on $`k^0`$ and provide some results for $`A(k_0,\stackrel{}{k})`$ in vacuum. In Section V we use the quark self-energy calculated here to obtain the nucleon self-energy in nuclear matter in a simple model. Finally, Section VI contains some further discussion and conclusions.

## II the quark self-energy

In our earlier work we obtained the quark self-energy from the solution of the equation depicted in Fig. 1 . There, the open circle is a momentum-dependent quark interaction, obtained by the replacement

$`G_Sf(kk^{})G_Sf(kk^{}),`$ (2.1)
where $`k`$ and $`k^{}`$ are the quark momenta entering (or leaving) the interaction. We have used

$`f(kk^{})=\text{exp}[(kk^{})^{2n}/2\beta ]`$ (2.2)
with $`n=4`$ and $`\beta =20`$ GeV8 in Ref. . The corresponding nonlocal Lagrangian is given in Ref. . On the right-hand side of Fig.1, the ’t Hooft interaction (third term) and the confinement interaction (fourth term) are neglected for the purposes of this work. In the second term we have contributions from the negative-energy states in vacuum as well as the positive-energy states at finite density, specified by the quark Fermi momentum, $`k_F`$, of the up and down quark Fermi seas.

In general, the quark propagator is

$`iS(k,k\eta )={\displaystyle \frac{i}{\text{ / }k\mathrm{\Sigma }(k^2,k\eta )+iϵ}},`$ (2.3)
which we will write as

$`iS(k^0,\stackrel{}{k})={\displaystyle \frac{i}{(k^0C(k^0,\stackrel{}{k}))\gamma ^0(1B(k^0,\stackrel{}{k}))\stackrel{}{\gamma }\stackrel{}{k}A(k^0,\stackrel{}{k})+iϵ}},`$ (2.4)
in the matter rest frame. In a first approximation we write

$`iS(k^0,\stackrel{}{k})={\displaystyle \frac{i}{(k^0C(\stackrel{}{k}))\gamma ^0(1B(\stackrel{}{k}))\stackrel{}{\gamma }\stackrel{}{k}A(\stackrel{}{k})+iϵ}},`$ (2.5)
with $`A(\stackrel{}{k})`$, $`B(\stackrel{}{k})`$ and $`C(\stackrel{}{k})`$ density-dependent, in general. (On occasion we will write $`A(\stackrel{}{k},\rho )`$, etc.) We see that the presence of $`C(\stackrel{}{k})`$ precludes the passage to Euclidean space that was made in Ref.. That is analogous to the problem created by the introduction of a chemical potential in the formalism. Note that

$`(k^0C(\stackrel{}{k}))^2(1B(\stackrel{}{k}))^2\stackrel{}{k}^2A^2(\stackrel{}{k})+iϵ`$
$`=[k^0E^+(\stackrel{}{k})+iϵ][k^0E^{}(\stackrel{}{k})iϵ]`$ (2.6)
with

$`E^+(\stackrel{}{k})=C(\stackrel{}{k})+\sqrt{\stackrel{}{k}^2(1B(\stackrel{}{k}))^2+A^2(\stackrel{}{k}^2)}`$ (2.7)
and

$`E^{}(\stackrel{}{k})=C(\stackrel{}{k})\sqrt{\stackrel{}{k}^2(1B(\stackrel{}{k}))^2+A^2(\stackrel{}{k}^2)}.`$ (2.8)
Here $`E^+(\stackrel{}{k})`$ may be interpreted as the (on-mass-shell) energy of the positive-energy states, while $`E^{}(\stackrel{}{k})`$ refers to the negative-energy states.

In order to simplify the notation somewhat, we write

$`iS(k^0,\stackrel{}{k})={\displaystyle \frac{i}{\text{ / }\mathrm{\Pi }(k^0,\stackrel{}{k})A(\stackrel{}{k})+iϵ}},`$ (2.9)
where we have defined a four-component quantity

$`\mathrm{\Pi }^\mu (k^0,\stackrel{}{k})=[k^0C(\stackrel{}{k}),(1B(\stackrel{}{k}))\stackrel{}{k}].`$ (2.10)
We also introduce a scalar quantity

$`\rho _S(\stackrel{}{k})=iN_c(1){\displaystyle \frac{d^4k^{}}{(2\pi )^4}\frac{4A(\stackrel{}{k}^{})f^2(kk^{})}{\mathrm{\Pi }^2(k^{\mathrm{\hspace{0.17em}0}},\stackrel{}{k}^{})A^2(\stackrel{}{k}^{})+iϵ}},`$ (2.11)
where the minus sign is due to the closed Fermion loop in Fig. 1 and the factor of 4 comes from forming the trace associated with the closed loop. We see that $`\rho _S(\stackrel{}{k})`$ does not depend upon $`k^0`$, since we are making an on-shell approximation in the calculation of $`f(kk^{})`$. In vacuum we have

$`iA(\stackrel{}{k})=im^0+(2G_Si)\rho _S(\stackrel{}{k}),`$ (2.12)
or

$`A(\stackrel{}{k})=m^02G_S\rho _S(\stackrel{}{k}),`$ (2.13)
where $`\rho _S(\stackrel{}{k})`$ is real and negative. We remark when $`f(kk^{})=1`$, $`A(\stackrel{}{k})m`$ and $`\rho _S(\stackrel{}{k})\overline{u}u`$, so that we regain the usual result \[6-8\]

$`m_u=m_u^02G_S\overline{u}u,`$ (2.14)
$`m_d=m_d^02G_S\overline{d}d.`$ (2.15)
(Note that our $`G_S`$ is one-half of the $`G_S`$ defined in Ref. .)

The evaluation of $`\rho _S(\stackrel{}{k})`$ proceeds by closing to contour in the complex $`k^0`$ plane. In the vacuum we find

$`\rho _S(\stackrel{}{k})=2N_c{\displaystyle \frac{d^3k^{}}{(2\pi )^3}\frac{f^2(kk^{})A(\stackrel{}{k}^{})}{\sqrt{\stackrel{}{k}^2(1B(\stackrel{}{k}^{}))^2+A^2(\stackrel{}{k}^{})}}}.`$ (2.16)
We define

$`E(\stackrel{}{k})=\sqrt{\stackrel{}{k}^2(1B(\stackrel{}{k}))^2+A^2(\stackrel{}{k})},`$ (2.17)
and

$`(kk^{})^2=(E(\stackrel{}{k})E(\stackrel{}{k}^{}))^2(\stackrel{}{k}\stackrel{}{k}^{})^2,`$ (2.18)
so that $`f(kk^{})`$ depends only upon $`|\stackrel{}{k}|`$, $`|\stackrel{}{k}^{}|`$ and the angle between $`\stackrel{}{k}`$ and $`\stackrel{}{k}^{}`$.

Integrals such as that in Eq. (2.16) require regularization. For our calculations we insert a factor $`\text{exp}[\stackrel{}{k}^2/\alpha ^2]`$ with $`\alpha =0.60`$ GeV. (We have used the same Gaussian regulator in our calculations of meson spectra, where we have used $`\alpha =0.605`$ GeV \[12-14\]. However, those calculations included a model of confinement, so that we can not directly take over the parameters $`G_S`$, $`G_V`$ and $`G_D`$ used in those works.) We remark that an expression for the density-dependent condensate $`\overline{u}u_\rho `$ may be obtained by using Eq. (2.16) with $`f(kk^{})=1`$.

In the presence of matter it is useful to separate the propagator into two parts, one of which will give rise to the explicitly density-dependent terms. In this regard, it is useful to generalize Eq. (5.8) of Ref. . We define

$`\mathrm{\Lambda }^{(+)}(\stackrel{}{k})={\displaystyle \frac{E(\stackrel{}{k})\gamma ^0\stackrel{}{\gamma }\stackrel{}{k}(1B(\stackrel{}{k}))+A(\stackrel{}{k})}{2A(\stackrel{}{k})}},`$ (2.19)
$`\mathrm{\Lambda }^{()}(\stackrel{}{k})={\displaystyle \frac{E(\stackrel{}{k})\gamma ^0\stackrel{}{\gamma }\stackrel{}{k}(1B(\stackrel{}{k}))+A(\stackrel{}{k})}{2A(\stackrel{}{k})}},`$ (2.20)
where $`E(\stackrel{}{k})=[\stackrel{}{k}^2(1B(\stackrel{}{k}))^2+A^2(\stackrel{}{k})]^{1/2}`$. Then

$`S(k)=`$ $`{\displaystyle \frac{A(\stackrel{}{k})}{E(\stackrel{}{k})}}\left[{\displaystyle \frac{\mathrm{\Lambda }^{(+)}(\stackrel{}{k})}{k^0E^+(\stackrel{}{k})}}{\displaystyle \frac{\mathrm{\Lambda }^{()}(\stackrel{}{k})}{k^0E^{}(\stackrel{}{k})}}\right]`$
$`+`$ $`2\pi i{\displaystyle \frac{A(\stackrel{}{k})}{E(\stackrel{}{k})}}\mathrm{\Lambda }^{(+)}(\stackrel{}{k})\theta (k_F|\stackrel{}{k}|)\delta (k^0E^+(\stackrel{}{k}))`$
In the limit that $`A(\stackrel{}{k})m^{}`$ and $`B(\stackrel{}{k})=C(\stackrel{}{k})=0`$, we have, with

$`E^{}(\stackrel{}{k})=\sqrt{\stackrel{}{k}^2+m^2},`$ (2.22)
$`S(k)`$ $`=`$ $`{\displaystyle \frac{m^{}}{E^{}(\stackrel{}{k})}}\left[{\displaystyle \frac{\mathrm{\Lambda }^{(+)}(\stackrel{}{k})}{k^0E^{}(\stackrel{}{k})}}{\displaystyle \frac{\mathrm{\Lambda }^{()}(\stackrel{}{k})}{k^0+E^{}(\stackrel{}{k})}}\right]`$
$`+`$ $`2\pi i{\displaystyle \frac{m^{}}{E^{}(\stackrel{}{k})}}\mathrm{\Lambda }^{(+)}(\stackrel{}{k})\theta (k_F|\stackrel{}{k}|)\delta (k^0E^{}(\stackrel{}{k}))`$
$`=`$ $`{\displaystyle \frac{\text{ / }k+m^{}}{k^2m^2}}+{\displaystyle \frac{i\pi }{E^{}(\stackrel{}{k})}}(\text{ / }k+m^{})\theta (k_F|\stackrel{}{k}|)\delta (k^0E^{}(\stackrel{}{k})),`$ (2.24)
which agrees with Eq. (5.8) of Ref. .

In the presence of matter, we have

$`A(\stackrel{}{k},\rho )=m^02G_S[\rho _S^{\text{vac}}(\stackrel{}{k})\rho _S^{\text{mat}}(\stackrel{}{k})],`$ (2.25)
where $`\rho _S^{\text{mat}}(\stackrel{}{k})`$ is calculated in the same manner as $`\rho _S^{\text{vac}}(\stackrel{}{k})`$, except that the upper limit of the integral over $`|\stackrel{}{k}^{}|`$ is $`k_F`$. Equation (2.21) is a generalization of a corresponding equation that may be found in Klevansky’s review. (See Eqs. (5.18) of Ref. .)

If we neglect $`B(\stackrel{}{k})`$, Eqs. (2.11) and (2.16) provide a nonlinear equation for $`A(\stackrel{}{k})`$ which may be solved by iteration. The results of such a calculation are reported in Table I where, for nuclear matter, we put $`k_F=0.268`$ GeV. In Figs. 2 and 3 we show values of the condensate $`\overline{u}u`$ and $`A(0,\rho )`$ as a function of the density. These results may be usefully discussed in terms of the relation

$`\overline{u}u_\rho =\overline{u}u_0\left(1{\displaystyle \frac{\sigma _N\rho _N}{f_\pi ^2m_\pi ^2}}+\mathrm{}\right).`$ (2.26)
Here $`\sigma _N`$ is the pion-nucleon sigma term and $`\rho _N`$ is the density of nucleons. If we put $`\sigma _N=0.050`$ GeV, $`\rho _N=(0.109\text{GeV})^3`$, $`f_\pi =0.0942`$ GeV and $`m_\pi =0.138`$ GeV, we find a 38% reduction of the condensate at nuclear matter density, which agrees with our results given in Table I and Fig. 2. It is of interest to note that the linear dependence on the density implied by Eq. (2.26) appears to be valid up to about twice nuclear matter density. However, one may be concerned that, since we study quark matter rather than nuclear matter, Eq. (2.26) may not be appropriate. Consider, however, the relation

$`\overline{u}u_\rho =\overline{u}u_0\left(1{\displaystyle \frac{\sigma _q\rho _q}{f_\pi ^2m_\pi ^2}}+\mathrm{}\right),`$ (2.27)
where $`\sigma _q`$ is the quark sigma term and $`\rho _q`$ is the number density of the quarks, which we may put equal to $`3\rho _N`$. Thus, if $`3\sigma _q=\sigma _N`$, we may use Eq. (2.26). The fact that $`\sigma _q15`$ MeV has been discussed by Vogl and Weise . We have also discussed this matter in great detail in Ref. , where we calculated similar values of $`\sigma _q`$ using the standard version of the NJL model. We conclude that the use of Eq. (2.26), or Eq. (2.27), with an appropriate value of $`\sigma _q`$, is satisfactory. For example, if $`\sigma _q=\sigma _N/3`$, as suggested in Ref. , the two relations imply the same density dependence of the condensate.

In Fig. 3 we show $`A(0,\rho )`$ which is the density-dependent mass parameter of the theory when $`B(\stackrel{}{k},\rho )=0`$. We see that $`A(0,\rho )`$ follows the trend seen in Fig.2 for the density dependence of the condensate.

## III lorentz-vector terms of the quark self-energy

We now write $`\mathrm{\Sigma }(\stackrel{}{k})=\mathrm{\Sigma }_S(\stackrel{}{k})+\mathrm{\Sigma }_V(\stackrel{}{k})`$ where

$`\mathrm{\Sigma }_V(\stackrel{}{k})=\stackrel{}{\gamma }\stackrel{}{k}B(\stackrel{}{k})+\gamma ^0C(\stackrel{}{k}).`$ (3.1)
For the calculation of $`C(\stackrel{}{k})`$ we obtain the contribution from the last term in Eq. (2.21), with the result that

$`C(\stackrel{}{k})=2G_V\rho _2^V(\stackrel{}{k})`$ (3.2)
with

$`\rho _2^V(\stackrel{}{k})=2N_c{\displaystyle ^{k_F}}{\displaystyle \frac{d^3k^{}}{(2\pi )^3}}f^{\mathrm{\hspace{0.17em}2}}(kk^{}).`$ (3.3)
An expression for $`B(\stackrel{}{k})`$ may be found from the relation

$`i[\stackrel{}{\gamma }\stackrel{}{k}B(\stackrel{}{k})]=(2G_Vi)(1)i{\displaystyle \frac{d^4k^{}}{(2\pi )^4}S(k^{})f^{\mathrm{\hspace{0.17em}2}}(kk^{})}`$ (3.4)
if we only keep the term proportional to $`\stackrel{}{\gamma }\stackrel{}{k}^{}`$ in the expression for the quark propagator. In vacuum we may compare corresponding terms in Eq. (3.4)

$`\stackrel{}{\gamma }\stackrel{}{k}B(\stackrel{}{k})=2G_V{\displaystyle \frac{d^4k^{}}{(2\pi )^4}\frac{\stackrel{}{\gamma }\stackrel{}{k}^{}[1B(\stackrel{}{k}^{})]f^{\mathrm{\hspace{0.17em}2}}(kk^{})}{[k_0^{}E^+(\stackrel{}{k}^{})+iϵ][k_0^{}E^{}(\stackrel{}{k}^{})iϵ]}}.`$ (3.5)
Thus,

$`B^{\text{vac}}(\stackrel{}{k})=2G_V\rho _1^{\text{vac}}(\stackrel{}{k}),`$ (3.6)
with

$`|\stackrel{}{k}|\rho _1^{\text{vac}}(\stackrel{}{k})=2N_c{\displaystyle \frac{d^3k^{}}{(2\pi )^3}\frac{|\stackrel{}{k}^{}|(\widehat{k}\widehat{k}^{})f^{\mathrm{\hspace{0.17em}2}}(kk^{})[1B(\stackrel{}{k}^{})]}{\sqrt{\stackrel{}{k}^2[1B(\stackrel{}{k}^{})]^2+A^2(\stackrel{}{k}^{})}}}.`$ (3.7)
Here, $`\widehat{k}`$ and $`\widehat{k}^{}`$ are unit vectors. As in the calculation of $`A(\stackrel{}{k})`$, Eq. (3.6) is generalized to read

$`B(\stackrel{}{k})=2G_V[\rho _1^{\text{vac}}(\stackrel{}{k})\rho _1^{\text{mat}}(\stackrel{}{k})],`$ (3.8)
where $`\rho _1^{\text{mat}}(\stackrel{}{k})`$ is calculated using Eq. (3.7) with an upper limit on $`|\stackrel{}{k}^{}|`$ of $`k_F`$.

In Table II we present results of our calculation of the condensate, $`A(0,\rho )`$, $`A(0,\rho )/[1B(0,\rho )]`$, $`C(0,\rho )`$, $`B(0,\rho )`$ and $`A(0,\rho )A(0,0)`$. We also define

$`U_S(\stackrel{}{k},\rho )=A(\stackrel{}{k},\rho )A(\stackrel{}{k},0),`$ (3.9)
where $`U_S(\stackrel{}{k},\rho )`$ is the density-dependent modification of $`A(\stackrel{}{k},0)`$ in matter.

It may be seen from the values given in Table II that there is a thirty percent reduction of the condensate at the density of nuclear matter, while the value of $`A(0)`$ is reduced by thirty-nine percent. We would obtain a thirty percent reduction of the condensate if $`\sigma _N=39`$ MeV.

In Figs. 4 and 5 we exhibit values of $`A(\stackrel{}{k},\rho )`$ and $`A(\stackrel{}{k},\rho )/[1B(\stackrel{}{k},\rho )]`$ for various densities and in Fig. 6 we present values of $`C(\stackrel{}{k},\rho )`$. Figure 7 shows the values of $`U_S(\stackrel{}{k},\rho _{NM})`$ and $`C(\stackrel{}{k},\rho _{NM})`$.

## IV off-mass-shell effects

In our calculations we have neglected the $`k_0`$ dependence of A, B and C. In this Section we provide some justification of that approximation. We can combine Eqs. (2.9) and (2.11) to yield a nonlinear equation for $`A(k^0,\stackrel{}{k})`$ in vacuum

$`A(k^0,\stackrel{}{k})=(2G_Si)N_c(1){\displaystyle \frac{d^4k^{}}{(2\pi )^4}\frac{4A(k^0,\stackrel{}{k}^{})f^{\mathrm{\hspace{0.17em}2}}(kk^{})}{\text{ / }k^{}A(k^0,\stackrel{}{k}^{})+iϵ}}.`$ (4.1)
Here we have neglected $`m^0`$ and $`B(k^0,\stackrel{}{k})`$. We see that $`k^0`$ dependence arises from $`f^{\mathrm{\hspace{0.17em}2}}(kk^{})`$ which is given by

$`f^{\mathrm{\hspace{0.17em}2}}(kk^{})=\text{exp}[(kk^{})^{\mathrm{\hspace{0.17em}2}n}/\beta ],`$ (4.2)
with

$`(kk^{})^2=(k^0k^0)^2(\stackrel{}{k}\stackrel{}{k}^{})^2,`$ (4.3)
which differs from the on-mass-shell version given in Eq. (2.18).

The solution of Eq. (4.1) for $`A(k^0,\stackrel{}{k})`$ is shown in Fig. 8 for various values of $`k^0`$. It is seen that the dependence on $`k^0`$ is weak as was assumed in this work.

## V the nucleon self-energy in matter

If one uses the Dirac equation to describe the interaction of a nucleon with a nucleus, or with nuclear matter, it is found that a strong scalar attraction is needed as well as a strong vector repulsion . The scalar field is of the order of -400 MeV and the vector field is about 300 MeV. It is of interest to see if the nucleon self-energy, $`\mathrm{\Sigma }_N=V_S+\gamma ^0V_V`$, can be calculated in terms of the quark self-energy obtained in this work. To carry out this program we use a simple model of the nucleon in which a quark is coupled to a scalar diquark. The full complexity of the wave function, including vector diquarks and various relativistic effects, is discussed in Ref. .

We can calculate the nucleon self-energy in nuclear matter using a triangle diagram in which one of the lower two vertices of the triangle represents a vertex function for a zero-momentum nucleon to emit a quark of momentum $`\stackrel{}{k}`$ leaving a spectator (on-mass-shell) diquark of momentum $`\stackrel{}{k}`$. The other lower vertex represents the inverse process. At the upper vertex we insert the quark self-energy, $`\mathrm{\Sigma }(\stackrel{}{k})=U_S(\stackrel{}{k})+\gamma ^0C(\stackrel{}{k})`$ calculated in this work and integrate over $`\stackrel{}{k}`$. We make use of Fig. 4 of Ref.  and parametrize the product of the vertex function and the quark Greens function by the quark-diquark wave function

$`\psi (\stackrel{}{k})={\displaystyle \frac{1}{\sqrt{N}}}e^{\stackrel{}{k}^2/\lambda ^2}u(\stackrel{}{k},s),`$ (5.1)
where $`u(\stackrel{}{k},s)`$ is a spinor for a quark of momentum $`\stackrel{}{k}`$ and spin projection $`s`$ .

The normalization for a single quark is obtained from the relation

$`{\displaystyle \frac{1}{N}}{\displaystyle \frac{d^3k}{(2\pi )^3}e^{2\stackrel{}{k}^2/\lambda ^2}\left(1+\frac{\stackrel{}{k}^2}{[E_q(\stackrel{}{k})+m_q]^2}\right)}=1,`$ (5.2)
where we put $`\lambda =0.18`$ GeV to correspond to the results of Ref. . We may then relate the nucleon self-energy to the quark self-energy. For a nucleon of momentum $`\stackrel{}{P}=0`$, we have, with $`E_q(\stackrel{}{k})=[\stackrel{}{k}^2+m_q^2]^{1/2}`$, and $`m_q=0.364`$ GeV,

$`V_V={\displaystyle \frac{3}{N}}{\displaystyle \frac{d^3k}{(2\pi )^3}C(\stackrel{}{k})e^{2\stackrel{}{k}^2/\lambda ^2}\left(1+\frac{\stackrel{}{k}^2}{[E_q(\stackrel{}{k})+m_q]^2}\right)},`$ (5.3)
and

$`V_S={\displaystyle \frac{3}{N}}{\displaystyle \frac{d^3k}{(2\pi )^3}U_S(\stackrel{}{k})e^{2\stackrel{}{k}^2/\lambda ^2}\left(1\frac{\stackrel{}{k}^2}{[E_q(\stackrel{}{k})+m_q]^2}\right)}.`$ (5.4)
The difference of sign in the brackets appearing in Eqs. (5.3) and (5.4) is due to the different behavior of the Dirac matrices, $`\mathrm{𝟏}`$ and $`\gamma ^0`$, at the upper vertex of the triangle.

Since the momentum content of the quark-diquark wave function is small , we expect that $`V_S3U_S(0)`$ and $`V_V3C(0)`$, so that $`V_V296`$ MeV and $`V_S426`$ MeV. A more careful evaluation of the integrals in Eqs. (5.3) and (5.4) yields $`V_V=295`$ MeV and $`V_S=392`$ MeV, which is in general accord with the values given in Refs.  and .

## VI discussion

The behavior of matter at high density and low temperature has received a good deal of attention in the last few years \[1-5\]. A large part of the work in this area has made use of the NJL model. In our work we have attempted to modify the NJL model so that its predictions are in greater accord with QCD. As a first step, in Ref.  we introduced a momentum-dependent $`q\overline{q}`$ interaction which allowed us to reproduce the Euclidean-space behavior for the constituent quark mass obtained in lattice simulation of QCD . In the present study we have considered the important vector interactions of an extended NJL model. Of particular interest is the behavior of the quark condensate at finite density. We find that the linear behavior in the density exhibited in Eq. (2.26) holds in our model up to about twice nuclear matter density. Also, we note that the use of a nonzero current quark mass is important for that result, since in the absence of an explicit chiral symmetry breaking term, the model exhibits a first-order phase transition at about 1.25 times nuclear matter density.

Another point of interest are the results shown in Fig. 7. The quark self-energy is similar to that found in relativistic nuclear physics, with strong scalar attraction and strong vector repulsion. The simple calculation reported in Section V suggests that the quark self-energy, when multiplied by 3, provides a satisfactory estimate of the nucleon self-energy which is in accord with results of relativistic nuclear physics .

There is a body of work based upon the solutions of the Schwinger-Dyson and Bethe-Salpeter equations with a phenomenological form for the gluon propagator \[20-23\]. This body of work is reviewed in Ref. . In contrast to the results of our work, it is found that the quark condensate increases with increasing chemical potential. The authors of Ref.  argue that the baryon density is zero up to the critical chemical potential, $`\mu _c`$, for deconfinement. They state that “This result is an expected consequence of confinement which entails that each additional quark must be locally paired with an antiquark thereby increasing the density of condensate pairs as $`\mu `$ is increased. For this reason, as long as $`\mu <\mu _c`$, there is no excess of particles over antiparticles in the vacuum and hence the baryon number density remains zero.” We note, however, that the baryon density usually considered in such calculations is due to the presence of nucleons, which have a finite baryon density. As noted earlier, in our work the baryon density is due to the presence of Fermi seas of up and down quarks, which serve to provide a model of the baryon density that would arise due to the presence of nucleons. We suggest that the conclusions presented in Ref.  do not refer to the situation of physical interest in which one studies the value of the quark condensate in nuclei or in nuclear matter.

References