text:
# The Black Hole Masses and Host Galaxies of BL Lac Objects

## 1 Introduction

One of the goals of AGN research is to develop a unified framework in which the diversity of the AGN family might be understood in terms of variations in a few fundamental parameters such as the black hole mass, the ratio of the accretion rate to the Eddington accretion rate, and the orientation relative to our line of sight (e.g., Laor, 2000; Boroson, 2002). Determining the masses of the central black holes in different classes of AGNs is an important step toward this goal. An key advance in this area was the recent discovery of the $`M_{\mathrm{BH}}\sigma `$ relation, a tight correlation between black hole mass and stellar velocity dispersion in the host galaxy bulge (Ferrarese & Merritt, 2000; Gebhardt et al., 2000a). The $`M_{\mathrm{BH}}\sigma `$ relation has provided an important consistency check for black hole masses determined from reverberation mapping of Seyfert nuclei (Nelson, 2000; Gebhardt et al., 2000b; Ferrarese et al., 2001). The tightness of the correlation over a wide range of host galaxy types, both elliptical and spiral, makes it tremendously useful as a means to determine the black hole masses in AGNs. While stellar-dynamical or gas-dynamical black hole mass measurements are particularly difficult for AGNs, stellar velocity dispersions can be measured relatively easily in some classes of active galaxies. Furthermore, black hole mass estimates based on velocity dispersions are expected to be much more accurate than estimates derived from the loose correlation between $`M_{\mathrm{BH}}`$ and bulge luminosity (Magorrian et al., 1998; Kormendy & Gebhart, 2001).

We recently began a program to measure the stellar velocity dispersions in the host galaxies of low-redshift BL Lac objects in order to apply the $`M_{\mathrm{BH}}\sigma `$ relation and determine their black hole masses. The first measurement of a velocity dispersion in a BL Lac object, for the TeV $`\gamma `$-ray source Markarian 501, indicated a likely black hole mass of $`(0.93.4)\times 10^9`$ $`M_{}`$ (Barth, Ho, & Sargent, 2002a, Paper I). In this paper, we report on measurements for a sample of 11 objects and discuss the results in the context of unification models for radio-loud AGNs. We also compare our results with measurements of stellar velocity dispersions for BL Lac objects reported recently by Falomo, Kotilainen, & Treves (2002, hereinafter FKT).

## 2 Observations and Reductions

The spectra were obtained at the Hale 5m telescope at Palomar Observatory, and at the 2.5m du Pont and 6.5m Baade (Magellan 1) telescopes at Las Campanas Observatory. We used the Double Spectrograph at Palomar (Oke & Gunn, 1982), the ModSpec Spectrograph at the du Pont telescope, and the Boller & Chivens Spectrograph at Magellan. Details of the instrumental setups are given in Table 1; Table 2 describes the exposures for each object. Conditions were mediocre during the Palomar nights, with typical seeing of 1$`\stackrel{}{\mathrm{.}}`$5 to 2$`\stackrel{}{\mathrm{.}}`$0 and thick clouds at times, while the Las Campanas nights were clear. The gratings and grating angles were chosen to cover the most useful spectral regions for measuring velocity dispersions: the Mg I*b* region in the blue, and the Ca II near-infrared triplet in the red. The values for the seeing listed in Table 2 were measured from observations of template and standard stars, but these were mostly observed during twilight and the seeing during the BL Lac observations may have been somewhat different. For each galaxy, the slit was oriented close to the parallactic angle at the midpoint of the observation, and the total exposure time was divided into individual exposures of typically 1800 s.

Bias subtraction, flat-fielding, spectral extractions, and wavelength calibration were done using standard routines in IRAF, and the flux calibration and subsequent measurements were performed in IDL. Spectra were extracted using the optimal algorithm of Horne (1986) as implemented in IRAF; the extraction widths are listed in Table 2. The strongest telluric absorption bands were removed by division by the normalized spectrum of a nearly featureless star following methods similar to those of Wade & Horne (1988). After wavelength calibration, a final shift to the wavelength scale was applied to each exposure based on the wavelengths of night sky emission lines. This step ensured that the stellar absorption features were not artificially broadened by summing exposures with slightly offset wavelength scales.

During each observing run, several stars were observed for use as velocity templates. The template stars were chosen to be primarily giants of spectral type between G8 and K5, but we also observed some subgiants and dwarfs and also a few F and early-G stars. In the following analysis, we consider only the G and K giant templates, since they give the best match to the dominant stellar populations in the blazar host galaxies.

In addition to the sample listed in Table 1, we also observed the object BL Lac at Palomar, but were unable to measure its velocity dispersion. The stellar absorption lines were only barely detected in a total exposure of 7200 s due to overwhelming AGN contamination. Also, BL Lac has a Galactic extinction of $`A_B=1.42`$ mag (Schlegel, Finkbeiner, & Davis, 1998), and the dominant features over much of its optical spectrum are the Galactic diffuse interstellar bands (Vermeulen et al., 1995).

## 3 Measurement of Velocity Dispersion

Measurement of the velocity dispersion of a BL Lac host galaxy is straightforward in principle, but there are some complications arising from the dilution of the host galaxy starlight by nonthermal emission. For the observations presented here, the detected flux within the spectroscopic aperture is $`5080\%`$ nonthermal continuum emission. As a result, a BL Lac object spectrum yields a more uncertain velocity dispersion than a normal galaxy spectrum of the same signal-to-noise (S/N) ratio. Furthermore, Galactic or telluric absorption features, if not corrected for, will have a relatively larger impact on the velocity dispersion measurements for the BL Lac objects.

We measured the velocity dispersions of the BL Lac hosts by fitting broadened stellar spectra to the object spectra, following a method similar to those described by van der Marel (1994), Rix et al. (1995), and Kelson et al. (2000). The object and template spectra are first rebinned to a wavelength scale that is linear in $`\mathrm{log}\lambda `$, and the object spectrum is transformed to zero redshift. The template spectrum is broadened by convolution with a Gaussian, and adjusted to match the galaxy spectrum by adding a linear featureless continuum and multiplying by a low-order polynomial. The fit includes six free parameters: the systemic velocity of the galaxy, the width of the Gaussian broadening function, the normalization and slope of the featureless continuum, and the three coefficients of a quadratic polynomial. Multiplication by the polynomial allows the fit to account for any reddening or other differences in overall spectral shape between the galaxy and the template, without introducing any high-frequency structure that might affect the dispersion measurement. The fit is optimized using a downhill simplex algorithm (Press et al., 1988). For each template, the fit is repeated several times with a range of initial parameter values to ensure that the global minimum of $`\chi ^2`$ is found. The $`1\sigma `$ uncertainty for each template is determined by finding how far the velocity dispersion must be varied from its best-fit value (with all other parameters allowed to float freely) in order to increase $`\chi ^2`$ by unity relative to its minimum value.

The wavelength range for the fits was selected to include the regions most sensitive to the stellar velocity dispersion. For each object, we attempted to use the largest wavelength range over which a good fit could be achieved, excluding emission lines or regions that consistently gave poor fits with a range of template stars. In general, velocity dispersions measured from the Ca triplet lines are preferred, if they can be observed, because the Ca triplet is less sensitive to template mismatch than the blue spectral region (Dressler, 1984). Also, the Ca triplet lines are intrinsically strong and well isolated from neighboring lines, and can be detected even in highly diluted spectra. For the lowest-redshift objects (Mrk 421, Mrk 501, and 2201+044), we used the Ca triplet lines to measure $`\sigma `$.

At recession velocities beyond $`9000`$ km s-1, the Ca II $`\lambda 8662`$ line is redshifted into deep telluric water vapor absorption bands that begin at about 8925 Å. Reliable velocity dispersions can still be measured using the Ca lines at 8498 and 8542 Å out to $`cz12000`$ km s-1. For the objects having larger redshifts, we used the region around Mg I*b* to measure the velocity dispersions. Best results were obtained for the region $`52055450`$ Å, which contains reasonably strong features due to Ca and Fe. We also included a wavelength range just blueward of Mg I*b* in the fit, typically $`50305150`$ Å, for four objects having sufficient S/N in this region. The wavelength region was adjusted for some objects to exclude Galactic absorption features or bad pixels, or to include a wider wavelength range when a good fit could be obtained.

It proved impossible to fit both the $`52055450`$ Å region (which contains the Fe 5270 and Fe 5335 blends) and the Mg I*b* line simultaneously. This is a consequence of the correlation between the \[Mg/Fe\] ratio and velocity dispersion in elliptical galaxies (Worthey, Faber, & Gonzalez, 1992; Trager et al., 1998; Kuntschner et al., 2001). Since the Mg lines increase in strength more steeply than do the Fe lines as a function of $`\sigma `$, it is difficult to find local template stars that will provide an exact match to the \[Mg/Fe\] ratio of any individual giant elliptical galaxy. This source of template mismatch is a particular problem for direct fitting measurements, since the calculation of $`\chi ^2`$ is sensitive to both the width and the depth of absorption lines. Thus, we excluded the Mg I*b* line from the fits, preferring to fit the larger surrounding region. Barth et al. (2002b) discuss the \[Mg/Fe\] template matching problem in further detail, using spectra of nearby galaxies obtained during the same observing runs. The best fits for each object are shown in Figure 1.

The dominant contribution to the measurement uncertainty in most cases is template mismatch, since the scatter in $`\sigma `$ measurements for different templates is generally larger than the fitting uncertainty for a single template. The final velocity dispersions listed in Table 2 are the weighted mean of values measured with several templates of type G8–K5. To account for template matching errors, the measurement error for each galaxy is the sum in quadrature of the fitting error for the best-fitting template and the standard deviation of the measurements for all templates.

There is one additional source of systematic uncertainty that contributes significantly to the error budget. As discussed by Barth et al. (2002b), tests on high S/N spectra of nearby galaxies demonstrate that our fitting results for the Mg I*b* region are somewhat sensitive to the order of the polynomial that is used to match the continuum shape of the template to the galaxy. For polynomial orders of 2, 3, or 4, we find variations in the derived value of $`\sigma `$ at the 5% level. To account for this systematic uncertainty in the model fits, we add $`0.05\sigma `$ in quadrature to the uncertainties on each measurement of $`\sigma `$ derived from the blue spectral region. The Ca triplet measurements are much less sensitive to the polynomial order, with variations at the level of $`1\%`$.

Another potential source of systematic error is the variation of the instrumental linewidth across the detector, since a given spectral feature will fall on a different portion of the CCD for different redshifts. This would have the largest effect on the Palomar blue side measurements with the 600 line grating, which was the setting with the lowest spectral resolution, and also with the largest focus variations across the detector. From examination of comparison lamp spectra, we estimate that the maximum variation in instrumental dispersion over the region used for the measurements is $`15`$ km s-1 in this setting. Since the instrumental dispersion contributes in quadrature to the total observed linewidth, these variations will have a very small impact on the final measurements. For the higher-dispersion settings the impact of focus variations on the measurements will be negligible.

### 3.1 Tests of the Measurement Routine

At Palomar, we also observed six velocity dispersion standard galaxies from the catalog of McElroy (1995), with $`\sigma `$ ranging over $`100300`$ km s-1, as a check on the fitting routine. The results are listed in Table 3. The blue and red side measurements for each galaxy agree within their $`1\sigma `$ uncertainty ranges. We find good agreement between our measurements and the average velocity dispersion data compiled by McElroy for four of the six galaxies, while our dispersions for NGC 4579 and NGC 4736 are significantly lower than the McElroy average values. However, in both cases our measurements for these objects do fall in the range of previously measured values as listed in the Hypercat database,111The Hypercat database can be found at http://www-obs.univ-lyon1.fr/hypercat/ . and the close agreement between the blue and red measurements suggests that our values of $`\sigma `$ do give an accurate measure of the velocity dispersion within the observed aperture. We believe that the differences between our results and the mean catalog results for these galaxies are due to individual, discrepant measurements that are included in the catalog average velocity dispersions. A full description of the nearby galaxy measurements is given by Barth et al. (2002b).

Our observations have a wide range in S/N, and one possible concern is that at low S/N the fitting routine could be biased in a systematic way. To test the accuracy of the results, we performed Monte Carlo simulations to determine how well we could recover the correct velocity dispersion as a function of S/N. Starting with the Palomar blue side spectrum of NGC 2841, which was observed with the 1200 line grating, we added a featureless $`f_\lambda \lambda ^1`$ continuum normalized to contribute 2/3 of the total flux at 5500 Å. Then, Poisson noise was added to the spectrum to achieve a desired S/N ratio. For each value of the S/N ratio (ranging from 500 to 50 per pixel), we created 100 realizations of the diluted, noise-added spectrum and measured their velocity dispersions using the same measurement routine and fitting over the wavelength ranges 5050–5150 and 5205–5450 Å. Ten template stars were used for each measurement.

Measurement of the original spectrum gives $`\sigma =229\pm 5`$ km s-1, in good agreement with previous results from the literature (McElroy, 1995). The results of the Monte Carlo calculations are displayed in Figure 2. At each S/N ratio, we computed the weighted average and the standard deviation of the velocity dispersion measurements for the 100 realizations of the spectrum. Over the full range of S/N used in the calculations, the mean velocity dispersion is within $`\pm 3`$ km s-1 of the expected result (229 km s-1), demonstrating that there is no systematic offset as a function of S/N. As expected, the scatter in the measurements rises at low S/N, but for S/N ¿ 100 per pixel the standard deviation is less than 15 km s-1. All but one of our observations have S/N ¿ 100 per pixel, so the expected measurement error is not large. The actual expected error for any measurement is a function of the S/N, dilution, spectral resolution, and velocity dispersion, and we do not attempt to calculate it in more detail, but the simulations demonstrate that the template fitting routine should work well for our data.

### 3.2 Notes on Individual Objects

*Mrk 421:* The Ca triplet lines are easily detected and we found a good fit to the $`\lambda 8542`$ Å profile for $`\sigma =216`$ km s-1. The $`\lambda 8662`$ line was not used in the fit as the telluric H2O absorption beyond 8900 Å could not be adequately removed from the spectrum by standard star division. The blue side spectrum was more strongly diluted by nonthermal light than the other galaxies in the sample, and it proved difficult to obtain an adequate fit. For five template stars, the blue side data yields $`\sigma =213\pm 47`$ km s-1. Thus, we use the red side measurement as our final result.

*Mrk 501:* The Ca II triplet region gives $`\sigma =372\pm 18`$ km s-1, where the uncertainty range includes the scatter among 24 template stars ranging between F3 and K5 and including giants, subgiants, and dwarfs (Paper I). In the blue spectral range, the best-fitting template (HD 125560, a K3III star) gives $`\sigma =386\pm 9`$ km s-1 over 5200-5600 Å, and the standard deviation of results with *all* templates is 81 km s-1. However, if we restrict the fitting to only G and K giants, the scatter is reduced to $`\pm 11`$ km s-1.

During the night that Mrk 501 was observed, the blue side flatfield structure shifted on the detector along the dispersion direction, presumably due to mechanical flexure in the spectrograph. We observed dome and internal flats in the afternoon. Objects observed during the first half of the night flattened well, but in objects observed during the second half of the night (including Mrk 501) there were some residual bumps in the flattened images at the $`10\%`$ level at $`\lambda <5000`$ Å. Although this did not visibly affect the region used to measure $`\sigma `$, and the blue and red side measurements are in agreement, we consider the red side measurement to be more reliable and we use it as our best estimate of $`\sigma `$.

A *Hubble Space Telescope* (*HST*) image of Mrk 501 was taken as part of the snapshot survey of Scarpa et al. (2000) and Urry et al. (2000), but it was not included in the analysis of their *HST* sample. To measure the structural parameters of the host galaxy, we downloaded the images from the *HST* data archive. There are 5 separate exposures with exposure times of 2 s, $`2\times 30`$ s, and $`2\times 120`$ s, all with the F702W filter (similar to ground-based $`R`$). The nucleus is saturated in the 30 s and 120 s exposures. We combined the 120 s exposures, eliminating cosmic-ray hits, and used the galfit program (Peng et al., 2002) to perform a two-dimensional fit to the host galaxy profile. The model used in the fit was a de Vaucouleurs (1948) profile combined with an unresolved nuclear point source, convolved with a synthetic point-spread function (PSF) generated by the Tiny Tim program (Krist & Hook, 2001). Saturated pixels at the nucleus were masked out in the fits to the longer exposure. For consistency with the Urry et al. (2000) sample, we converted the F702W filter magnitude to the Cousins $`R`$ band with the synphot package in IRAF, using a redshifted elliptical galaxy spectrum to calculate the correction. The total host galaxy magnitude is $`m_R=13.07`$ mag, in good agreement with the ground-based measurement of Wurtz, Stocke, & Yee (1996). The fit gives an effective radius of $`r_e=4\stackrel{}{\mathrm{.}}4`$, substantially smaller than previous ground-based measurements (Nilsson et al., 1999); the comparison between ground-based and *HST* measurements of $`r_e`$ is discussed further in §5.

*0521–365 and 0548–322:* During the Magellan run, only one K-giant stellar template was observed (HR 1629, K4III). To allow for template mismatch error, we add 20 km s-1 in quadrature to the uncertainty in the velocity dispersion measurement from this template.

*0706+591:* At $`z=0.125`$, this is the most distant object in the sample. Galactic Na D absorption falls in the region used to measure $`\sigma `$ and is excluded from the fit. The broad H2O telluric absorption band between 5870 and 6000 Å also falls in this region. We could not entirely exclude this band from the fit as it covers the most useful portion of the spectrum; inclusion of part of the telluric band in the fitting region is probably responsible for the relatively poor fit in the neighborhood of the Mg I*b* line. The measurement yields $`\sigma =216\pm 23`$ km s-1. As a check on this result, we also performed a fit over the region 4920–5150 Å (excluding a 20 Å region around \[O III\] $`\lambda 5007`$ although this line was not clearly detected), and found $`\sigma =204\pm 27`$ km s-1.

*2201+044:* The core of the Ca II $`\lambda 8542`$ absorption line falls on a night sky emission line and the extracted spectrum has a sharp negative “spike” at this wavelength. We masked out a 5 Å-wide region surrounding this spike in the final fit, but the measured value of $`\sigma `$ changes by less than 10% if the entire line is included in the fit. Over the wavelength range used for the fit there are several extremely weak telluric absorption features with expected equivalent widths less than 0.02 Å. They were not masked out in the fit, since their amplitudes are smaller than the noise level in the data.

### 3.3 Comparison with the results of Falomo et al.

FKT reported velocity dispersion measurements for seven BL Lac objects, six of which are in our sample. Comparison of the results shows close agreement for some galaxies such as Mrk 421, but also some surprisingly large differences. The worst disagreement is found for Mrk 501 and 1 Zw 187. For these two galaxies, FKT find $`\sigma _c`$ = 291 and 253 km s-1, respectively, while we obtain $`\sigma `$ = 372 and 171 km s-1. FKT suggest that the different results for Mrk 501 could be the result of telluric absorption features in our red spectrum. However, this is not the case because we excluded these telluric lines from the fitting region, as described in Paper I. Different aperture sizes cannot be the cause of such large discrepancies either. FKT used a 1″ slit and extraction widths of $`3\mathrm{}5\mathrm{}`$, similar to the aperture sizes we used, and the resulting variations in $`\sigma `$ should be below the 10% level in most cases (Gebhardt et al., 2000a; Merritt & Ferrarese, 2001).

FKT used the Fourier quotient routine in IRAF which, in principle, should yield the same results on the same data as a direct fitting routine. However, their measurements were carried out over a much wider wavelength range containing different spectral features. They used two grating settings for their measurements: 4800–5800 Å and 5700–8000 Å (setups “A” and “B”, respectively). These ranges include H$`\alpha `$ and H$`\beta `$, which are among the strongest absorption features in the spectrum but which also appear in emission in some objects. They also contain other possible emission features including \[O III\] $`\lambda 5007`$. (In our sample of 11 objects, \[O III\] emission is clearly visible in five, and H$`\beta `$ in three.) FKT removed emission lines from their spectral by connecting a linear continuum across the wavelength range spanned by the emission lines. This is a very poor approximation to the true starlight spectrum of the galaxy, particularly for the Balmer lines, and may adversely effect the results. It would have been preferable to replace the emission lines with the spectra of appropriately broadened template stars (van der Marel & Franx, 1993) or with spectra of galaxies having similar velocity dispersion and line strength (Dalle Ore et al., 1991), or instead to perform the fit over a wavelength range that is devoid of emission lines. Even when H$`\alpha `$ and H$`\beta `$ emission are not clearly visible in the spectra, they could still partially fill in the Balmer absorption lines and affect the measurements.

The strongest absorption feature in the wavelength range used by FKT is the Na D doublet at 5890, 5896 Å, but this could have a contribution from interstellar absorption in the host galaxies. There is Galactic interstellar Na D absorption in some objects as well; our spectra reveal strong Galactic Na D absorption in Mrk 421, and weak absorption in 0706+591. With the Fourier quotient routine there is no entirely satisfactory way to exclude such a feature from the measurement. Since there are few other strong stellar features within a few hundred Å of Na D, this spectral region is considered to be one of the least useful portions of the optical spectrum for measuring velocity dispersions (Dressler, 1984). The weak telluric O2 and H2O absorption bands between 5800 and 6600 Å might also have a small effect on the measurements, especially in objects with a very high degree of dilution by nonthermal light.222See Fig. 16 of Matheson et al. (2000) for a plot of the telluric absorption spectrum over this wavelength range.

Whenever possible, our measurements were performed on regions selected to have good sensitivity to velocity dispersion without having significant contamination from nebular emission, interstellar absorption, or telluric absorption. A major advantage of the direct fitting method over the Fourier quotient method for these objects is that it provides direct visual feedback on the quality of the template match. This makes it possible to mask out regions which are not fit well by any template star, as well as telluric absorption bands and emission lines. We note, also, that our spectra have considerably higher S/N than those obtained by FKT; this enables us to perform our velocity dispersion measurements on a smaller spectral region. In summary, we believe that our velocity dispersion determinations are more accurate than those of FKT.

## 4 Black Hole Masses

To determine $`M_{\mathrm{BH}}`$ from the velocity dispersion, we use the recent fit to the $`M_{\mathrm{BH}}\sigma `$ relation from Tremaine et al. (2002), which is given by

$$\mathrm{log}(M_{\mathrm{BH}}/M_{})=(8.13\pm 0.06)+(4.02\pm 0.32)\mathrm{log}\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right).$$

This relation was derived using a carefully culled sample of the best-quality BH mass measurements, and using a fitting technique that properly accounts for uncertainties both in $`M_{\mathrm{BH}}`$ and $`\sigma `$.

The measurements used to calibrate the above relation were luminosity-weighted velocity dispersions $`\sigma _e`$ measured in a slit aperture of length 2$`r_e`$. However, our values of $`\sigma `$ do not exactly correspond to this aperture size, and our spectra generally do not have sufficient sensitivity to measure $`\sigma `$ outside the inner few arcseconds. Velocity dispersions can be corrected to a standard aperture size using the relations given by Jørgensen, Franx, & Kjærgaard (1995a), but these relations are not calibrated for long-slit apertures such as those used by Tremaine et al. (2002). Another difficulty is that the corrections depend on the effective radii, which are not known accurately for these objects. Rather than apply an uncertain correction to the measured velocity dispersions, we choose to use our measured dispersions directly to estimate $`M_{\mathrm{BH}}`$. As shown by Gebhardt et al. (2000a), the value of $`\sigma `$ measured from slit apertures of different length will differ from $`\sigma _e`$ by less than 10% for typical, nearby elliptical galaxies.

The results are listed in Table 2, with uncertainties in $`M_{\mathrm{BH}}`$ determined by propagating the measurement uncertainty in $`\sigma `$ as well as the uncertainties in the coefficients of the $`M_{\mathrm{BH}}\sigma `$ relation from Tremaine et al. (2002). We do not include the additional possible error due to intrinsic scatter in the $`M_{\mathrm{BH}}\sigma `$ relation, because the amount of scatter is not yet well determined. Tremaine et al. (2002) have determined that this intrinsic scatter is $`0.3`$ dex; it is hoped that future BH mass measurements with *HST* will provide a more definitive result. The masses range from $`7\times 10^7`$ to $`1.6\times 10^9`$ $`M_{}`$, with a median value of $`1.8\times 10^8`$ $`M_{}`$.

We emphasize that these black hole masses are indirect estimates, not dynamical measurements of the mass. Nevertheless, there is no strong reason to believe that the $`M_{\mathrm{BH}}\sigma `$ correlation would not apply to BL Lac objects. Imaging of BL Lacs with *HST* has demonstrated that the host galaxies are all apparently normal giant elliptical galaxies (Scarpa et al., 2000; Urry et al., 2000), in contrast to the disturbed morphologies seen in some QSO host galaxies (e.g., Bahcall et al., 1997). Thus, the stellar velocity dispersions in BL Lac objects should not be affected by recent mergers or other peculiarities. No nearby elliptical galaxies have yet been found to deviate strongly from the $`M_{\mathrm{BH}}\sigma `$ relation, so it would be very surprising if the BL Lac hosts did not follow the correlation. Furthermore, the upper end of the $`M_{\mathrm{BH}}\sigma `$ relation is largely calibrated with FR I radio galaxies, including M87 and NGC 4261, which presumably would appear as BL Lac objects if viewed along the jet, and there is no evidence for any offset in the $`M_{\mathrm{BH}}\sigma `$ correlation of radio galaxies in comparison with inactive ellipticals. Still, it would be a useful consistency check to perform stellar-dynamical measurements of $`M_{\mathrm{BH}}`$ in a few of the lowest-redshift BL Lac objects to confirm that these assumptions are indeed valid.

## 5 Discussion

According to unified models for radio-loud AGNs, BL Lac objects are FR I radio galaxies viewed along the jet axis (e.g., Urry & Padovani, 1995; Urry et al., 2000). If radio galaxies and BL Lac objects were derived from identical parent populations, then the distribution of black hole masses in BL Lac objects should be the same as that in nearby FR I radio galaxies. However, different selection effects for BL Lacs and radio galaxies might be expected to bias the BL Lac objects toward less massive and less luminous host galaxies. For example, Urry et al. (2000) point out that radio-loud AGNs with jets beamed in our direction might not be identified as BL Lac objects if they occur in very bright host galaxies so that their stellar continua are not diluted enough by nonthermal emission to match the standard defining criteria for classification as BL Lac objects. Thus, the naive expectation that BL Lac objects and radio galaxies should have identical distributions of black hole masses or other host galaxy properties is likely to be an oversimplification. Such differences do not violate the basic premise of unified models, as long as objects classified as BL Lacs would still appear to be normal radio galaxies when viewed from an oblique angle.

As a comparison sample of nearby radio galaxies, we use the compilation by Bettoni et al. (2001, 2002). While this is not a complete sample in any sense, it is still the largest sample of nearby radio galaxies for which velocity dispersions, effective radii, and surface brightnesses are available. Figure 3 shows the distribution of velocity dispersions for the BL Lac objects and the radio galaxy sample of Bettoni et al. (2001). As demonstrated previously by FKT, the radio galaxies and BL Lac objects occupy a similar range of values of $`\sigma `$, and therefore $`M_{\mathrm{BH}}`$. The two samples have nearly the same mean velocity dispersion, 247 km s-1 for the BL Lac objects and 256 km s-1 for the radio galaxies. A Kolmogorov-Smirnov test (Press et al., 1988) on the two distributions shows that the hypothesis of identical parent populations for the two distributions can only be rejected at the 9% confidence level. Thus, the two distributions appear to be indistinguishable within the limited statistics, but since sample incompleteness and selection effects are not accounted for in any way, we caution that this comparison does not lead to any firm conclusions regarding the nature of the two populations.

Our sample consists of three low-frequency-peaked (LBL) and eight high-frequency-peaked (HBL) objects. Figure 3 illustrates that, in agreement with FKT, there is no large systematic offset between the black hole masses in LBL and HBL objects. The LBL and HBL objects in our sample have ¡$`\sigma `$¿ = 238 and 235 km s-1, respectively. Similarly, *HST* images do not reveal any differences in the host galaxy magnitudes of these blazar subtypes either (Urry et al., 2000).

The *HST* imaging survey of BL Lac objects by Scarpa et al. (2000) and Urry et al. (2000) demonstrated that the host galaxies appear to be morphologically normal ellipticals. The morphological and dynamical properties of elliptical galaxies are linked via the fundamental plane (FP) relations (Djorgovski & Davis, 1987; Dressler et al., 1987). Therefore, by combining the velocity dispersions with magnitudes and effective radii, we can test whether the BL Lac host galaxies fall on the same FP as normal ellipticals. For consistency with Urry et al. (2000) and Bettoni et al. (2001), we use $`H_0=50`$ km s-1 Mpc-1.

One source of uncertainty in such a comparison comes from the fact that the effective radii measured from ground-based images and from the *HST* snapshot survey of Scarpa et al. (2000) and Urry et al. (2000) are in severe disagreement for some of the galaxies in our sample. Table 4 lists the effective radii from the *HST* survey, and from ground-based imaging from various sources in the literature. The ground-based measurements of $`r_e`$ are nearly always larger than the *HST* measurements. In the worst case (3C 371), the disagreement is a factor of 5.

The main advantage of the *HST* images is the narrow PSF core which makes it possible to determine the point-source flux accurately, but the *HST* measurements are hindered by the small field of view of the WFPC2 Planetary Camera (PC) detector (35″ square), and brief exposures of 200 to 1300 s. High-quality ground based imaging (e.g., Nilsson et al., 1999) demonstrates that the host galaxies of the nearest objects are indeed large enough to completely fill the PC chip. As a result, the sky background level determined from the outer portions of the detector will be contaminated by galaxy light, and measurements of the galaxy radial profile after subtraction of an artificially high sky background would lead to an underestimate of $`r_e`$. This systematic effect could be largely responsible for the differences between the ground-based and *HST* measurements. The total host galaxy magnitudes are in better agreement, with differences of typically less than a few tenths of a magnitude. The problem of determining the host galaxy parameters for BL Lac objects from shallow *HST* images is somewhat analogous to the difficulty of detecting QSO host galaxies in short WFPC2 exposures, as discussed by McLeod & Rieke (1995). Nevertheless, it is clear that the ground-based measurements of $`r_e`$ are quite uncertain as well. The case of Mrk 501 illustrates the difficulty of measuring $`r_e`$ from the ground. Even though this is one of the nearest and best-resolved BL Lac objects, ground-based measurements of $`r_e`$ range from 9″ to 20″.

Figure 4 plots the FP for the BL Lac objects using the parameterization given by Bettoni et al. (2001) for low-redshift radio galaxies. Given the uncertainty in the measurements of $`r_e`$, and the consequent uncertainty in $`\mu _e`$ (the average $`R`$-band surface brightness enclosed within $`r_e`$), we plot the BL Lac objects using both the ground-based and *HST* measurements. The host galaxy magnitudes have been corrected for Galactic extinction and $`k`$-corrected as described by Urry et al. (2000). The ground-based data included in the plot are from Abraham, McHardy, & Crawford (1991), Wurtz, Stocke, & Yee (1996), and Nilsson et al. (1999). For comparison, the radio galaxy sample of Bettoni et al. (2001) and the Coma cluster ellipticals observed by Jørgensen, Franx, & Kjærgaard (1995a, b) are also shown.

The overall agreement is very reasonable, considering the relatively large uncertainties in the measurements of $`\sigma `$ and $`r_e`$ for the BL Lac objects. Changes in $`r_e`$ at fixed total magnitude result in displacement along a vector that is nearly parallel to the FP sequence; as a result both the ground-based and *HST*-based measurements appear consistent with the FP. In either case, most of the BL Lac objects fall within the scatter of the FP. The largest outlier is 3C 371 when plotted with the effective radius of 2.9 kpc from Urry et al. (2000); it deviates from the FP sequence by $`0.3`$ dex. However, ground-based measurements indicate an effective radius that is larger by a factor of five, and with the ground-based value of $`r_e`$ this galaxy falls nicely on the FP. Thus, despite the uncertainty in the effective radii, the host galaxies do appear to be normal elliptical galaxies.

The contrast between the ground-based and *HST*-based measurements of $`r_e`$ and $`\mu _e`$ can be seen more clearly by plotting projections of the FP onto the $`r_e`$ and $`\mu _e`$ axes, as shown in Figure 5. With the *HST* measurements, the BL Lac objects have effective radii comparable to the smallest radio galaxies, and most of the BL Lac hosts have $`\mu _e`$ brighter than 19 mag arcsec-2, in a region where very few radio galaxy hosts lie. The ground-based measurements place the BL Lac hosts in a region that overlaps much more closely with the radio galaxies, better matching expectations from unified models. However, the BL Lac objects are still displaced toward smaller radii than the largest radio galaxies; this is also apparent in Figure 4, where the BL Lac objects do not extend as far up the FP sequence as the radio galaxies. Similarly, Urry et al. (2000) noted that no BL Lac hosts are found to be as luminous as the most luminous hosts of FR I radio galaxies. This offset could be a consequence of selection effects, or might indicate an environmental effect disfavoring the formation of blazar jets in rich clusters (Urry et al., 2000).

The results presented here, and by FKT, strongly disagree with several recent estimates of $`M_{\mathrm{BH}}`$ in blazars based on other methods. Using variability timescales to estimate the size of the emission region, Xie et al. (1998) and Fan, Xie, & Bacon (1999) derive estimates of the black hole masses for several blazars that are in the range $`10^610^8`$ $`M_{}`$, much smaller than the masses obtained from the $`M_{\mathrm{BH}}\sigma `$ relation. Wang, Xue, & Wang (2002) have proposed a method to derive $`M_{\mathrm{BH}}`$ from the peak luminosity and frequency of the blazar emission. They conclude that HBL and LBL blazars have a bimodal distribution of $`M_{\mathrm{BH}}`$, with masses of $`10^510^8`$ $`M_{}`$ for HBL and $`10^910^{11}`$ $`M_{}`$ for LBL objects. The stellar velocity dispersions imply black hole masses that are are completely inconsistent with this conclusion. Models that interpret periodicities in blazar light curves in terms of orbital motion of black hole binary systems (Rieger & Mannheim, 2000; De Paolis, Ingrosso, & Nucita, 2002) also give much smaller masses than those derived from the $`M_{\mathrm{BH}}\sigma `$ relation.

Xie et al. (2002) propose a new method to estimate black hole masses in BL Lac objects, based on the assumption that the observed variability timescale is determined by the orbital timescale near the innermost stable orbit around the black hole. They find masses ranging from $`10^{7.4}`$ to $`10^{9.2}`$ for 13 objects, in the same range as the masses we derive. Xie et al. show that their $`M_{\mathrm{BH}}`$ estimates are strongly correlated with the beamed nonthermal luminosity but not with the host galaxy luminosity. However, the opposite would be expected if BL Lac host galaxies follow the $`M_{\mathrm{BH}}L_{\mathrm{bulge}}`$ correlation of normal galaxies. Furthermore, Urry et al. (2000) show that there is no correlation between the nuclear and host galaxy luminosities for BL Lac objects, implying that the nonthermal luminosity should not be tightly correlated with $`M_{\mathrm{BH}}`$ if $`M_{\mathrm{BH}}`$ scales with host galaxy luminosity. Only two objects are in common between our sample and Xie et al., so it is not possible to perform a detailed comparison of the mass estimates, but we note that their upper limit to the black hole mass in 3C 371 is an order of magnitude lower than the mass we derive from the $`M_{\mathrm{BH}}\sigma `$ correlation.

Wu, Liu, & Zhang (2002) and Woo & Urry (2002) have recently published estimates of the BH masses in BL Lac objects by using the *HST* imaging survey results to obtain approximate velocity dispersions via the FP relation, and then applying the $`M_{\mathrm{BH}}\sigma `$ relation to estimate $`M_{\mathrm{BH}}`$. The advantage of this method is that it can be easily applied to a large number of objects, as there are over 60 BL Lac objects with morphological parameters measured from the *HST* survey. However, there are some potential disadvantages as well, since the $`\sigma `$ estimates are limited by the accuracy of the host galaxy decompositions and by the intrinsic scatter in the FP. In both cases, the FP-derived estimates of $`\sigma `$ are systematically larger than our measured values. For objects in common with our sample, the Wu et al. estimates of $`\sigma `$ are higher than our values by 30% on average, and the Woo & Urry estimates are higher by 74%. The most severe disagreement is for 3C 371; for this object Wu et al. predict $`\sigma =524`$ km s-1 and Woo & Urry predict 618 km s-1, in comparison with our measured value of $`249\pm 25`$ km s-1. This illustrates the need for accurate values of $`r_e`$ and $`\mu _e`$ when applying this method.

Determining the Eddington ratio $`L/L_{\mathrm{Edd}}`$ in BL Lac objects can give important clues to the physical structure of their accretion flows. It would be interesting to search for systematic differences between LBL and HBL objects; this could yield new constraints for models of the spectral energy distributions (e.g., Ghisellini et al., 1998; Böttcher & Dermer, 2002). Unfortunately, it is particularly difficult to measure the bolometric luminosity of the unbeamed AGN core in BL Lac objects that are dominated by the beamed nonthermal emission from the jet. One alternative is to use the emission lines as a measure of the ionizing luminosity, because the emission-line flux is the only unbeamed component of the AGN that can be directly observed. This method is not ideal since there is no single emission line that can be observed in all objects; in practice the only viable alternative is to use the lines that are detected and perform a rough conversion either to a total broad-line region luminosity (Celotti, Padovani, & Ghisellini, 1997) or to a single line such as H$`\beta `$ (Cao, 2002). For our sample, broad-line fluxes are available only for a few objects, while a few others have \[O III\] emission detected. Given this small number of objects and the large uncertainties involved in comparing the luminosities of different emission lines, we are unable to perform any meaningful comparison of the Eddington ratios of the HBL and LBL objects in our sample. This would be a useful objective for future work if stellar velocity dispersions and emission-line luminosities can be measured for larger samples of BL Lac objects.

## 6 Conclusions

Our conclusions are summarized as follows:

1. Stellar velocity dispersions in low-redshift BL Lac objects are in the range $`170370`$ km s-1. Using the Tremaine et al. (2002) fit to the $`M_{\mathrm{BH}}\sigma `$ relation, the corresponding black hole masses are $`M_{\mathrm{BH}}10^{7.9}`$ to $`10^{9.2}`$ $`M_{}`$.

2. Measurement of velocity dispersions in BL Lac objects requires particular care in selecting spectral regions that are sensitive to $`\sigma `$ but are not severely affected by interstellar absorption lines, telluric absorption bands, or emission lines. Due to the dilution by nonthermal emission, high S/N (typically $`100`$ per pixel in the extracted spectrum) is required for accurate measurements with a direct template-fitting method.

3. The distribution of velocity dispersions in BL Lac objects appears superficially similar to that of the radio galaxy sample of Bettoni et al. (2001), indicating a similar distribution of black hole masses for the two samples. However, neither sample is statistically well defined or complete in any sense and we caution against drawing conclusions about radio-loud AGNs in general from this limited comparison.

4. There does not appear to be any systematic difference between the black hole masses in HBL and LBL objects.

5. The host galaxies of BL Lac objects lie on the fundamental plane of nearby elliptical galaxies. We do not find any BL Lac objects near the top of the fundamental plane sequence in the region occupied by the most luminous hosts of FR I radio galaxies; this may be a consequence of selection effects in the identification of BL Lac objects.

6. For nearby, well-resolved BL Lac objects, the host galaxy effective radii measured in the *HST* WFPC2 snapshot survey of Scarpa et al. (2000) and Urry et al. (2000) are systematically smaller than effective radii measured from deeper ground-based images. We suggest that the measurements of $`r_e`$ from the *HST* snapshot survey images may be biased toward underestimates of the true radii, because they were derived from shallow exposures that did not detect the faint outer envelopes of the host galaxies that are seen in deep ground-based images.

We thank Swara Ravindranath for assistance with the observing runs at Las Campanas and with the *HST* imaging analysis, Tom Matheson for providing some of the spectroscopic reduction software, and Meg Urry for helpful discussions. Research by A.J.B. is supported by NASA through Hubble Fellowship grant #HST-HF-01134.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. Research by W.L.W.S. is supported by NSF grant AST-9900733. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research was partly based on observations made with the NASA/ESA *Hubble Space Telescope*, obtained from the data archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.


format: text

text:
# Quantum Corrections to Dilute Bose Liquids

## Abstract

It was recently shown \[A. Bulgac. Phys. Rev. Lett. 89, 050402 (2002)\] that an entirely new class of quantum liquids with widely tunable properties could be manufactured from bosons, fermions and their mixtures by controlling their interaction properties by means of a Feshbach resonance. Since the quantum fluctuations in this regime could in principle destabilize these objects, we extend the previous mean–field analysis of these quantum liquids by computing the lowest order quantum corrections to the ground state energy and the depletion of the Bose–Einstein condensate and by estimating higher order corrections as well. We show that the quantum corrections are relatively small, do not lead to the destabilization of the droplets and are controlled by the diluteness parameter $`\sqrt{n|a|^3}1`$, even though strictly speaking in this case there is no low density expansion.

Among the quantum fluids observed so far only 3He and 4He and atomic nuclei are self–bound/liquid systems. It is very likely that this relatively short list could become extremely long by the addition of boson droplets (boselets), fermion droplets (fermilets) and their mixtures (ferbolets) , which owe their existence to the possibility of tuning two–atom interactions and to some unusual properties of the three–body system.

It was predicted about ten years ago  and was confirmed experimentally about five years ago  that by immersing atoms in a magnetic field (and, in principle, in electric or laser fields as well), one can alter essentially at will the scattering length $`a`$ between two atoms. Under such circumstances, it is possible to create a situation where the scattering length is negative and significantly larger than the interaction radius. In the regime where the magnitude of the scattering length is much larger than the typical scale of atomic interactions $`r_0`$, the three–atom system exhibits some very unusual and interesting phenomena . The three–body observables depend on a new scale $`a_3`$ which is not determined by any two-body observable. In the language of effective field theory, this corresponds to the fact that a three-body force appears at leading order of the low energy expansion . Also, three–body observables are periodic functions of $`\mathrm{ln}(a/a_3)`$ and the three–particle scattering amplitude can be made arbitrarily large by fine–tuning $`a`$. By choosing $`a`$ negative, large ($`|a|r_0`$) and with a three–body bound state close to the threshold, a rather unique situation is created: two low energy atoms experience an effective attraction, but at the same time three atoms will experience an effective repulsion. This situation is very well understood theoretically and it is known as the Efimov effect . A large ensemble of such atoms can condense in the same way as water droplets condense from vapor. Since such condensation can occur for any atom species for which a Feshbach tunable resonance exist, an entirely new class of (bose, fermi, or mixtures) quantum liquids  can be created. Their equilibrium density is determined by the interplay between the two–body attraction and the three–body repulsion between particles in the Efimov regime described above. Since the magnitudes of the two–body attraction and of the three–body repulsion are to a significant extent under experimental control, the basic properties of these new quantum liquids are widely tunable .

In Ref. , the basic properties of this new class of quantum fluids were established in the mean field approximation. It is not at all obvious that quantum fluctuations could not destabilize a boselet. It is generally expected that the mean field approximation is rather accurate for dilute systems . Under normal circumstances (when $`a>0`$ and $`ar_0`$) it is sufficient to take into account only the binary atom–atom collisions and the effects of quantum fluctuations, whereas the triple atom and higher collisions are rather small and controlled by the diluteness parameter $`\sqrt{na^3}1`$ . The regime $`a<0`$ was considered until recently as intrinsically unstable towards collapse, unless the number of bosons is smaller than about 1,500 or so in a trap . In the specific Efimov regime, we are interested in the fact that the two–body scattering amplitude is large ($`a<0,|a|r_0`$), but the three–body amplitude is even larger. Even though the system is dilute with respect to the two–body collisions ($`n|a|^31`$), it is not obvious that diluteness with respect to the three–body collisions is achieved and one might suspect that the mean–field approximation could be violated. It is imperative to understand the character and the magnitude of these corrections to the ground state properties in the case of these new quantum liquids and determine whether quantum fluctuations can lead to instabilities. One can naturally expect that such effects would perhaps manifest themselves particularly strongly in the case of bosons and less so in the case of fermions. For this reason, we shall focus our attention here on Bose systems only.

We start with the Lagrangian density

$`+\mu 𝒩=\psi ^{}\left(i\mathrm{}{\displaystyle \frac{}{t}}+{\displaystyle \frac{\mathrm{}^2\mathbf{}^2}{2m}}+\mu \right)\psi `$ (1)
$`{\displaystyle \frac{g_2}{2}}\psi ^{}\psi ^{}\psi \psi {\displaystyle \frac{g_3}{6}}\psi ^{}\psi ^{}\psi ^{}\psi \psi \psi ,`$ (2)
where $`\psi ^{}`$ and $`\psi `$ are the creation and annihilation operators for a boson, $`\mu `$ is the chemical potential, and the coupling constants are

$`g_2`$ $`=`$ $`{\displaystyle \frac{4\pi \mathrm{}^2a}{m}},`$ (3)
$`g_3`$ $`=`$ $`{\displaystyle \frac{12\pi \mathrm{}^2a^4}{m}}\left[d_1+d_2\mathrm{tan}\left(s_0\mathrm{ln}{\displaystyle \frac{a}{a_3}}+{\displaystyle \frac{\pi }{2}}\right)\right]`$ (4)
$`=`$ $`{\displaystyle \frac{6\pi \mathrm{}^2a^4}{m}}\mathrm{{\rm Y}}.`$ (5)
Here $`a`$ is the two–body scattering length and $`s_01.00624`$. $`d_1`$ and $`d_2<0`$ are universal constants, whose numerical values have been determined recently . The values of the coupling constant are determined by considering the scattering of two and three particles at zero momentum. $`a_3`$ is the value of the two–body scattering length for which a three–body bound state has exactly zero energy. Unlike $`d_{1,2}`$ and $`s_0`$, the parameter $`a_3`$ is system dependent and is also a genuine three–body characteristic. The rest of the symbols have their usual meaning. We have introduced a new dimensionless quantity $`\mathrm{{\rm Y}}1`$ in Eq. (5), which will prove very convenient in the power counting. We will only use the Lagrangian in Eq. (2) for particle momenta much smaller than $`\mathrm{}/|a|`$ so that it is legitimate to subsume the complicated dynamics on the scale $`\mathrm{}/a`$, leading to the Efimov effect into a contact three–body interaction. We show below that the typical loop momenta are of order $`Q\mathrm{}/(|a|\mathrm{{\rm Y}}^{1/2})\mathrm{}/|a|`$. More precision can be systematically attained by including terms in Eq. (2) with more derivatives or fields, but their effect is strongly suppressed. As we mentioned above, the interesting regime is when $`g_2<0`$ and $`g_3>0`$ .

Since we are interested in the condensed state of a system of bosons, it is natural to split the field $`\psi `$ into a classical ($`c`$–number) condensate $`\varphi `$ and the fluctuations $`\psi `$: $`\psi \psi +\varphi `$. ($`\varphi `$ can be chosen real in this case.) In terms of the new variables, the Lagrangian density becomes

$`+\mu 𝒩={\displaystyle \frac{1}{2}}g_2\varphi ^4{\displaystyle \frac{1}{6}}g_3\varphi ^6+\mu \varphi ^2+`$ (6)
$`\varphi (\psi +\psi ^{})\left(\mu g_2\varphi ^2{\displaystyle \frac{g_3}{2}}\varphi ^4\right)+`$ (7)
$`{\displaystyle \frac{1}{2}}(\psi ^{}\psi )\left(\begin{array}{cc}i\mathrm{}\frac{}{t}+\frac{\mathrm{}^2\mathbf{}^2}{2m}+\chi & g_2\varphi ^2g_3\varphi ^4\\ g_2\varphi ^2g_3\varphi ^4& i\mathrm{}\frac{}{t}+\frac{^2}{2m}+\chi \end{array}\right)\left(\begin{array}{c}\psi \\ \psi ^{}\end{array}\right)`$ (12)
$`+(\mathrm{interactions}),`$ (13)
where $`\chi =g_2\varphi ^2+g_3\varphi ^4`$. The interaction terms include vertices with three $`\psi `$’s ($`g_2\varphi `$ and $`g_3\varphi ^3`$), four $`\psi `$’s ($`g_2`$ and $`g_3\varphi ^2`$), five $`\psi `$’s ($`g_3\varphi `$), and six $`\psi `$’s ($`g_3`$) and normal ordering is implied as well. At tree level/mean field, the particle number density is $`n=\varphi ^2`$ and the energy density and chemical potential (determined by imposing the condition $`\psi =0`$) are given by

$`={\displaystyle \frac{g_2}{2}}n^2+{\displaystyle \frac{g_3}{6}}n^3,`$ (14)
$`\mu ={\displaystyle \frac{d}{dn}}=g_2n+{\displaystyle \frac{g_3}{2}}n^2.`$ (15)
In writing down the contribution quadratic in $`\psi `$’s in Eq. (13) we have implicitly taken into account the condition (15). At equilibrium/zero pressure ($`P=nd/dn=0`$) the values for the particle number density, energy density, and chemical potential are easily determined to be

$`n_0={\displaystyle \frac{1}{|a|^3}}{\displaystyle \frac{1}{\mathrm{{\rm Y}}}},`$ (16)
$`_0=\mu _0n_0={\displaystyle \frac{\mathrm{}^2}{ma^2}}{\displaystyle \frac{\pi }{\mathrm{{\rm Y}}}}n_0,`$ (17)
$`\mu _0={\displaystyle \frac{d}{dn}}|_{n=n_0}={\displaystyle \frac{_0}{n_0}}={\displaystyle \frac{\mathrm{}^2}{ma^2}}{\displaystyle \frac{\pi }{\mathrm{{\rm Y}}}}.`$ (18)
Notice that only at equilibrium the chemical potential becomes equal to the energy per particle, see Eq. (18). The tree level/mean field approximation is valid when the system is dilute, which holds as long as $`\mathrm{{\rm Y}}1`$. At and near equilibrium, the contribution of the two–body and three–body collisions are comparable in magnitude. The three–body collision term, which in dilute systems is typically negligible, is important in this case, because the three particle bound state is close to the threshold and the particles are in the Efimov regime.

From the Lagrangian (13) we can determine the normal and anomalous propagators for the Bogoliubov quasiparticle, which can be conveniently represented as a matrix

$$D_k=\frac{1}{\omega ^2\omega _k^2}\left(\begin{array}{cc}\omega +\frac{\mathrm{}^2𝒌^2}{2m}+\chi & \chi \\ \chi & \omega +\frac{\mathrm{}^2𝒌^2}{2m}+\chi \end{array}\right),$$
(19)
where $`\omega _k`$ is the Bogoliubov dispersion relation for quasiparticle excitations,

$$\omega _k=\sqrt{\frac{\mathrm{}^2𝒌^2}{2m}\left(2\chi +\frac{\mathrm{}^2𝒌^2}{2m}\right)}\mathrm{}k\sqrt{\frac{\chi }{m}}.$$
(20)
Since $`g_2<0`$, the sound velocity $`s=\sqrt{\chi /m}`$ for long wavelengths is imaginary for small particle number densities, namely when

$$n<n_s=\frac{|g_2|}{g_3}=\frac{2}{3}n_0,$$
(21)
and the homogeneous matter at these densities ($`n<n_s`$) is unstable towards collapse and condensation into denser droplets. For densities $`n_snn_0`$, the pressure is negative and such a system tends to increase its density and binding energy. If the density is $`n>n_0`$, the internal pressure is positive and the system tends to expand towards equlibrium and lowers its energy, unless external walls exert inward pressure. Unlike the phonons in liquid 4He, the curvature of the quasiparticle excitation $`\omega _k`$ as a function of the wave vector $`k`$ is positive. As a result these quasiparticles have a finite lifetime due to decay into two quasiparticles, corresponding to a width $`\mathrm{\Gamma }k^5`$ , since the velocity of a quasiparticle with finite wave vector $`k`$ is larger than the sound velocity.

The quantum fluctuations can now be evaluated at the one–loop level with usual techniques. We have chosen to use dimensional regularization for the evaluation of various diverging integrals . Old fashioned techniques , such as diagonalization of the quadratic part of the Lagrangian (13) with the a Bogoliubov transformation and subsequent evaluation of the corrections to the ground state density and energy would lead to identical results. The particle density and the energy density at one–loop level are given by the following expressions

$`n=\varphi ^2+{\displaystyle \frac{1}{3\pi ^2}}\left({\displaystyle \frac{m\chi }{\mathrm{}^2}}\right)^{3/2},`$ (22)
$`={\displaystyle \frac{g_2\varphi ^4}{2}}+{\displaystyle \frac{g_3\varphi ^6}{6}}+{\displaystyle \frac{8\mathrm{}^2}{15\pi ^2m}}\left({\displaystyle \frac{m\chi }{\mathrm{}^2}}\right)^{5/2}+\mu (n\varphi ^2).`$ (23)
The first correction term to the energy density (and similarly to the number density as well) is formally identical to the Lee and Yang term , which was computed, however, for a Bose gas with repulsive two–body interaction. The equilibrium values for the particle number density (both condensed and noncondensed), the energy density, and energy per particle, after including first–order quantum corrections become

$`n={\displaystyle \frac{1}{|a|^3}}{\displaystyle \frac{1}{\mathrm{{\rm Y}}}}\left(1{\displaystyle \frac{72}{5}}\alpha \right),\mathrm{where}\alpha ={\displaystyle \frac{2^{3/2}}{3\pi ^{1/2}\mathrm{{\rm Y}}^{1/2}}},`$ (24)
$`=_0\left(1{\displaystyle \frac{88}{5}}\alpha \right),{\displaystyle \frac{}{n}}=\mu _0\left(1{\displaystyle \frac{16}{5}}\alpha \right).`$ (25)
Surprisingly, Eqs. (22) and (23) determine correctly both the linear and the quadratic corrections in $`\alpha `$ for energy per particle $`/n`$. Therefore, the first order quantum corrections are controlled by the same diluteness parameter $`\sqrt{n|a|^3}=\mathrm{{\rm Y}}^{1/2}1`$ as in the case of a Bose gas with repulsive interaction. In spite of this formal apparent similarity, the meaning of the present result for the condensed Bose liquid when $`a<0`$ is qualitatively different. There is strictly no low–density expansion in the present case, since for densities $`n<n_s=2n_0/3`$, the system is unstable towards long–wavelength density fluctuations as the sound velocity is imaginary, see discussion of Eqs. (20) and (21).

A general argument can be given for the suppression of higher loop corrections. Consider a generic diagram contributing to the energy density (that is, without external legs) containing I propagators, L loops and $`n_i`$ vertices with $`i`$ legs. It can be estimated by

$`(\mathrm{Energy}\mathrm{Diagram})\left({\displaystyle \frac{m}{Q^2}}\right)^I\left(\sqrt{{\displaystyle \frac{\chi }{m}}}QQ^3\right)^L`$ (26)
$`\times \left({\displaystyle \frac{\chi }{n^{1/2}}}\right)^{n_3}\left({\displaystyle \frac{\chi }{n}}\right)^{n_4}\left({\displaystyle \frac{\chi }{n^{3/2}}}\right)^{n_5}\left({\displaystyle \frac{\chi }{n^2}}\right)^{n_6}`$ (27)
$`m^{3L/2}\chi ^{I+5L/2+n_3+n_3+n_5+n_6}`$ (28)
$`\times n^{n_3/2n_43n_5/22n_6},`$ (29)
where $`Q\sqrt{m\chi }1/(|a|\mathrm{{\rm Y}}^{1/2})`$ is the typical loop momentum . When using dimensional regularization no powers of the cutoff are present and the only momentum scale in these diagrams is set by $`\sqrt{m\chi }`$. Combining the above estimate with the relations

$`2I`$ $`=`$ $`3n_3+4n_4+5n_5+6n_6,`$ (30)
$`I`$ $`=`$ $`n_3+n_4+n_5+n_6+L1,`$ (31)
we find

$$(\mathrm{Energy}\mathrm{Diagram})m^{\frac{3L}{2}}\chi ^{\frac{3L}{2}+1}n^{1L}_0\mathrm{{\rm Y}}^{L/2}.$$
(32)
Consequently, up to logarithmic corrections, all higher–order loop corrections are suppressed in the dilute limit by powers of $`\sqrt{n|a|^3}=\mathrm{{\rm Y}}^{1/2}1`$. As in the case of dilute Bose gases with repulsive interactions  this parameter is proportional to the ground–state fraction of non–BEC particles. In this respect, this new class of quantum liquids, in particular the boselets, is qualitatively different from liquid 4He, where the fraction of non–BEC particles is close to 90%.

For large $`\mathrm{{\rm Y}}1`$, the self–bound system is dilute. The states with positive pressure, however, can be dense. These states could still be dilute with respect to two–body collisions, however of such a density that $`|g_2|g_3n`$. In this regime, the loop expansion breaks down and a further resummation is required. The energy per particle of the system is large and the system can exist only if external pressure is applied. It is very likely that the system undergoes a transition to a coherent mixture of atoms and trimers, similar to the situation put in evidence in gaseous BEC, where the coexistence of an atomic BEC with a molecular BEC was predicted  and experimentally observed . The formation of trimers will have two effects, lowering the particle density and lowering the energy. The spatial size of a trimer when $`a<0`$ is smaller than $`|a|`$, even if the trimer binding energy is vanishing. The wave function of an Efimov state of zero energy is normalizable and the wave function is concentrated predominantly at interparticle distances $`|𝒓_i𝒓_j|<|a|`$, where $`i,j=1,2,3`$ are particle labels in a trimer . Since the average separation between atoms is presumed to be larger than the scattering length, when all atoms convert into trimers, the trimer number density decreases roughly by a factor of 3 when compared to the atomic number density . The trimer phase is most likely stable, provided the trimer–trimer scattering length is positive, which is most probable . When a trimer is formed, a significant amount of energy is released as well. These two factors, the significant drop in density and energy, should be explicitly taken into account when evaluating the energy density at higher particle densities, where triple collisions might dominate. The self–bound state we have considered is only meta–stable. It can decay either by a four particle process, leading to the formation of a trimer and one energetic atom, or through a three–body process with the formation of a deeper dimer and one energetic atom (the formation of a shallow dimer is excluded since $`a<0`$). The rate for the first process is proportional to $`\mathrm{}n^4a^7/m`$ and likely very small. The formation of a deep dimer was considered in Ref.  and it was found to be of the order of $`\mathrm{}n^3a^4/m`$. This decay rate is, however, enhanced when a three–body state is close to the threshold, where the correct rate can be determined only through a non–perturbative calculation, but this is not yet attempted. This recombination rate is related to the widths of the Efimov states due to the existence of deep two–body bound states, which follow the same exponential behavior as the energies of the Efimov states .

We thank G.F. Bertsch and D.B. Kaplan for discussions.

Note added. Boselets made of spin polarized tritium atoms are the likely candidates, since the regime we study here, when $`g_2<0`$ and $`g_3>0`$, could apparently be easily reached according to Ref. . In this particular case, three–body recombination processes are absent.

Note added in proof. Recently Braaten and Hammer  have performed a nonperturbative calculation of three–body recommbination rates in dilute gases.


text:
# Two-Component Approach to 𝐽/Ψ Production in High-Energy Heavy-Ion CollisionsThis work was supported by the U.S. Department of Energy under Grant No. DE-FG02-88ER40388.

## 1 Introduction

In the standard picture of $`J/\mathrm{\Psi }`$ production in heavy-ion collisions , charmonium states are (pre-) formed in primordial $`N`$-$`N`$ collisions and subject to subsequent nuclear absorption, Quark-Gluon Plasma (QGP) and hadronic dissociation. Hence, the magnitude of $`J/\mathrm{\Psi }`$ suppression reflects on the plasma effect provided that nuclear and hadronic suppression are known and/or small. Recently, an alternative mechanism of $`J/\mathrm{\Psi }`$ production in relativistic heavy-ion collisions has been suggested . The underlying idea is that charmonium states are formed by statistical coalescence of $`c`$ and $`\overline{c}`$ quarks at the hadronization transition according to thermal weights, with the total $`c\overline{c}`$ abundance being determined by primordial parton fusion.

In the following, we present a combined approach of direct and statistical $`J/\mathrm{\Psi }`$ production , constructed within a thermal framework that is consistent with basic hadro-chemistry and expansion dynamics at both SPS and RHIC, and also reproduces electromagnetic observables at the SPS. In particular, we do not invoke any “anomalous” open-charm enhancement beyond expectations from $`N`$-$`N`$ collision scaling.

## 2 Two-component model of $`J/\mathrm{\Psi }`$ production

The first (“direct”) component of $`J/\mathrm{\Psi }`$ yield is attributed to primordial $`J/\mathrm{\Psi }`$’s arising from early (hard) parton-parton collisions. A first suppression factor, $`𝒮_{nuc}`$, is due to nuclear absorption evaluated in the Glauber picture with a phenomenological (constant) cross-section $`\sigma _{\mathrm{\Psi }N}6.4`$ mb. Subsequent QGP suppression is based on inelastic scattering of thermal partons off the $`c`$ or $`\overline{c}`$ quark within the bound state, using a quasifree approximation with account for in-medium modified charmonium binding energies . Hadronic dissolution cross sections are estimated within a $`SU(4)`$-symmetric effective theory . The pertinent dissociation rates for both QGP and hadron gas (HG) are convoluted over a thermal fireball evolution  for a heavy-ion collision at given impact parameter, yielding charmonium suppression factors $`𝒮_{QG}^X`$ and $`𝒮_H^X`$ ($`X`$= $`J/\mathrm{\Psi }`$, $`\psi ^{}`$, $`\chi `$), respectively. This leads to the number of direct $`J/\mathrm{\Psi }`$’s from primordial production (including feeddown),

$$𝒩_{J/\mathrm{\Psi }}^{dir}=\sigma _{pp}^{J/\mathrm{\Psi }}ABT_{AB}(b)𝒮_{nuc}\left[0.6𝒮_{QG+H}^{J/\mathrm{\Psi }}+0.032𝒮_{QG+H}^\chi +0.08𝒮_{QG+H}^\mathrm{\Psi }^{}\right].$$
(1)
The second (“statistical”) component of $`J/\mathrm{\Psi }`$ yield originates from coalescence of $`c`$ and $`\overline{c}`$ quarks at hadronization. The abundances of charmed particles follow from hadronic thermal weights at the critical temperature $`T_c`$, with the absolute number matched to the amount of charm quarks created in primordial $`N`$-$`N`$ collisions. This necessitates the introduction of an effective fugacity factor $`\gamma _c`$ which is determined enforcing exact (local) charm conservation. The statistical $`J/\mathrm{\Psi }`$ production, including feeddown, is then given by

$$𝒩_{J/\mathrm{\Psi }}^{th}=\gamma _c^2V_H\left[n_{J/\mathrm{\Psi }}𝒮_H^{J/\mathrm{\Psi }}+(XJ/\mathrm{\Psi })n_X𝒮_H^X\right],$$
(2)
where the volume at hadronization, $`V_H`$, is taken from our fireball evolution. The factor $``$$``$1 approximates the effects of incomplete charm-quark thermalization (calculated in a relaxation time approach), inducing a relative reduction of statistical production.

The total number of $`J/\mathrm{\Psi }`$’s observed per heavy-ion collision follows as $`𝒩_{J/\mathrm{\Psi }}=𝒩_{J/\mathrm{\Psi }}^{dir}+𝒩_{J/\mathrm{\Psi }}^{th}`$. The model involves essentially two parameters: the strong coupling constant figuring into the QGP dissociation cross section, as well as the $`c`$-quark thermalization time.

## 3 SPS

Our results for the centrality dependence of the $`J/\mathrm{\Psi }`$ yield in the $`Pb`$(158 AGeV)-$`Pb`$ system are compared to NA50 data  in Fig. 2. Direct production (dashed line) always prevails over the thermal component (dot-dashed line). The latter sets in once a QGP starts forming, which, coupled with the rather strong QGP suppression, mimics the first drop at around $`E_T`$$``$40 GeV. The overall good description of the NA50 data extends beyond $`E_T>100`$ GeV (not shown here) when incorporating transverse energy fluctuations  and a careful assessment of the minimum bias analysis , cf. ref. . In addition, the model also reproduces the NA38 results for the $`S`$-$`U`$ system reasonably well. We note that the relative smallness of statistical production stems from not invoking any “anomalous” open-charm enhancement.

Within our two-component model, the $`\mathrm{\Psi }^{}/\mathrm{\Psi }`$ ratio can be extracted without further assumptions. The results (Fig. 2) for both $`S`$-$`U`$ (dashed) and $`Pb`$-$`Pb`$ (solid) indicate significant discrepancy with NA38/NA50 data. We believe the deficiency in our description to reside in the hadronic phase, cf. also ref. . Indeed, an artificial increase of the hadronic $`\mathrm{\Psi }^{}`$ dissociation cross-section by a factor of 5 clearly improves the agreement with the data. Further analysis is required by addressing, e.g., in-medium effects on $`D`$-mesons.

## 4 Predictions for RHIC and excitation function

At RHIC, due to much more abundant primordial charm-quark production and due to a stronger plasma suppression, the statistical $`J/\mathrm{\Psi }`$ production prevails over the direct component. This is apparent from Fig. 4 where the number of $`J/\mathrm{\Psi }`$’s per binary collision is plotted as a function of centrality. The full curve represents our prediction of ref. , with a minor adjustment of the input parameters from $`pp`$ collisions according to the (central values of the) recent PHENIX measurements  which reported $`\sigma _{J/\mathrm{\Psi }}^{pp}(\sqrt{s}=200\text{GeV})=3.8\pm 0.6\text{(stat)}\pm 1.3\text{(sys)}\mu `$b and $`\sigma _{c,\overline{c}}^{pp}(\sqrt{s}=200\text{GeV})650\mu `$b. Despite the large error bars, it appears that scenarios involving standard $`J/\mathrm{\Psi }`$ suppression only, as well as more extreme recombination scenarios, are disfavored by the data.

We also note that the probabilistic treatment of nuclear absorption, valid in the low energy regime, is not justified at RHIC energies. However, evaluations of nuclear absorption at collider energies using a more rigorous quantum-field theoretical approach  find that the suppression effects are quantitatively very close to the ones in the probabilistic treatment.

The transition between the regimes of predominantly direct (SPS) to statistical (RHIC) charmonium production can be further scrutinized experimentally at RHIC in terms of an excitation function for the ratio $`N_{J/\mathrm{\Psi }}/N_{c\overline{c}}`$, cf. Fig. 4. This ratio exhibits a nontrivial shallow minimum structure, a direct consequence of the interplay between suppressed and thermal production (assuming no anomalies in open-charm production).

## 5 Conclusions

We have shown that a two-component model for charmonium production in heavy-ion collisions, including both “direct” $`J/\mathrm{\Psi }`$’s arising from primordial $`N`$-$`N`$ collisions subject to nuclear, QGP and hadronic suppression, as well as “thermal” $`J/\mathrm{\Psi }`$’s emerging from coalescence of $`c`$ and $`\overline{c}`$ quarks at hadronization, accounts well for $`J/\mathrm{\Psi }`$ centrality dependencies measured at SPS. A potential discrepancy has been identified in the $`\mathrm{\Psi }^{}/\mathrm{\Psi }`$ ratio, which might hint at shortcomings in the evaluation of hadronic $`\mathrm{\Psi }^{}`$ dissociation. Extrapolating our approach to higher energies, the centrality dependence of $`J/\mathrm{\Psi }`$ production at RHIC is in line with preliminary PHENIX data. We furthermore emphasized the importance of an excitation function to determine a potential transition from (predominantly) direct to (mostly) statistical production, as indicated by an almost flat (but nontrivial) structure of the $`N_{J/\mathrm{\Psi }}/N_{c\overline{c}}`$ ratio as a function collision energy.


text:
# Dynamical Ordering of Driven Stripe Phases in Quenched Disorder

\[

## Abstract

We examine the dynamics and stripe formation in a system with competing short and long range interactions in the presence of both an applied dc drive and quenched disorder. Without disorder, the system forms stripes organized in a labyrinth state. We find that, when the disorder strength exceeds a critical value, an applied dc drive can induce a dynamical stripe ordering transition to a state that is more ordered than the originating undriven, unpinned pattern. We show that signatures in the structure factor and transport properties correspond to this dynamical reordering transition, and we present the dynamic phase diagram as a function of strengths of disorder and dc drive.

\]

Recently, there has been considerable interest in the dynamics of driven elastic media in the presence of quenched disorder and an applied drive. Physical systems that fall into this category include vortex lattices in disordered superconductors , sliding charge-density waves , driven electron crystals in the presence of random impurities , sliding friction , and domain wall motion . One of the most studied systems in this class is a vortex lattice driven by an applied current through disordered superconductors. Higgins and Bhattacharya  used transport measurements to map out a dynamical phase diagram for driven vortex matter based on transport signatures, and proposed the existence of three dynamical phases: a low drive pinned phase where the vortices do not move, a plastic phase where inhomogeneous flow and tearing of the highly disordered vortex lattice occurs, and an elastic flow regime where the lattice slides as a whole . Koshelev and Vinokur  investigated the driven vortex lattice system theoretically and numerically, observed three phases as a function of increasing applied drive, and found that the disordered lattice in the plastic flow regime can undergo a striking dynamical freezing transition in which the vortex lattice reorders at high drives.

Further theoretical work has shown that the recrystallized state is not fully ordered but is still strongly affected by transverse modes from the pinning. Thus the reordered state may form a moving smectic lattice with anisotropic ordering, where the vortices move in one-dimensional (1D) partially coupled channels aligned with the drive . This reordering transition to an aligned moving smectic state has been experimentally confirmed by transport  and direct imaging  experiments. In addition, numerical work has confirmed the presence of a field-driven plastic to ordered or elastic transition . Dynamical reordering has also been studied in sliding charge density waves  and driven Wigner crystals .

A natural question is whether these dynamical phases and reordering transitions can occur in other systems which do not form triangular lattices in the absence of quenched disorder. For example, many systems form “stripe” or “labyrinth” states , including diblock-copolymers, magnetic domain walls , flux in type-I superconductors, water-oil mixtures , and charge ordered or electron liquid crystal states in 2D electron gases  or superconductors . These stripe-like phases may be disordered or destroyed by quenched disorder. With an applied drive, however, it may be possible to return to a partially ordered state. One intriguing possibility is that the reordered state may be fully aligned with the drive, meaning that the moving ordered state would be more ordered than the stationary state observed without quenched disorder. This effect could be useful for straightening a labyrinth forming system when 1D arrays of aligned stripes are the desired pattern.

Here, we model a labyrinth or stripe forming system in 2D by conducting molecular dynamics simulations of interacting particles that have a long-range Coulomb repulsion and a short range exponential attraction. Such a system has previously been shown to produce stripe, bubble and crystalline phases depending on the particle density and the relative strength of the attractive interaction . We use overdamped dynamics which should be appropriate for systems such as colloids and magnetic domains. We do not take into account possible hydrodynamic effects. The equation of motion for a particle $`i`$ is $`𝐟=𝐟_{ij}+𝐟_p+𝐟_d=\eta 𝐯_i`$. Here, $`𝐯_i`$ is the particle velocity, $`\eta `$ is the damping term which we set to unity, and the force from the other particles is $`𝐟_{ij}=_jU(r)`$, where

$$U(r)=\frac{1}{r}Be^{\kappa r}.$$

For small $`r`$ the repulsive Coulomb term dominates. To avoid a divergence of the Coulomb term we cut off the force at very short distances ($`r0.1`$) which does not affect any of the results presented here. We also use a numerically efficient summation method for the long range Coulomb interaction . The force from the quenched disorder $`𝐟_p`$ comes from $`N_p`$ randomly placed attractive parabolic pins of maximum strength $`f_p`$ and radius $`r_p=0.3`$. The driving term $`𝐟_d`$ is applied in small increments of $`0.02`$ up to $`F_d=6.0`$. We measure the resulting particle velocities $`V_x=𝐯_i\widehat{𝐱}`$ by averaging the velocities at each $`f_d`$ increment value for 1000 time steps.

We fix $`B=1.25`$ and choose a particle density $`n_i=0.128`$ at which the equilibrium state in the absence of pinning is a labyrinth phase, rather than oriented stripes, which is consistent with the idea that stripe forming systems contain self-generated randomness . The initial state is prepared by annealing from a high temperature molten state to $`T=0`$. We conduct a series of simulations for varied $`f_p`$ and $`N_p`$, with the number of particles $`N_i=1280`$ and the system size $`L=100`$. We have also run simulations for fixed $`n_i`$ at various $`L`$ and $`N_i`$, and find that only the transient responses are affected.

We first consider the case for fewer pinning sites than particles, $`N_p=400`$. After annealing, we observe a distorted labyrinth pattern. Upon applying a drive, we find that there is no pinning threshold within our resolution and that the motion initially occurs by interstitial flow of unpinned particles around the particles that are trapped at the pinning sites. The overall velocity $`V_x`$ is still less than the free flow velocity since a portion of the particles remain pinned. As the drive increases, the system becomes increasingly disordered and signatures of the labyrinth ordering are further lost. The flow in this regime consists of river structures where the particles flow in individual channels, similar to those observed in vortex systems . Increasing the drive, we find that the system can reorder into two distinct moving states: (i) For weak disorder, the system reorders into a labyrinth state similar to the pattern that forms if the system is annealed and without quenched disorder. (ii) For sufficiently strong disorder, the particles reorder into a stripe state that is aligned in the direction of the drive. In Fig. 1 we show the real space images and the corresponding structure factor $`S(k)`$ for these phases.

In Fig. 1(a), the moving labyrinth phase shows stripe like features which have no long-range orientation. The structure factor has an inner ring of width corresponding to the spacing between the stripes. The ring shape reflects the fact that the stripes do not have a particular orientation. The outer ring in $`S(k)`$ corresponds to the length scale of the lattice of individual particles that forms within each stripe. This tendency to form a lattice is reflected by some remnant of six-fold order in the outer ring. We find the same $`S(k)`$ characteristics for both the moving and stationary labyrinth states. In Fig. 1(c), we show the real space image of the particles in the interstitial flow region. Here the system is much more strongly disordered than the labyrinth regime. There are areas where the particles form moving channels, one particle wide, aligned in the direction of drive. The structure factor (Fig. 1(d)) is much more diffuse, reflecting the stronger disordering. There are, however, now two peaks for small $`k`$, which is due the smectic ordering of the chain structures. In Fig. 1(e) we show the high drive state for strong disorder, where an ordered stripe state forms aligned in the direction of the drive. The corresponding structure factor in Fig. 1(f) shows two peaks for small $`k`$ corresponding to the smectic ordering of the stripes, in contrast to the ring structure observed in the labyrinth phase. Additionally the outer ring shows considerable six-fold modulations, due to the individual particles within the stripe forming a crystal structure.

The oriented stripe state for strong disorder is consistent with studies for vortex matter in 2D where, for strong disorder, the vortices form a smectic state with ordering occurring in only the transverse direction. For weaker disorder the vortex lattice anisotropy is reduced and the moving vortex lattice resembles the unpinned equilibrium state. The stripe ordering can also be viewed as a shear induced ordering. Shear forces are known to cause alignments into smectic phases . In our case there is no global shear, but there can be local shear if a portion of the particles remain pinned while other particles slide past. For weak disorder the local shear forces

are too weak to cause a realignment of the stripes.

We next show that the onset of the different phases strongly affects the transport properties. Additionally we show that a convenient order parameter is the fraction $`P_6`$ of six-fold coordinated particles. For a perfect triangular lattice $`P_6=1.0`$. We find that the stripe and labyrinth phases have distinct $`P_6`$ values near $`0.6\pm 0.015`$ and $`0.5\pm 0.3`$, respectively, while the disordered phase has $`P_6<0.5`$. In Fig. 2(a) we plot the velocity $`V_x`$ vs $`F_d`$ for a sample with strong disorder, $`f_p=6.5`$, along with the corresponding $`dV_x/dF_d`$ curve (solid line). In Fig. 2(b) we show the corresponding $`P_6`$ vs $`F_d`$. There is a prominent peak in $`dV_x/dF_d`$ near $`F_d=1.4`$ and a second prominent peak near $`F_d=2.5`$, each of which corresponds to a rise in $`P_6`$. For $`F_d<1.5`$ the system is disordered and only interstitial particles move through individual channels. For $`1.5<F_d<2.5`$, the system is still highly disordered, as seen from the low value of $`P_60.35`$. The flow now occurs throughout the sample, with individual particles intermittently escaping briefly from pins, and then being repinned. The first peak in $`dV_x/dF_d`$ corresponds to the increase in $`V_x`$ caused by the previously pinned particles depinning and taking part in the motion. The additional small peaks in $`dV_x/dF_d`$ correspond to specific portions of the pinned particles breaking free and becoming mobile in 1D chains. Similar small peaks in $`dV_x/dF`$ curves observed for driven vortices are believed to correspond to the depinning of large clumps of vortices . The peak near $`F_d=2.5`$ corresponds to the initial formation and alignment of the stripe state, which is also seen as a sharp rise in $`P_6`$ to around $`0.6`$. The stripe phase is fully formed for $`F_d=2.5`$ where $`P_6`$ plateaus. The peak in $`dV_x/dF`$ appears because the effectiveness of the pinning decreases as the stripes align. The $`dV_x/dF_d`$ curve

saturates to a constant value, reflecting free flow ohmic behavior. We find the same features for other simulations for $`2.5<f_p`$. For weaker pinning where a moving labyrinth forms, the $`dV_x/dF_d`$ curves shows only one peak \[shown in the inset of Fig. 2(b)\] corresponding to the initial formation of the moving labyrinth. For driven vortex lattices only one prominent peak is observed in the $`dV_x/dF`$ curves when the vortex lattice reorders and makes a transition from the plastic flow phase to the ohmic regime. We also find considerable hysteresis in the $`V_xF_d`$ curves in the oriented stripe phase, while there is no hysteresis for the system with weak disorder.

In Fig. 3 we show the dynamic phase diagram for $`F_d`$ vs $`f_p`$ for $`N_p=400`$. We delineate four phases: the interstitial phase, determined from the $`dV_x/dF_d`$ curves as the region before the first peak; the disordered flow region which is between the first peak in the $`dV_x/dF_d`$ curve and the plateau in $`P_6`$; the moving stripe phase which is identified from the real space images and from $`P_60.57`$; and the moving labyrinth phase, identified from the real space images and from $`P_60.5`$. Figure 3 shows that the oriented stripe phase appears only above a critical disorder strength. In Fig. 4(a) we show $`P_6`$ at a fixed high drive of $`F_d=5.0`$ for a series of simulations at the various $`f_p`$ used in creating Fig. 3. Fig. 4(a) clearly shows a jump near $`f_p=2.5`$ indicating the transition from the moving labyrinth to the moving stripe state, as a function of $`f_p`$.

We have also performed simulations at fixed $`f_p`$ and varied $`N_p`$, and find the same general phases. The moving stripe phase can be still realized with as few as $`N_p=50`$ pins. For $`N_p<50`$, the flowing state consists of coexisting aligned stripes and labyrinths. In Fig. 4(b) we show the $`V_xF_d`$ curves for varied $`N_p`$ of $`50`$ to $`2400`$. The onset of the different phases as a function of drive changes little

with increasing $`N_p`$. The value of $`V_x`$ in the interstitial region continuously decreases since more particles can be pinned. At very high pin density, the interstitial region is lost and is replaced with a pinned region. The corresponding phase diagram would have the interstitial flow regime replaced by a pinned flow regime. For the high drives $`F_d>2.6`$ all the curves become ohmic in the free flow regime. We have also simulated systems with weaker pinning $`f_p=1.0`$ and various $`N_p`$ and find that the high drive state again becomes a moving labyrinth phase.

In summary, we have investigated the dynamic phases in a driven system with quenched disorder that forms a disordered stripe or labyrinth phase in equilibrium without pinning. For weak pinning the system is disordered at low drives and reorders into a moving labyrinth phase at higher drives. Above a critical disorder strength, however, the system can reorder into a moving aligned stripe state. This stripe state is more ordered than the unpinned equilibrium state. Our results confirm that the dynamical reordering phenomenon studied in vortex lattices can be applied to other systems that do not form a triangular lattice in equilibrium. We have shown how the phases can be identified through transport properties and structure factors, and have mapped out the dynamical phase diagram as a function of disorder strength. Our results suggest that quenched disorder or other “obstacles” may be a useful way to engineer aligned domains in a labyrinth forming system. Additionally our results will be useful for the interpretation of transport studies in, e.g., driven electron-liquid crystal systems.

Acknowledgments This work was supported by DOE Contract No. W-7405-ENG-36.