text:
# Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 1footnote 11footnote 1Based 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.

## 1 Introduction

The Sombrero galaxy (M104 = NGC 4594) is well known for its dominating bulge/halo component which contributes with about 80% of the total luminosity (Kent, 1988; Baggett, Baggett & Anderson, 1998). This makes the Sombrero an important intermediary case between later-type spirals and elliptical/S0 galaxies, and an ideal place to study the relation between spiral bulges and early-type galaxies. One interesting aspect of such studies is to compare the globular cluster (GC) systems of early- and late-type galaxies. The Sombrero has a very rich GC system with an estimated total of $`2000`$ clusters (Harris et al., 1984; Bridges & Hanes, 1992; Forbes, Grillmair & Smith, 1997), about an order of magnitude more than the Milky Way and M31 (Harris, 1996; Fusi Pecci et al., 1993), but comparable to many early-type galaxies of the same luminosity. Using archive WFPC2 images, Larsen, Forbes & Brodie (2001, hereafter Paper I) detected a clearly bimodal color distribution for the Sombrero GCs, confirming hints by Forbes, Grillmair & Smith (1997). The number of red (metal-rich) and blue GCs were found to be roughly equal, in contrast to M31 and the Milky Way where the metal-poor GCs constitute about 2/3 – 3/4 of the total population, but again similar to the situation in early-type galaxies. Thus, we suggested in Paper I that the metal-rich GCs in the Sombrero are associated with the bulge rather than with its disk, in which case a much more modest number of metal-rich GCs would have been expected. The WFPC2 data presented in Paper I also showed that the mean colors of the two GC populations in the Sombrero are very similar to those in early-type galaxies as well as the Milky Way, M31 and other spiral galaxies (Forbes, Brodie & Larsen, 2001), with inferred mean metallicities of \[Fe/H\]$`1.40`$ and $`0.54`$.

Although there may not be a one-to-one correspondence between star formation and formation of globular clusters, information about the age distribution of GCs would presumably provide valuable insight about major star formation episodes in their host galaxies. An important, unresolved issue is whether age differences exist between the GC subpopulations in galaxies, and how large such age differences might be. This question is still controversial even in our Galaxy, where GC ages can be derived from color-magnitude diagrams (Stetson, Vandenberg, & Bolte, 1996; Buonanno et al., 1998; Rosenberg et al., 1999), so it is not surprising that it has proven to be very difficult to obtain reliable information about ages of extragalactic GCs from integrated photometry and/or spectroscopy. As shown by de Freitas Pacheco & Barbuy (1995), blue horizontal-branch morphology can have a significant impact on observed H$`\beta `$ indices, the most commonly used age indicator. More recently, Schiavon et al. (2002) have discussed a number of problems in the modeling of globular cluster spectra, and concluded that uncertainties in the luminosity function of the red and asymptotic giant branches, $`\alpha `$-element enhancement, horizontal branch morphology and other factors can affect spectroscopic age estimates by at least 3–4 Gyr. Observationally, useful constraints on typical GC ages require very accurate measurements of Balmer line indices, as the separation between e.g. the 10 Gyr and 11 Gyr isochrones is only about 0.1 Å in the equivalent width of H$`\beta `$. Furthermore, at lower metallicities ($`[\mathrm{Fe}/\mathrm{H}]1`$), recent models (Maraston 2001, priv. comm., Bruzual & Charlot 2001, priv. comm.) suggest that the Balmer lines no longer depend on age in a simple way for ages greater than about 10 Gyr.

A more indirect way to constrain the formation history of stellar populations in general, including GCs, is by studying their abundance patterns. The Milky Way bulge and halo populations, as well as many early-type galaxies, show an enhancement in the abundances of $`\alpha `$-elements relative to the Fe-group compared to the Sun (van den Bergh, 1995, and references therein). The enhancement is typically around $`[\alpha /\mathrm{Fe}]+0.3`$ (Carney, 1996) and is also seen in the Sombrero bulge stars (Peletier et al., 1999). The most commonly adopted explanation is that different elements are produced by stars with different lifetimes, of which supernovae of type Ia are generally assumed to have the longest lifetimes and make the strongest contribution to Fe. Thus, populations that formed on time-scales too short for SN Ia to have made any significant contribution to the abundance patterns show “enhancement” of other elements relative to Fe (e.g. Tinsley, 1979). Therefore, studies of abundance patterns may provide an alternative chronometer for the chemical evolutionary history of galaxies, although the timescale for development of SN Ia remains very uncertain.

Because of its relatively small discance and rich GC system, the Sombrero galaxy is an attractive target for detailed studies of individual GCs. With 8–10 m class telescopes, high-quality spectra can be obtained for many of the brighter clusters, and by co-adding spectra of several GCs it is possible to achieve sufficient S/N to bring the errors down to a level where model uncertainties rather than observational errors dominate. This makes it possible to test whether the models yield internally consistent results, e.g. by comparing age estimates based on different Balmer lines (H$`\beta `$, H$`\gamma `$, H$`\delta `$).

Spectroscopically, the GC system of the Sombrero has been studied by Bridges et al. (1997) who obtained rather low S/N spectra for 34 confirmed GCs. By co-adding the spectra they found a mean metallicity of \[Fe/H\] = $`0.70\pm 0.3`$, but the S/N was too low to constrain metallicities of individual clusters. In this paper we present new spectra for GCs in the Sombrero, selected from the Paper I sample. Of 17 GC candidates, 14 are confirmed as GCs. Combining our data with those of Bridges et al. (1997), the total number of confirmed GCs in the Sombrero is increased to 48. The best of our spectra have a S/N of about 50 per pixel, allowing us to accurately measure the metallicities of individual clusters and compare with previous photometric estimates. By co-adding the spectra of blue (metal-poor) and red (metal-rich) clusters, we also put constraints on the ages (Section 4.2) and abundance ratios (Section 4.3) of the Sombrero GC subpopulations. Finally, we use all available GC radial velocities in the Sombrero to constrain the mass of the Sombrero galaxy (Section 5).

## 2 Data

The data were obtained with the LRIS spectrograph (Oke et al., 1995) on the Keck I telescope during two observing runs on 2001, April 28 and May 19–20. We utilized the newly installed blue side on LRIS (“LRIS-B”) and a dichroic splitting at 5600Å to simultaneously obtain blue and red spectra, covering the regions $``$3500–5500Å and $``$5700–8500Å, respectively. On the blue side a 400 l/mm grism was used, providing a dispersion of 1.74Å pix-1 and a resolution of $`8`$Å. On the red side we used a 600 l/mm grating blazed at 5000Å, with a dispersion of 1.28Å pix-1 and a resolution of $`6`$Å. In addition to the multislit spectra of Sombrero GC candidates, we also obtained long-slit spectra for a number of radial velocity standard stars, flux standards, and Lick/IDS standard stars from Worthey et al. (1994).

Because of anticipated problems with “ghosts” in LRIS-B, we applied 2″ offsets to the telescope pointing for some exposures. Although these offsets were intended to be parallel to the slitlets, we found that they inevitably included some component perpendicular to the slits as well, resulting in greatly reduced signal from our science targets. Thus, we quickly decided to abandon this procedure. Another problem was related to the input coordinates. For the April run we used coordinates measured directly on the “central” and “inner halo” WFPC2 images from Paper I, using the METRIC task in the STSDAS package. Unfortunately, there were systematic differences between the two WFPC2 pointings which made it impossible to align all objects in the slitlets at the same time. For the May run we used direct images from the April run to alleviate this problem, but as a result of these slitmask alignment problems, we achieved a significantly lower S/N for some of the objects than could have been possible under optimal circumstances.

Following initial processing of the images (bias subtraction, flatfielding), the spectra were extracted using standard tools in the IRAF SPECRED package222IRAF is distributed by the National Optical Astronomical Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc. under contract with the National Science Foundation. As the spectra were obtained over several hours, significant shifts in the wavelength scale due to mechanical flexure were found to be present between individual integrations, amounting to up to $`5`$ Å. Assuming that all spectra in a given exposure were subject to the same shifts, the shifts were determined by cross-correlating the spectra of the brightest objects in each exposure with the spectra in a reference exposure, taken immediately before the arc spectra. For the red spectra, the wavelength zero-point was further fine-tuned using the NaD night sky lines at $`\lambda \lambda `$ 5890 and 5896 Å. Finally, the individual spectra of each object were co-added, using a sigma-clipping algorithm to reject cosmic ray events. The S/N at each pixel in the combined spectra was estimated from the variance of the combined pixels. The total integration time was $`5\frac{1}{2}`$ hours.

Table 1 lists all science objects for which spectra were obtained. The object IDs refer to Paper I, where the coordinates are listed. For convenience, we have also included the $`V,I`$ photometry from Paper I. Note that the brightest cluster, H2–22, was saturated on the WFPC2 data used in Paper I, making its $`VI`$ color unreliable. The radial velocities in Table 1 are generally measured on the red spectra, for which the wavelength calibration zero-point could be checked using the night sky lines. The errors on the radial velocities in Table 1 are random errors, based on the uncertainties in the location of the cross-correlation peak reported by the FXCOR task. These do not take into account various systematic errors, mostly due to uncertainties in wavelength calibration. We estimate that such errors are on the order $`0.5`$Å, or $`25`$ km/s. No skylines were present in the blue spectra, and differences of 100–150 km/s between radial velocities measured on the blue and red spectra were not uncommon. Another problem with the wavelength calibration in the blue stems from the fact that the calibration lamps in LRIS only had a few (5–7) suitable lines. The best fitting 2nd order polynomial wavelength solutions in the blue generally had a significant residual scatter (1–2Å), presumably indicating that higher-order terms ought to be included in the wavelength solution, but the small number of calibration lines made it impossible to include higher-order terms.

The radial velocity of the Sombrero galaxy is listed as 1024 km/s in the RC3 catalogue (de Vaucouleurs et al., 1991). If we exclude the three objects in Table 1 with RV$`<`$600 km/s, we are left with a sample of 14 globular clusters. These have a mean radial velocity of 1019 km/s and a dispersion of 246 km/s, in excellent agreement with the radial velocity of the Sombrero itself. Of the three objects with radial velocities less than 600 km/s, two (H2-06 and H2-10) are very compact and essentially unresolved on the HST images with effective radii $``$ 0.3 pc (Larsen, Forbes & Brodie, 2001) and are most likely foreground stars. The third (H2-27) is quite well-resolved on the HST images and has an effective radius of 3.4 pc, typical for globular clusters. It therefore seems likely that this object is indeed a globular cluster. The S/N of this spectrum is very poor and the cross-correlation signal was not very strong, so the radial velocity may be in error. In fact, using the blue spectrum of this object, we get a radial velocity of 1069 km/s. Because of the uncertain radial velocity and poor S/N we exclude this object from further analysis in this paper.

Before spectrophotometric indices were measured, the spectra must be corrected for radial velocity. Rather than blindly applying the radial velocities in Table 1 to the blue spectra, we used a number of easily identifiable features such as the H$`\beta `$, H$`\gamma `$, H$`\delta `$, Ca II H$`+`$K lines and the G-band to determine small zero-point shifts and scale corrections to the dispersion solutions. Thus, the wavelength scale of the final spectra used for the spectrophotometric analysis was typically accurate to better than 1Å in the region of interest (i.e. $`3700`$–5400 Å).

## 3 Measuring Lick/IDS indices on LRIS-B spectra

In order to test if there might be any systematic differences between indices measured on LRIS-B spectra and the standard Lick/IDS system, we observed a number of stars from Worthey et al. (1994) in longslit mode. In Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. we compare our LRIS-B measurements of all Lick/IDS indices that were covered by our blue spectra with the standard values. Generally, our measurements agree well with the standard values, although there are slight systematic offsets for some of the narrower indices like Ca4227 and Ca4455. This may be attributed to the slightly higher resolution of LRIS-B compared to the Lick/IDS spectrograph.

Unfortunately, the wavelength range covered in longslit mode by LRIS is slightly different from the one covered in multislit mode and does not extend beyond 5300Å. Therefore, several key indices (such as Mg2, Fe5270 and Fe5335) redward of this limit could not be compared. Furthermore, since this part of the spectrum is not covered by the flux standards, the flux calibration of the multislit spectra beyond 5300Å also becomes uncertain. This is particularly problematic because the blue reflection efficiency of the dichroic already begins to drop around $``$5300Å.

## 4 Constraining ages and metallicities

### 4.1 Metallicities

We derived metallicities for each GC using the calibration by Brodie and Huchra (1990), but adopting the modification to the error on the combined metallicity estimate described in Larsen and Brodie (2000). Metallicities for each of the individual calibrators, as well as the combined metallicity estimates are listed in Table 2. For comparison we have also included photometric metallicities based on the $`VI`$ colors in Table 1 and the calibration in Kissler-Patig et al. (1998). The spectroscopic and photometric metallicity estimates are compared in Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. and generally agree quite well. Three representative globular cluster spectra with low, intermediate and high metallicities are shown in Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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..

### 4.2 Comparison with population synthesis models

High S/N spectra have the potential to remove the age-metallicity degeneracy inherent in photometric studies. At least in principle, ages are well constrained by Balmer line strengths, while metallicity information is provided by indices such as Fe5270, Fe5335 and Mg2 (Worthey, 1994).

Here we compare our Sombrero spectra with two sets of recent population synthesis models, kindly provided by C. Maraston (see Maraston & Thomas, 2000) and R. Schiavon (Schiavon et al., 2002b). Like most of the current models, the Maraston models are based on the “fitting” functions by Worthey et al. (1994). The Schiavon models are based on a new set of fitting functions derived from the spectral library of Jones (1999), and also include empirical corrections to the luminosity function of red giants (Schiavon et al., 2002). We use a preliminary version of the Schiavon models where only the Balmer line indices are based on the new fitting functions, i.e. other indices such as Mg$`b`$, Mg2 and the various Fe indices are still based on the Worthey fitting functions. The Schiavon models also cover a more limited age/metallicity range than the Maraston models, but the comparison nevertheless provides some useful insight into the differences between various available models. As discussed in Sect 3, our instrumental system appears to be very similar to the Lick/IDS system so we have not applied any corrections to spectrophotometric indices measured on the LRIS-B spectra.

As mentioned earlier, some indices may be compromised by difficulties with the calibration of the blue spectra redward of 5300Å. The Mg2 index is particularly vulnerable to such problems because of the wide separation between the continuum passbands, and offsets are frequently found between the standard Lick/IDS indices and instrumental systems (e.g. Kuntschner, 2000). In Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. we compare our measurements of Mg2 with Mg$`b`$, which is slightly narrower but measures essentially the same feature. However, while the continuum passbands of Mg2 are located about 150Å from the feature passband, those of Mg$`b`$ are located adjacent to it, making Mg$`b`$ much less sensitive to errors in the continuum slope. Generally, Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. shows that the two Mg sensitive indices correlate well, but the slope of our observed Mg$`b`$Mg2 relation is shallower than the Maraston model predictions, which in turn agree well with Mg$`b`$Mg2 plots in e.g. Trager et al. (1998) and Beasley et al. (2000). One possible explanation for this discrepancy is that errors in the continuum level at the position of the red continuum passband cause us to measure too low values for the Mg2 index because of uncorrected reduced throughput by the dichroic. Another possibility is that peculiarities in the abundance patterns affect the continuum passbands of Mg$`b`$ and Mg2 in different ways (Tripicco & Bell, 1995).

In Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. we compare various index measurements with the Maraston and Schiavon model predictions. The Schiavon models are tabulated at 5, 7.9, 11.2 and 14.1 Gyr and at metallicities of $`[\mathrm{Fe}/\mathrm{H}]`$ = $`0.7`$, $`0.4`$, 0.0 and $`+0.2`$. While the Maraston models are tabulated at 1 Gyr intervals from 2 to 15 Gyrs, we have plotted only the 5, 8, 11 and 14 Gyr models, for easy comparison with the Schiavon models. Note that the metallicity divisions of the Maraston models are different from those of the Schiavon models ($`[\mathrm{Fe}/\mathrm{H}]`$= $`2.25`$, $`1.35`$, $`0.33`$, 0.00 and $`+0.35`$). Panels (a) and (b) show H$`\beta `$ vs. $`\mathrm{Fe}`$ and the Mg$`b`$ index, where $`\mathrm{Fe}`$= (Fe5270 + Fe5335)/2. The Maraston and Schiavon model grids are shown with solid lines and dashed lines, respectively. In panels (c) and (d) we show H$`\gamma _A`$ and H$`\delta _A`$ vs. $`\mathrm{Fe}`$, using the definitions in Worthey & Ottaviani (1997). Because of the difficulties with the calibration of Mg2 discussed earlier, we have chosen Mg$`b`$ instead of Mg2 even though the formal errors on Mg$`b`$ are somewhat larger. The errors on the indices were estimated directly from the S/N of the spectra.

To first order, Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. indicates high ages for most clusters, although the error bars are generally still too large to put strong constraints on the ages of individual clusters. In order to compare the mean properties of metal-rich and metal-poor GCs, we therefore co-added all spectra with \[Fe/H\]$`<1`$ and \[Fe/H\]$`>1`$, respectively. The co-added spectra are shown in Figure Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. and the corresponding datapoints are shown with filled circles in Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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.. In the H$`\beta `$ plots both the metal-poor and metal-rich co-added datapoints fall near the $`11`$ Gyr Maraston isochrones, with little evidence of any significant age differences between the two within the formal $`2`$ Gyr error bars on the co-added data. The Schiavon models tend to predict lower ages for a given H$`\beta `$ index than the Maraston models, especially at younger ages and lower metallicities. At the location of the metal-rich Sombrero GCs in the $`\mathrm{Fe}`$– H$`\beta `$ plane the difference is relatively modest and amounts to an age difference of only 1–2 Gyrs, but at lower metallicities the difference between the two sets of models becomes more pronounced and amounts to about 4 Gyrs at the lowest metallicity tabulated by Schiavon. Thus, the trend in these models is to decrease the ages of the metal-poor clusters more than those of the metal-rich ones, relative to the Maraston models.

The H$`\gamma _A`$ and H$`\delta _A`$ measurements basically confirm the ages estimated from H$`\beta `$, but note that both H$`\gamma _A`$ and H$`\delta _A`$ reach a minimum around 11 Gyr at low metallicities and then increase strongly at higher ages. At $`[\mathrm{Fe}/\mathrm{H}]=1.35`$, the Maraston models predict about the same H$`\gamma _A`$ and H$`\delta _A`$ strengths at 15 and 5 Gyrs. The Schiavon models do not include tabulations of H$`\gamma _A`$ and H$`\delta _A`$, but in Figure Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. we show their predictions for another H$`\gamma `$ index defined by Vazdekis et al. (2001). This H$`\gamma `$ index (in the following denoted H$`\gamma _{\mathrm{vaz}}`$) differs from H$`\gamma _A`$ in that the continuum passbands are located closer to the feature passband, thereby avoiding the strong metallicity dependence of H$`\gamma _A`$ which results from its blue continuum passband coinciding with the G–band. In fact, Figure Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. shows that H$`\gamma _{\mathrm{vaz}}`$ is practically independent of metallicity, at least within the range spanned by the models. Taken at face value, Figure Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. would suggest that the metal-rich clusters are slightly *older* than the metal-poor ones, but both sub-populations are compatible with the same (high) age within the errors. Thus, although it should be stressed that even relative ages may still be somewhat model-dependent, it seems reasonable to conclude that both GC sub-populations in the Sombrero are old, and probably coeval to within the measurement uncertainties, i.e. with age differences of 2–3 Gyrs at most.

### 4.3 Abundance ratios

Some attempts to model the effect of $`\alpha `$-element enhancement on Lick indices have been made by Milone, Barbuy, & Schiavon (2000), who modeled the effect of $`\alpha `$-element enhancement on the Lick Mg2 index and a broad index (TiO12.5) covering the TiO band centered at 7257 Å. The TiO12.5 index is included in our red LRIS spectra. Because of the problems with the calibration of Mg2, we converted the Milone et al. model predictions for Mg2 to Mg$`b`$, using the Maraston models. We found that the linear relation $`\text{Mg}b=14.03\text{Mg}\text{2}+\mathrm{\hspace{0.17em}0.361}`$ approximated the Maraston models to within 0.1 Å for the age range 10–15 Gyr and $`\text{Mg}\text{2}<0.3`$ (Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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.).

In addition to the strong terrestrial O2 absorption bands, which are located just outside the TiO12.5 continuum passbands, the wavelength region around 7200–7500Å is affected by terrestrial H2O absorption. We corrected our red spectra for this by normalizing the flux standard spectra with a smooth continuum and then dividing the GC spectra by the normalized flux standard spectra. This correction affected the TiO12.5 index downwards by about 7Å. Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. shows our measurements of TiO12.5 and Mg$`b`$ on the co-added spectra, compared to the alpha-enhanced models by Milone, Barbuy, & Schiavon (2000). In the $`\mathrm{Fe}\text{Mg}b`$ plot we also show the Maraston models for an age of 12 Gyr. Note that the most metal-poor Maraston models are actually rather similar to the $`[\alpha /\mathrm{Fe}]=+0.3`$ models by Milone, Barbuy, & Schiavon (2000), but approach the $`[\alpha /\mathrm{Fe}]=0`$ models at higher metallicities. This is probably because many of the metal-poor stars used to set up the Worthey fitting functions are alpha-enhanced (Maraston & Thomas, 2000). For the metal-rich Sombrero GCs, both Mg$`b`$ and TiO12.5 indicate clearly supersolar $`[\alpha /\mathrm{Fe}]`$ ratios around $`[\alpha /\mathrm{Fe}]+0.4`$. For the metal-poor clusters, the TiO12.5 index again indicates supersolar $`[\alpha /\mathrm{Fe}]`$ ratios similar to those of the metal-rich clusters, although the Mg$`b`$ plot actually suggests $`[\alpha /\mathrm{Fe}]`$ closer to 0 for the metal-poor clusters. However, at low metallicities the separation between the different $`[\alpha /\mathrm{Fe}]`$ models is smaller and it seems most likely that both the metal-poor and metal-rich GCs in the Sombrero have super-solar $`[\alpha /\mathrm{Fe}]`$ ratios similar to those observed in Milky Way GCs, the bulge of the Sombrero and in early-type galaxies.

## 5 Mass of the Sombrero

Globular clusters can be used as tracers of the mass of their host galaxy, as demonstrated e.g. by Huchra & Brodie (1987) (for M87) and Kissler-Patig et al. (1998) (NGC 1399). We combined our data in Table 1 with data for 34 Sombrero GCs from Bridges et al. (1997). Our sample does not overlap with that of Bridges et al. (1997), so by combining theirs and our data we get a total of 48 GCs.

Fig. Keck Spectroscopy of Globular Clusters in the Sombrero Galaxy 111Based 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. shows the radial velocities as a function of the offset in R.A. from the center of the Sombrero (the major axis is aligned roughly east-west). Bridges et al. (1997) found hints of rotation in the cluster system, but our data clearly do not add much to the significance of their result. For completeness we also plot the radial velocities for the metal-rich and metal-poor clusters from our sample separately, but the data are clearly insufficient to provide any constraints on rotation in either sample.

Following Huchra & Brodie (1987), we can use the GC radial velocities to obtain a mass estimate for the Sombrero in two ways. The virial mass is estimated as

$$M_{\mathrm{VT}}=\frac{3\pi N}{2G}\frac{\mathrm{\Sigma }_i^NV_i^2}{\mathrm{\Sigma }_{i<j}1/r_{ij}}$$
(1)
and the projected mass is

$$M_P=\frac{f_p}{GN}\left(\mathrm{\Sigma }_i^NV_i^2r_i\right)$$
(2)
In these equations, $`r_{ij}`$ is the distance between the $`i`$th and $`j`$th clusters, $`r_i`$ is the distance between the $`i`$th cluster and the galaxy center and $`V_i`$ is the velocity difference between the $`i`$th cluster and the mean system velocity. $`N`$ is the number of clusters. The value of $`f_p`$ depends on the cluster orbits; following Bridges et al. (1997) we adopt $`f_p=16/\pi `$ for isotropic orbits.

For the total sample of 48 clusters we get $`M_{\mathrm{VT}}=(8.1\pm 1.8)\times 10^{11}\text{M}\text{}`$ and $`M_P=(5.3\pm 1.0)\times 10^{11}\text{M}\text{}`$ within 17 kpc, adopting a distance of 8.7 Mpc (Larsen, Forbes & Brodie, 2001). Using only our 14 clusters, we get $`M_{\mathrm{VT}}=(3.2\pm 1.5)\times 10^{11}\text{M}\text{}`$ and $`M_P=(4.2\pm 2.1)\times 10^{11}\text{M}\text{}`$ (within 13 kpc). For comparison, Bridges et al. (1997) got $`M_P=\mathrm{\hspace{0.17em}5.2}_{3.9}^{6.7}\times 10^{11}\text{M}\text{}`$. The 12 clusters in the central HST pointing are useful as probes of the mass located within the central regions of the galaxy, and yield $`M_{\mathrm{VT}}=(2.1\pm 1.1)\times 10^{11}\text{M}\text{}`$ and $`M_P=(2.0\pm 1.3)\times 10^{11}\text{M}\text{}`$ within the central 4.5 kpc. The errors were estimated by a standard “jackknife” procedure.

## 6 Summary and conclusions

We have analyzed Keck / LRIS-B spectra for 17 globular cluster (GC) candidates in the Sombrero galaxy, of which 14 are confirmed as true GCs based on their radial velocities. One additional object has too low S/N for reliable radial velocity measurements, but based on the HST images presented by Larsen, Forbes & Brodie (2001) it is probably also a GC. We find generally good agreement between spectroscopic metallicities and previous photometric metallicites based on $`VI`$ colors. Comparison of Lick/IDS indices with population synthesis models indicates that both the metal-rich and metal-poor clusters are 10–15 Gyr old, with preferred ages around 11-12 Gyr. We see no detectable age differences between the two GC populations, although differences of a few Gyrs cannot be ruled out. Comparison of $`\mathrm{Fe}`$, Mg$`b`$ and TiO indicates that the Sombrero GCs probably have super-solar $`[\alpha /\mathrm{Fe}]`$ ratios $`+0.4`$, similar to or perhaps slightly larger than those observed in Milky Way globular clusters.

Overall, the Sombrero GC system shares many similarities with that of the Milky Way, as well as many early-type galaxies. With few exceptions, both GC populations generally appear to be old, and another common feature seems to be an enhancement in $`\alpha `$-element abundances, probably indicating that the gas out of which the clusters formed was enriched predominantly by type-II supernovae. This suggests that the majority of GCs formed early on, before type Ia SNe made a significant contribution to chemical enrichment, and presumably within a few Gyrs after star formation commenced. The metal-rich GCs in the Sombrero, like those in most galaxies, probably formed early on when galaxies were still assembling with high star formation rates and gas densities providing favorable conditions for the formation of large numbers of massive star clusters. The metal-poor clusters might have formed even earlier at high redshifts, triggered by ionization fronts (Cen, 2001).

Finally, we combine our sample with data for 34 GCs from Bridges et al. (1997) and obtain virial and projected masses for the Sombrero of $`M_{\mathrm{VT}}=(8.1\pm 1.8)\times 10^{11}\text{M}\text{}`$ and $`M_P=(5.3\pm 1.0)\times 10^{11}\text{M}\text{}`$, respectively. This is in good agreement with previous estimates of the mass of the Sombrero. For 12 clusters located within the central 4.5 kpc (projected) we estimate $`M_{\mathrm{VT}}=(2.1\pm 1.1)\times 10^{11}\text{M}\text{}`$ and $`M_P=(2.0\pm 1.3)\times 10^{11}\text{M}\text{}`$.

This work was supported by National Science Foundation grant number AST9900732. We are grateful to Dr. Andre Milone for providing his $`\alpha `$-enhanced models in electronic form and to Dr. Ricardo Schiavon for allowing us to use his population models prior to publication. Helpful comments by an anonymous referee are appreciated. MB thanks the Royal Society of London for their support.


format: text

text:
# Nearby quasar remnants and ultra-high energy cosmic rays

## 1 Introduction

Different experiments over the past few decades have detected several giant air showers, confirming the arrival of cosmic rays (CRs) with energies up to a few hundred EeV (1 EeV $`10^{18}`$ eV). The nature and origin of these energetic particles remain a mystery . CRs, however, can not travel unaffected through inter-galactic space. The thermal photon background becomes highly blueshifted for ultra-relativistic protons, and the reaction $`p\gamma \mathrm{\Delta }^+\pi ^0p`$, and similar others, effectively degrade the primary proton energy. This provides a strong constraint on the proximity of CR-sources, discovered early on by Greisen, and by Zatsepin & Kuz’min, and referred to as the GZK cutoff . Specifically, fewer than 20% of 300 EeV (100 EeV) protons can survive a trip of 18 (60) Mpc. For nuclei and photons the situation is, in general, more drastic .

The possible solutions to this puzzle appear to fall into three broad categories, viz: 1) There are many nearby sources (e.g., based on a substantial present-epoch population of galaxies hosting core supermassive black holes). 2) There are only a few nearby sources (e.g. ) and particles are isotropized by strong deflections in Galactic and/or extragalactic magnetic fields of micro-Gauss strength, close to existing upper limits (e.g., see , also ). 3) The particles somehow travel relatively unhindered through the Cosmic Microwave Background radiation field, either by virtue of being some exotic new kind of weakly photon-interacting entity, or by violating the Lorentz symmetry of special relativity (e.g. ), or by an as yet unknown other effect.

The arrival directions of the primary particles could be a useful source of information about the origin(s) of CRs. However, due to the small number of events, particularly at the highest energies, correlation studies are to be considered preliminary. For instance, the first five events observed with $`E>80`$ EeV did in fact point toward high redshift radio-loud quasars, astrophysical environments that could well accelerate CRs above the GZK energies via shock mechanisms (see Farrar and Biermann, in Ref. ). However, with the inclusion of subsequent data, this association now seems to have disappeared .

In this communication we explore whether UHECRs can be in any way related to present-epoch supermassive black holes. To that end, we analyze what we can learn from the arrival directions of cosmic rays concerning their possible correlation with nearby massive dark objects, candidate QRs (quasar remnants), and their underlying astrophysics.111The term “quasar remnants” was introduced by Chokshi and Turner  to describe the present-epoch population of dead quasars harboring supermassive black hole nuclei. The compact dynamo model has been proposed as a natural mechanism for accelerating cosmic rays in such environments . In this model UHECRs are produced in nearby dead quasars harboring spinning supermassive black holes. The required emf is generated by the black hole induced rotation of externally supplied magnetic field lines threading the horizon. The observed flux of CRs would apparently drain only a negligible amount of energy from the black hole dynamo, and particles up to at least 100 EeV are expected. It is then interesting to ask if we are able to see any correlation between the sky position of these QRs and those of the highest energy CRs, assuming different simple configurations for the intervening magnetic field.

## 2 Quasar remnant candidates

The Nearby Optical Galaxy (NOG) catalog of Giuricin et al.  is a complete magnitude-limited (corrected blue total magnitude $`B14`$), distance-limited (redshift $`z0.02`$) sample of several thousand galaxies of latitude $`\left|b\right|>20^{}`$, with their morphology $`T`$ also provided . We shall impose very restrictive selection criteria, those necessary for obtaining candidate objects providing the most favorable setting for a black hole based compact dynamo model of UHECR production . The key goal of our strategy is to select a priori an NOG subsample of galaxies that are likely to be quasar remnants. Therefore, we are seeking optically bright galaxies whose luminosities are bulge dominated (e.g., giant ellipticals). Imposing a GZK-related horizon of $`50`$ Mpc assures that there are no quasars in our sample. We already know that the best determined (and most massive) black hole nuclei tend to be associated with bulge luminosities corresponding to $`M_B21`$ and brighter (for tables, see Kormendy’s website at http://chandra.as.utexas.edu/$``$kormendy/wwwbhtable-tech or the second reference in ).

Firstly, then, we shall impose a cutoff in redshift, requiring all galaxies in our sample to be within $`z0.01`$, which for $`H_0=75`$ km/s/Mpc corresponds to 40 Mpc. Since the CRs detected are mostly in the northern hemisphere (except those coming from SUGAR, which are not used in the present analysis), we shall impose the restriction that all QRs have equatorial latitudes north of $`10`$ deg. We further require the absolute blue magnitude to be brighter than $`M_B=21`$ (for $`H_0=75`$ km/s/Mpc) and an RC3 (Third Reference Catalog of Bright Galaxies) morphological type $`T<3`$. A more negative $`T`$ indicates greater bulge prominence. In order to correlate the UHECR arrival directions with their putative QR origins, we require that these candidate QRs not lie within rich clusters (those having more than 50 members). Some of the most massive QRs could reside in rich clusters of galaxies (e.g., all four QRs considered by Boldt and Loewenstein in Ref.  are in rich clusters), but the magnetic field strength in those clusters is presumably several micro-Gauss ; hence, UHECRs originating in such QRs would be extremely deflected away from their sources. In addition, there are only a handful of galaxies within large clusters in the Giuricin et al.’s sample that survive all other constraints in order to be declared plausible quasar remnant candidates. Using these selection criteria, we obtain a sample of 12 candidate QRs. These candidate QRs are listed in Table 1, and their astrophysical properties are analyzed below.

The UHECR sample used is that obtained with AGASA above 40 EeV . There are 38 such events at $`\left|b\right|>20^{}`$; the angular precision of their arrival directions was estimated as a circle of radius 1.6 degrees. For energies $`>100`$ EeV we also consider the 7 events observed at $`\left|b\right|>20^{}`$ compiled by Sigl et al. , 4 of which were obtained with AGASA.222See Table 1 in  for this compilation and the associated error estimates. We note that the statistical test we shall report was blind, i.e. we did not know beforehand if any of the QR galaxies were coincident with high energy cosmic rays. The simulation technique is the same as that used in Refs. .

Looking for superposed QRs in the cosmic ray error circles, we find the results reported in Table 2. When we consider the AGASA sample with $`E>40`$ EeV, there is an excess in the real result which approaches $`3\sigma `$. This could argue in favor of a correlation between our sample of QRs and UHECRs at small angular scales. This is not the case for the CRs with energies above 100 EeV, for which we found only an apparent excess for the nominal error circles, not significant enough to give to it any confidence. A summary of these results is given in Figure 1. There we show the number of real coincidences in rings surrounding the position of the central QR as a function of the internal radius of the ring, $`A`$. Each ring has an angular width of 1.6 deg, and the number of coincidences is shown with angular offsets between $`A`$ and $`A+1.6`$. An excess of coincidences is clearly evident for the innermost radii. At larger radii, the real superpositions are compatible with the random coincidences. This result stands disregarding the value chosen for the size of the ring.

Note that when increasing the error circles, no additional QRs are involved in the coincidences, but instead, there are new CRs involved. This is what makes the real result increase. Details of the coinciding pairs are given in Table 3. It is risky, of course, to evaluate the level of significance that we should attribute to this correlation, especially when so few events/objects constitute the samples involved. The 12 galaxies, however, are selected a priori on a restrictive physical basis, aimed at identifying QRs harboring supermassive black holes. Furthermore, we can see that the correlated pairs identify an apparently preferred region of the sky, as Figure 2 shows. The expected number of random coincidences decreases from 0.40 to 0.14 when the region of the sky is restricted to an area given roughly by $`160<\alpha <200`$ and $`50<\delta <80`$. This is $`5\sigma `$ below the real result of 2 coincidences. This latter result, however, is obtained a posteriori, and its importance is thereby reduced.

If we consider that charged particles (protons) are accelerated in these QRs and then travel through the intergalactic and Galactic magnetic fields towards the Earth, we expect their trajectories to be deflected. When the Larmor radius of a particle ($`r_\mathrm{L}10^2`$ Mpc $`E_{20}/B_9`$) is much larger than the coherence length of the magnetic field $`\mathrm{}_{\mathrm{coh}}`$, the characteristic deflection angle $`\theta `$ from the direction of the source, located at a distance $`D`$, can be estimated assuming that the particle makes a random walk in the magnetic field (see e.g. )

$$\theta (E)3.8^{}\left(\frac{D}{50\mathrm{Mpc}}\right)^{1/2}\left(\frac{\mathrm{}_{\mathrm{coh}}}{1\mathrm{Mpc}}\right)^{1/2}\left(\frac{B_9}{E_{20}}\right),$$
(1)
where $`E_{20}`$ is the energy of the particle in units of $`10^{20}`$eV, and $`B_9`$ is the magnetic field in units of $`10^9`$ G. We can see that scattering in large scale magnetic irregularities $`𝒪`$ (nG)  are enough to bend the orbits of trans-GZK protons by about $`4`$ deg in a 50 Mpc traversal. As exhibited in Table 3, the CR angular offsets observed for these QRs are much smaller than $`\theta `$. For a variety of assumed magnetic field scenarios, $`\theta `$ is often substantially larger than the estimated AGASA measurement error of 1.6 degrees. If, in fact, due to the value of the inter-galactic magnetic field, such is the case, it would then indicate that

1. the number of apparent QR/UHECR associations is no longer above random expectations (see Figure 1 for large value of $`A`$), and
2. the a priori probability is relatively very small for offset angles as little as those actually observed (Table 3).

Should the apparent clustering of correlated pairs be supported by future data, what are the viable scenarios under which this could occur? One possibility is to consider that the intergalactic medium between Earth and the three apparently ‘clustered’ QRs is sufficiently different from the intergalactic medium in front of the remaining nine objects that are much more uniformly distributed on the accessible sky. For the deflection of an energetic (60 EeV, the mean of the two CR energies in Table 3) proton in traversing 34 Mpc (the mean of the 3 QR distances in Table 3) to be less than a degree, we would need $`B<2\times 10^{10}\mathrm{}_{\mathrm{coh}}^{0.5}`$ G, which appears to be not as drastic a difference from the canonical nano-Gauss $`B`$ field and coherence length of 1 Mpc that are usually assumed. Indeed, the IGM B field is likely to be an order of magnitude less than a nano-Gauss in voids comparable in size to the GZK horizon (Peter Biermann; 2002, personal communication). Independent estimates of the Intergalactic Medium (IGM) magnetic field towards these galaxies would prove very useful in assessing this explanation.

It is important to note that, in some directions, the magnetic field of our own galaxy could well lead to a deflection of up to several degrees for primaries with energies below 60 EeV; see, for instance, Table 1 of Ref. . A recent study of this issue was presented by Alvarez-Muñiz et al. . However, it is important to realize that the possible filamentary topology of the Galaxy’s magnetic field (Gerrit Verschuur; 2002, personal communication) would likely allow some directional windows, albeit narrow, where the deflection of an UHECR could be much less than typical. Deflections due to the magnetic fields would of course be avoided if the primary were a photon, generated in the neighborhood of the QR via an accelerated charged particle interaction. A photon primary is ruled out for the highest cosmic ray energy event detected with Fly’s Eye , but there is yet some discussion of the likelihood of having some photons as primaries in the AGASA sample . Protheroe and Johnson (see second reference in ) have studied the propagation of ultra-high energy gamma rays in the intergalactic magnetic field in detail. They have shown that if the intergalactic magnetic field is low enough, say $`B10^{10}`$ G, a photon can survive a distance comparable to those of the QRs in our sample and be a plausible primary for the observed showers. A new neutral hadron could obviously be another extreme possibility. However, as we have seen, the involved QRs are all within a relatively small region of the sky; their angular separation is clearly less than the average angular separation between neighbors in the remaining sample of nine. Why would sources of neutral hadronic CRs be restricted to a limited angular region of the sky? Albeit suggestive, we warn the reader that the effect observed here might be no more than an artifact resulting from the small number statistics available at this time. It is, however, interesting to study the underlying astrophysics of our candidates as possible cosmic rays emitters.

## 3 Underlying Astrophysics

Table 4 exhibits the physical attributes of our sample of QR candidates, viz: morphological type, mass of supermassive black hole nucleus ($`M_{\mathrm{BH}}`$), galactic bulge mass ($`M_{\mathrm{GB}}`$), black hole accretion rate $`\dot{m}`$ (in Eddington units), the maximum resulting proton energy ($`E_{\mathrm{max}}`$) expected, and lower limits to the proton’s radiation length ($`\mathrm{\Lambda }`$) for losses arising from photo-pion production in ultra-relativistic collisions with the ambient electromagnetic radiation (photon) field. That the nuclei of these galaxies are indeed dark as well as massive may be appreciated by noting that the radio luminosity ($`\nu L_\nu `$ at 6 cm) within a 5 arcsec beam  is in no instance greater than $`10^8`$ (in Eddington units); the ROSAT obtained X-ray luminosity , within a broader beam that extends beyond the nucleus, is in every case already less than $`10^5`$ (in Eddington units). The four QRs that appear to be correlated with UHECRs ($`>40`$ EeV) are identified here by their offset angles relative to the measured arrival direction of the primary particle initiating the associated air shower; this offset is never greater than the estimated experimental error in determining this arrival direction. For the worst case, NGC 4589, the associated Haverah Park UHECR direction measured has an angular uncertainty radius of 5.6 degrees.

Along with the six published values cited in Table 4, we estimate the black hole mass ($`M_{BH}`$) for all twelve QRs from the velocity dispersion of stars  within the bulge. The galactic bulge mass ($`M_{\mathrm{GB}}`$) is derived from its stellar luminosity using the mass–luminosity formula $`\left(M/M_{}\right)/\left(L/10^{10}L_{}\right)^{0.15}`$ . For the S0 galaxy NGC 2300, we assume a typical S0 bulge–to–disk ratio of $`m_V(\mathrm{bulge})m_V(\mathrm{total})=0.60`$ . We note that the bulge mass loss rate in giant elliptical galaxies is $`[M_{\mathrm{GB}}/(10^{12}M_{})]M_{}/\mathrm{year}`$  and assume that the rate of accretion onto the central black hole is an order of magnitude less, as typically found for the Bondi rate ; based on this, we estimate that the accretion rate is $`\mathrm{d}M/\mathrm{d}t0.1M_{12}(M_{}/\mathrm{year})`$, where $`M_{12}[(M_{\mathrm{GB}})/(10^{12}M_{})]`$. In Eddington units, this rate is then

$$\dot{m}c^2(\mathrm{d}M/\mathrm{d}t)/L_{\mathrm{Edd}}0.45(M_{12}/M_8),$$
(2)
where $`M_8[(M_{\mathrm{BH}})/(10^8M_{})]`$, and $`L_{\mathrm{Edd}}=1.3\times 10^{46}M_8`$ ergs/s is the Eddington luminosity limit.

From equations 2—4 in the last paper of Ref. , we obtain that, for losses dominated by curvature radiation, the maximum proton energy expected via the dynamo action considered is

$$E_{\mathrm{max}}=77(\mathrm{d}M/\mathrm{d}t)^{1/8}\left(M_8\right)^{1/4}\mathrm{EeV}\mathrm{\hspace{0.17em}58}\left(M_{12}\right)^{1/8}\left(M_8\right)^{1/4}\mathrm{EeV}.$$
(3)
We note that the values of $`E_{\mathrm{max}}`$ for all 12 galaxies in Table 4 lie above 40 EeV, the lower limit characterizing the AGASA sample considered.

A lower limit to the radiation length ($`\mathrm{\Lambda }_{\mathrm{min}}`$) for proton energy loss associated with photo-pion production is estimated by considering the population of target photons within the source region $`\left[R(\mathrm{source}\mathrm{radius})2GM/c^2\right]`$ at radio frequencies $`\nu 360(\gamma /10^{11})^1`$ GHz , where $`\gamma `$ is the Lorentz factor given by $`\gamma =(E_{\mathrm{max}})/(938\mathrm{M}\mathrm{e}\mathrm{V})`$. The radio estimates required are extrapolated from data at lower frequencies . Those QRs with $`\mathrm{\Lambda }_{\mathrm{min}}>R`$ are expected to successfully accelerate protons up to $`E_{\mathrm{max}}`$.

Apart from their apparent correlation with UHECRs, there is no obvious systematic difference between the first four QR candidates listed at the top of Table 4 and the eight remaining ones at the bottom. However, it’s important to note that we do not as yet know the spin states of the supermassive black holes associated with the nuclei of these galaxies; this key parameter might differ substantially among them. In each instance, the present state of spin depends on the specific prior history involved (e.g., accretion evolution, merger interactions, and earlier activity). Since the galactic nuclei considered here are X-ray dark, as are most (see first Ref. in ), they are not viable candidates for black hole spin determination by means of the accretion disk iron K-line x-ray florescence that appears to be so promising a spectroscopic tool for active galactic nuclei, particularly Seyferts . We emphasize that a black hole state of near maximal spin is a necessary condition for the realization of a compact black hole dynamo of the sort envisaged for accelerating UHECRs . In this sense, if confirmed, the correlation of UHECRs and QRs might well signal the introduction of a new means for identifying those nearby isolated non-active galactic nuclei that harbor highly spun-up black holes .

## 4 Concluding remarks

As recently suggested, quasar remnants are plausible sites of black-hole based compact dynamos that could be capable of accelerating protons up to ultra-relativistic energies. We have found that nearby quasar remnant candidates present an above-random positional correlation with the sample of UHECRs. The correlation appears on closer angular scales than those expected when taking into account the deflection caused by typically assumed intergalactic or Galactic magnetic fields. Possible scenarios producing this effect were discussed; if real, the plausible fine structure of the Galactic field may ultimately provide the basis for the most natural explanation. In order to substantiate and further investigate the apparent correlation reported here between QR candidates and CR arrival directions, we need a large, reliable sample of the most energetic UHECRs, those of the very highest magnetic rigidity. It is hard to claim a definitive correlation with few objects and CRs constituting the samples. Future experiments, such as the Pierre Auger observatory , EUSO  and the NASA space-borne OWL (Orbiting Wide-angle Light-collectors) mission , should vastly increase the availability of UHECRs $`100`$ EeV. As we have already noted, however, CR sources in rich galactic clusters (those with pervasive micro-Gauss fields) are not well suited for our present kind of correlative investigation, even though the nearest QR candidates of interest reside in such systems. We suggest that QR candidates located within the relatively nearby Virgo and Fornax clusters  might be best studied by means of the TeV curvature radiation expected from the putative compact black-hole dynamos associated with these objects (see second paper in ). This work suggests that QRs should also be taken into account when analyzing coincidences in the forthcoming Auger Observatory.

### Acknowledgments

DFT acknowledges L. Anchordoqui and G. Sigl for discussions and criticism, and support from Fundación Antorchas, CONICET, and Princeton University. Much appreciated questions posed by the referee have led to substantial improvement of our paper.


text:
# Mode-coupling theory for structural and conformational dynamics of polymer melts

## Abstract

A mode-coupling theory for dense polymeric systems is developed which unifyingly incorporates the segmental cage effect relevant for structural slowing down and polymer chain conformational degrees of freedom. An ideal glass transition of polymer melts is predicted which becomes molecular-weight independent for large molecules. The theory provides a microscopic justification for the use of the Rouse theory in polymer melts, and the results for Rouse-mode correlators and mean-squared displacements are in good agreement with computer simulation results.

Polymeric materials often find applications as amorphous solids or glasses. During their processing, the structural relaxation in polymeric melts plays an important role because it affects transport and viscoelastic properties. An understanding of conformational dynamics is also required for single-chain diffusional processes. Whereas a number of phenomenological concepts have been invented, no microscopic theory exists which successfully explains the slow structural and conformational relaxations of polymeric melts close to the glass transition, as for example studied in detailed computer simulations of model systems Binder et al. .

We present a first-principles theory for structural and conformational dynamics of unentangled polymer chains using the site formalism which has already provided a good description of their static packing Schweizer and Curro (1997). Building upon the mode-coupling theory (MCT) for simple and low-molecular systems Götze (1991); Götze and Sjögren (1992); Chong and Götze , our macromolecular extension captures the segmental “cage effect” which causes a (asymptotically) molecular-weight-independent ideal glass transition driven by steric hindrance on microscopic length scales. A decoupling of local collective structural relaxation from single-chain conformational fluctuations is obtained, leading to a first-principles derivation of the Rouse model Doi and Edwards (1986). This widely used model assumes a chain to be in a Markovian heat bath, and therefore its applicability requires a time-scale decoupling as our theory achieves. Entanglement effects Doi and Edwards (1986) are neglected in our work because we do not ensure uncrossability of chains and consider only local isotropic forces Schweizer (1989).

A system of $`n`$ chains, consisting of $`N`$ identical monomers or segments, distributed with density $`\rho `$ is considered. In the site formalism Chandler and Andersen (1972), the structural variables are the monomer-density fluctuations for wave vector $`\stackrel{}{q}`$, $`\rho _\stackrel{}{q}^a=_{i=1}^n\mathrm{exp}(i\stackrel{}{q}\stackrel{}{r}_i^a)`$, where $`\stackrel{}{r}_i^a`$ denotes the position of $`a`$th monomer in $`i`$th chain. Structural dynamics shall be described by coherent density correlators, $`F_q^{ab}(t)=\rho _\stackrel{}{q}^a(t)^{}\rho _\stackrel{}{q}^b/n`$. Here $`\mathrm{}`$ denotes canonical averaging for temperature $`T`$. A single (or tagged) polymer (labeled $`s`$) exhibits density fluctuations $`\rho _{\stackrel{}{q},s}^a=\mathrm{exp}(i\stackrel{}{q}\stackrel{}{r}_s^a)`$, and correlators, $`F_{q,s}^{ab}(t)=\rho _{\stackrel{}{q},s}^a(t)^{}\rho _{\stackrel{}{q},s}^b`$, characterize the single-chain dynamics.

We apply MCT equations in the site representation Chong and Götze  to flexible macromolecules expecting them to capture intermolecular caging in polymeric melts. The required inputs are the static structure factors $`S_q^{ab}=F_q^{ab}(0)`$, $`w_q^{ab}=F_{q,s}^{ab}(0)`$, and the direct correlation functions $`c_q^{ab}`$ Chandler and Andersen (1972). There are severe difficulties to solve the $`N\times N`$-matrix MCT equations for polymers because the degrees of polymerization $`N`$ of interest are large. The simplifications adopted here are to neglect chain-end effects for $`c_q^{ab}`$ Schweizer and Curro (1997), $`c_q^{ab}=c_q`$, and to consider the site-averaged correlator, $`F_q(t)=(1/N)_{a,b=1}^NF_q^{ab}(t)`$, which deals with the total monomer-density fluctuations. This mean-field like approximation replaces the site-specific surroundings of a segment by an averaged one. It is supported by the observation that $`S_q=F_q(0)`$ captures the static correlations on the segmental length scale Schweizer and Curro (1997).

We find a set of scalar equations for the normalized coherent correlator $`\varphi _q(t)=F_q(t)/S_q`$:

$`_t^2\varphi _q(t)+\mathrm{\Omega }_q^2\varphi _q(t)`$
$`+\mathrm{\Omega }_q^2{\displaystyle _0^t}𝑑t^{}m_q(tt^{})_t^{}\varphi _q(t^{})=0,`$ (1)
$`m_q(t)={\displaystyle \frac{1}{2}}{\displaystyle 𝑑\stackrel{}{k}V(\stackrel{}{q};\stackrel{}{k},\stackrel{}{p})\varphi _k(t)\varphi _p(t)}.`$ (2)
Here $`\mathrm{\Omega }_q^2=q^2v^2/S_q`$ and $`V=\rho _mS_qS_kS_p\{\stackrel{}{q}[\stackrel{}{k}c_k+\stackrel{}{p}c_p]\}^2/(2\pi )^3q^4`$ with $`v`$ denoting the monomer thermal velocity, $`\rho _m=N\rho `$, and $`\stackrel{}{p}=\stackrel{}{q}\stackrel{}{k}`$. These equations are formally identical to the ones for simple systems Götze (1991). On the other hand, one finds for the single-chain dynamics:

$`_t^2𝐅_{q,s}(t)+𝛀_{q,s}^2𝐅_{q,s}(t)`$
$`+𝛀_{q,s}^2{\displaystyle _0^t}𝑑t^{}𝐦_{q,s}(tt^{})_t^{}𝐅_{q,s}(t^{})=\mathrm{𝟎},`$ (3)
$`m_{q,s}^{ab}(t)={\displaystyle \underset{c}{}}{\displaystyle \frac{w_q^{ac}}{q^2}}{\displaystyle 𝑑\stackrel{}{k}V_s(\stackrel{}{q};\stackrel{}{k},\stackrel{}{p})F_{k,s}^{cb}(t)\varphi _p(t)}.`$ (4)
Here $`𝛀_{q,s}^2=q^2v^2𝐰_q^1`$ and $`V_s=\rho _m(\stackrel{}{q}\stackrel{}{p}/q)^2S_pc_p^2/(2\pi )^3`$.

A traditional description of the single-chain dynamics is in terms of Rouse modes Doi and Edwards (1986). For discrete chains, Rouse-mode correlators are defined by $`C_{pp^{}}(t)=\stackrel{}{X}_p(t)\stackrel{}{X}_p^{}/3N`$ with $`\stackrel{}{X}_p=\sqrt{2/N}_{a=1}^N\stackrel{}{r}_s^a\mathrm{cos}[(a1/2)p\pi /N]`$. Since $`\rho _{\stackrel{}{q},s}^a1+i\stackrel{}{q}\stackrel{}{r}_s^a`$ for small $`\stackrel{}{q}`$, $`C_{pp^{}}(t)`$ can be expressed as a linear combination of $`F_{q,s}^{ab}(t)`$ for $`q0`$. Here, it is neither assumed that $`C_{pp^{}}(t)`$ are diagonal, nor that the decay is exponential. Therefore, our microscopic theory can test the validity of the Rouse theory.

The polymers considered here shall be modeled as follows. First, chains are assumed to be Gaussian, for which $`w_q^{ac}=\mathrm{exp}[q^2|ac|b^2/6]`$ with the statistical segment length $`b`$. Second, each monomer is modeled as a hard sphere of diameter $`d`$ and mass $`m`$, and we set $`b=d`$. All equilibrium properties are then specified by the packing fraction $`\phi =\pi \rho _md^3/6`$ and the degree of polymerization $`N`$. Third, $`S_q`$ and $`c_q`$ are evaluated from the polymer RISM theory Schweizer and Curro (1997). Let us add that the value of $`\phi `$ can become larger than 1 for very dense systems since corrections due to unphysical intrapolymer monomer overlap Schweizer and Curro (1988) are not taken into account in the present work. From here on, the units will be chosen so that $`d=v=m=1`$.

The derived MCT equations (1) and (2) exhibit an ideal liquid-glass transition upon increasing $`\phi `$ at a critical value $`\phi _c`$ as discussed for simple systems Götze and Sjögren (1992); Franosch et al. (1997). The solid line in Fig. 1 exhibits the resulting glass-form factors $`f_q^c=\varphi _q(t\mathrm{})`$ at $`\phi _c`$ for $`N=10`$, and the dashed line $`S_q`$ multiplied by 0.15. The $`f_q^c`$ also measure the amplitude of the slow structural relaxation in the liquid state. Inspection of Eq. (2) shows that the liquid-glass transition is driven mainly by the changes of $`S_q`$ for $`q_{\mathrm{max}}7`$ where $`f_q^c`$ is maximal, i.e., by local fluctuations connected to the average nearest-neighbor-monomer distance (the cage effect). With varying $`q`$, $`f_q^c`$ oscillates in phase with $`S_q`$ as found in simple systems Götze and Sjögren (1992); Franosch et al. (1997). The circles denote the result of molecular-dynamics simulation performed for a similar system Bennemann et al. (1998); Bennemann et al. (1999a, b); Aichele and Baschnagel (2001). The agreement is semiquantitative especially for $`qq_{\mathrm{max}}`$, which is the relevant $`q`$ range for the ideal glass transition dis . The inset demonstrates that $`\phi _c(N)`$ becomes independent of $`N`$ for large $`N`$. This is because the glassy arrest is driven by the local cage effect, and the global chain size plays a gradually smaller role for larger $`N`$ McKenna (1989).

For packing fractions $`\phi `$ close to but below $`\phi _c`$, correlator $`\varphi _A(t)`$ of any variable $`A`$ coupling to density fluctuations exhibits slow structural relaxation and decays in two steps: the decay towards the plateau $`f_A^c`$, followed by the decay to zero (so-called $`\alpha `$-decay) Götze (1991); Götze and Sjögren (1992). Detailed analysis of $`\varphi _q(t)`$ and $`F_{q,s}^{ab}(t)`$ will be presented in later publications. Here, only the Rouse-mode correlators $`C_{pp^{}}(t)`$ shall be studied. Typically, their off-diagonal elements are found to be much smaller (at most only a few %) compared to the diagonal ones at all times. Therefore, the normalized diagonal elements, $`c_p(t)=C_{pp}(t)/C_{pp}(0)`$, shall be considered. Figure 2 exhibits representative results for $`N=10`$ for a reduced packing fraction $`(\phi \phi _c)/\phi _c=10^2`$.

The curves shown in Fig. 2 do not clearly exhibit the two-step-relaxation scenario. This is because the plateaus of $`c_p(t)`$ are so large, $`f_p^c>0.9`$, that only less than 10% of the decay are left for the relaxation towards the plateau. Thus, most of the relaxation of $`c_p(t)`$ occurs in the $`\alpha `$ regime. Let us characterize the $`\alpha `$-time scale by $`c_p(\tau _p)=f_p^c/20`$. The corresponding $`\alpha `$-time scale $`\tau _{\mathrm{qmax}}`$ shall be introduced for the coherent correlator $`\varphi _q(t)`$ for $`q=q_{\mathrm{max}}`$, characterizing the local dynamics of the surrounding medium. The ratio $`\tau _p/\tau _{\mathrm{qmax}}`$ is shown in the inset as a function of $`p`$. For small $`p`$, the scales $`\tau _p`$ are much larger than the local scale $`\tau _{\mathrm{qmax}}`$.

The Rouse theory assumes that all dynamical correlations in the surroundings are much faster than the single-chain dynamics Doi and Edwards (1986). Since a polymer is surrounded by identical polymers, the assumption of the time-scale separation cannot be justified a priori. Our microscopic theory verifies this central assumption. The Rouse theory predicts, within our units, $`\tau _p=\zeta /12\mathrm{sin}^2(p\pi /2N)`$ where $`\zeta `$ denotes the monomer friction. This formula, using $`\zeta `$ as a fit parameter, is shown as the dotted line, and we found $`\zeta /\tau _{\mathrm{qmax}}=16.5`$.

The shape of $`c_p(t)`$ in the $`\alpha `$ regime is often characterized by the stretching exponent $`\beta _p`$ of the Kohlrausch-law fit: $`c_p(t)\mathrm{exp}[(t/\tau _p^{})^{\beta _p}]`$. We found $`\beta _p>0.9`$ for all $`p`$ as shown in the inset. For small $`p`$, $`\beta _p`$ is close to 1 due to the large separation of the scales $`\tau _p`$ and $`\tau _{\mathrm{qmax}}`$ as discussed above. $`\beta _p`$ decreases as $`p`$ increases since for large $`p`$ the scales $`\tau _p`$ become comparable to $`\tau _{\mathrm{qmax}}`$. However, $`\beta _p`$ remains close to 1 even for large $`p`$. This is because their plateaus $`f_p^c`$ are high as explained in Ref. Götze et al., 2000. Thus, all our $`c_p(t)`$ exhibit nearly Debye-relaxation as assumed by the Rouse theory. Let us note that the found features for $`c_p(t)`$ hold also for $`N=100`$, and are in agreement with simulation results Bennemann et al. (1999a).

Of interest are the monomer mean-squared displacement (MSD) $`g_a(t)=[\stackrel{}{r}_s^a(t)\stackrel{}{r}_s^a(0)]^2`$ and the center-of-mass MSD, $`g_C(t)`$. Let us also introduce the monomer-averaged MSD, $`g_M(t)=(1/N)_{a=1}^Ng_a(t)`$. Figure 3 exhibits $`g_M(t)`$ and $`g_C(t)`$ for $`N=10`$ and 100 for $`(\phi \phi _c)/\phi _c=10^2`$. For short times, the MSD exhibits ballistic motion: $`g_X(t)t^2`$ ($`X=M`$ or $`C`$). As the time increases, the MSD begins to be suppressed due to the cage effect, and there appears the so-called $`\beta `$ regime where $`g_X(t)`$ is close to the square of the critical localization length, $`g_X(t)6(r_X^c)^2`$. For this regime, there holds Götze (1991); Götze and Sjögren (1992); Fuchs et al. (1998)

$$g_X(t)=6(r_X^c)^26h_XG(t),X=M,C,$$
(5)
where $`G(t)`$ denotes the $`\beta `$ correlator and $`h_X`$ the critical amplitude. Thus, $`g_M(t)`$ and $`g_C(t)`$ cross their plateaus $`6(r_X^c)^2`$ at the same time. The value $`r_M^c0.077`$, quantifying the monomer localization, for both $`N=10`$ and 100 is consistent with Lindemann’s melting criterion. $`r_C^c`$ is reduced by about $`1/\sqrt{N}`$ compared to $`r_M^c`$, as expected for independent motions of constituent monomers. The asymptotic law (5) for each MSD is drawn as dashed line, and its range of validity is indicated by diamonds. The trends in $`g_X(t)`$ up to the $`\beta `$ regime are qualitatively similar to the ones for a sphere in a simple system Fuchs et al. (1998).

The increase of $`g_X(t)`$ above the plateau $`6(r_X^c)^2`$ towards the diffusion asymptote, $`g_X(t)6Dt`$ with the diffusivity $`D`$, is the $`\alpha `$ process of the MSD. In contrast to $`g_C(t)`$, $`g_M(t)`$ in this regime is significantly affected by chain connectivity since the segments participate in the conformational motion and most of the relaxation of $`c_p(t)`$ occurs here as explained above. As a result, there appears a subdiffusive ($`t^x`$) regime in $`g_M(t)`$. For $`N=10`$, the exponent $`x`$ is 0.60, which is close to the value (0.63) found in the mentioned simulations Bennemann et al. (1999a). For $`N=100`$, we find $`x=0.5`$ as predicted by the Rouse model and by the asymptotic evaluation of our theory, and the different value for $`N=10`$ can be attributed to finite $`N`$ effects. Thus, we find a strong polymer-specific effect for the beginning of the $`\alpha `$ process, while no such effect is reflected in $`g_M(t)`$ up to the $`\beta `$ process.

The insets of Fig. 3 show the ratio of the MSD for end and central monomers, $`g_1(t)/g_{N/2}(t)`$. The ratio starts from 1 in the ballistic regime, exhibits a maximum for intermediate times, and tends to 1 in the diffusion regime. It is seen that the ratio remains close to 1 also within the $`\beta `$ regime. This is because the dynamics here does not reflect the chain connectivity as explained above. The Rouse theory predicts the maximum to be 2 within the time regime where the monomer MSD exhibits the $`t^{1/2}`$ law, i.e. in the $`\alpha `$ regime. The result for $`N=100`$ clearly indicates this behavior, while the maximum for $`N=10`$ is slightly smaller (1.8) due to finite $`N`$ effects. The shape of the ratio for $`N=10`$ in the $`\beta `$ and $`\alpha `$ regimes is in semiquantitative agreement with the mentioned simulation results Baschnagel . Thus, the chain-end effect is properly taken into account in our theory even with neglecting that effect for $`c_q^{ab}`$. This is because the matrix structure of Eqs. (3) and (4) are preserved for describing the single-polymer dynamics.

In summary, a first-principles theory for structural slowing down of dense polymer systems has been presented which also provides a microscopic derivation of the Rouse model for unentangled chain melts. Chain connectivity is seen to cause polymer specific long-time anomalies of the $`\alpha `$-process. We use the concept of an ideal MCT glass transition, familiar for simple liquids and colloidal suspensions Götze (1999), and our results agree well with simulation studies of simple (coarse grained) polymer models. Our theory thus provides insights into dynamical aspects typical for polymer melts which are also observed in models van Zon and de Leeuw (1998, 1999) with a more realistic local chemistry.

###### Acknowledgements.

We thank J. Baschnagel, M. Aichele, and W. Götze for discussions. M.F. was supported by the Deutsche Forschungsgemeinschaft, grant Fu 309/3.


text:
# Incommensuration Effects and Dynamics in Vortex Chains

\[

## Abstract

We examine the motion of one-dimensional (1D) vortex matter embedded in a 2D vortex system with weak pinning using numerical simulations. We confirm the conjecture of Matsuda et al. (Science 294, 2136 (2001)) that the onset of the temperature induced motion of the chain is due to an incommensuration effect of the chain with the periodic potential created by the bulk vortices. In addition, under an applied driving force we find a two stage depinning transition, where the initial depinning of the vortex chain occurs through soliton like pulses. When an ac drive is added to the dc drive, we observe phase locking of the moving vortex chain.

\]

The static and dynamical behavior of elastic media interacting with periodic or random substrates is a topic of intense study as it relates to a wide variety of physical systems, such as friction where atoms move over a periodic potential created by the underlying stationary atoms, sliding charge-density waves , and vortex motion in superconductors with random quenched disorder  or periodic hole or dot arrays . A simplified case occurs when the motion is confined in 1D or quasi 1D systems. An example of this is vortex motion in superconductors with nano fabricated arrays of channels . There are also naturally occurring defects that can act as quasi-1D channels for vortex motion, such as grain  and twin boundaries . The vortices in fabricated channels experience two types of pinning: the first from the intrinsic random disorder in the sample, and the second from the periodic potential created by the immobile vortex lattice surrounding the channels. Experiments and simulations in these systems have observed commensuration effects and interesting dynamical states .

Another intriguing example of quasi 1D vortex matter is the vortex chain state which occurs in layered superconductors with tilted magnetic fields as first imaged by Bolle et al. . Here a portion of the vortices align in a string with a smaller lattice spacing than that of the bulk vortices. Recently scanning Hall probe  and Lorentz microscopy  images have shown that the vortex chain states exhibit a remarkably rich variety of behaviors and structures for various regimes. The recent Lorentz microscopy experiments found that for low temperatures, both the vortices along the chain state and the bulk vortices are stationary. For increasing temperatures the vortices along the chain disappear, which is speculated to be due to the onset of motion or oscillation of the vortex chains, so that the chain vortex images are blurred while the bulk vortices are still stationary . It was conjectured in  that this effect is due to an incommensuration effect where vortices in the chain move in a periodic potential created by the bulk vortices. Since the vortices in the chain have a different spacing than the periodic potential, the chain is incommensurate and hence more weakly pinned than if the chain were commensurate. These images also showed that the blurring first occurs where there a mismatch accompanied by a topological defect in the chain. Additional evidence for this scenario is that as the field is increased, the periodicity of the potential decreases and the incommensuration effects are enhanced, causing the disappearance of the chains to shift to lower temperatures, as observed. Dodgson  has proposed an alternative explanation in which the image disappears due the sudden formation of tilted chains. In addition to the intriguing physics of the vortex matter in confined geometries it is also of paramount importance to understand vortices in the chain state or along grain and twin boundaries, as these can act as weak lines where flux motion occurs, leading to the breakdown of superconductivity in the chains prior to the bulk.

In order to examine the conjecture that the disappearance of the vortex chain is due to an incommensuration effect, as well as to explore dynamical consequences of these incommensuration effects, we have performed numerical simulations of vortex chains coexisting with bulk vortices in the presence of weak quenched disorder and thermal noise. We also investigate the effects of applied dc and combined dc and ac external drives on the chain motion. Our results should be relevant for the vortex chain state where the vortices remain aligned in the $`z`$-direction so that much of the physics can be considered 2D. In addition our results are relevant for vortex motion in artificial channels and grain or twin boundaries where there are additional vortices outside the channels.

We consider a 2D sample with periodic boundary conditions in the $`x`$ and $`y`$ directions. We model the vortices as repulsive point particles which interact with weak random quenched disorder, applied driving forces, and thermal noise. A portion of the vortices are confined to move along a quasi-1D channel. The overdamped equation of motion for a vortex $`i`$ is

$$\eta \frac{\mathrm{𝐝𝐫}}{dt}=𝐟_i=𝐟_i^{vv}+𝐟_i^{rp}+𝐟_i^d+𝐟_i^T$$
(1)
Here $`\eta `$ is the damping constant that is set equal to unity. The vortex-vortex interaction force is

$`𝐟_i^{vv}=_{ji}^{N_v}_iU_v(r_{ij})`$, where the vortex interaction potential is $`U_v=A_v\mathrm{ln}(r)`$. The force from random quenched disorder $`𝐟_i^{rp}`$ comes from randomly placed attractive parabolic traps of range $`r_p`$ with $`r_p/a=0.1`$, where $`a`$ is the vortex lattice constant and the traps have a maximum force of $`f_p=1.0`$. The driving forces $`𝐟_i^d`$ include a dc force $`𝐟_i^{DC}`$ which is modeled as a uniform drive on the vortices in the $`x`$-direction. We start from zero driving and with small increments increase the dc driving force. At each increment we average the vortex velocity $`V_x`$ over a fixed time interval. An ac drive can also be applied in the form $`𝐟_i^{AC}=A\mathrm{sin}(\omega t)\widehat{𝐱}`$, where $`A`$ is the ac amplitude and $`\omega `$ is the frequency. The forces from the thermal noise $`𝐟_i^T`$ come from random Langevin kicks with $`<f^T(t)>=0`$ and $`<f^T(t)f^T(t^{^{}})>=2\eta k_BT\delta (tt^{})`$. We measure temperature in units of $`T_m`$, the temperature at which the bulk vortices melt. For high-$`T_c`$ materials at fields where the chain state is observed, $`T_m`$ is around $`7080K`$ and for the low-$`T_c`$ materials with nano-channels, $`T_m`$ would correspond to $`T_c`$. Lengths are measured in units of $`a`$, and forces in $`A_v`$, the vortex-vortex interaction prefactor.

We first consider the effects of temperature in systems with no external driving force to test the conjecture in . We start the system at $`T/T_m=0.0`$ where the bulk vortices form a triangular lattice with lattice constant $`a`$. Additional vortices are added along the channel, giving a vortex spacing in the channel of $`a^{}<a`$. We monitor the vortex displacements for fixed time intervals of $`10^4`$ MD steps, and measure the vortex positions and trajectories.

For low $`T`$ there is little motion of the vortices in the chain or the bulk as seen in Figure 1(a). In addition the incommensuration in the chain is stationary or pinned in a well defined location. For increasing temperatures there is a transition to a state in which the vortex chain exhibits a higher mobility than the bulk vortices, as seen in the vortex trajectory traces of Fig. 1(b) at $`T/T_m=0.25`$.

To elucidate the exact nature of the motion along the chain, in Fig. 2(a) we plot the position of a single vortex in the chain for the states in Fig. 1(a) and Fig. 1(b). For the stationary chain the vortex remains immobile with only small thermal fluctuations in position of magnitude less than $`a`$. For the mobile chain state the vortex does not move continuously, but instead moves in discrete jumps of magnitude nearly $`a`$. These jumps do not occur throughout the chain simultaneously, but run through the chain as a soliton like pulse with the vortices jumping sequentially. The discrete jumps occur due to the depinning of the incommensurations, which move randomly through the sample, causing each vortex to jump one lattice constant as the incommensuration passes over it. The discreteness of the jumps arises from the periodic potential in which the chain vortices sit, created by the stationary ordered bulk vortices. Since the incommensurations are much more weakly pinned than the bulk vortices, they depin at a much lower temperature.

Figure 2(a) also shows that although there is motion of vortices along the chain, a single vortex does not diffuse but only moves back and forth by a single lattice constant. This can be understood by considering a single randomly moving incommensuration. When the incommensuration traveling in the positive $`x`$-direction passes over a vortex, the vortex moves by $`a`$ in the positive $`x`$-direction. The same incommensuration cannot, however, later cause the same vortex to move by an additional distance in the positive $`x`$-direction, but only in the negative $`x`$-direction when the incommensuration returns. Thus a single randomly moving incommensuration produces no vortex diffusion along the chain but rather an oscillation. If a series of incommensurations all moving in the same direction pass over a vortex, the vortex could move several lattice constants in one direction. This is prevented because there is no net shear in the system. The moving incommensuration sets up a strain field that prevents the incommensuration from traveling large distances in the same direction. As a result, the incommensuration does not diffuse freely through the chain but tends to drift back and forth about its original starting position.

Using a Voronoi construction, we find that the incommensuration in the vortex chain coincides with a local defect composed of one 5-fold and two 7-fold coordinated vortices. It is the 1D motion of this defect that occurs in Fig. 1(b). The images in Ref.  also indicate that the initial blurring in the chain occurs where there are localized topological defects in the vortex lattice. This agrees well with our finding that it is the motion of defects that is occurring along the chains.

In Fig. 2(b) we show the power spectra for the mobile and immobile chain. For the immobile chain, the noise spectrum is white, while for the mobile chain, the low frequency noise power is several orders of magnitude higher, due to the correlated low frequency motions of the incommensuration, and there is now a $`1/f^\alpha `$ signature, with $`\alpha =1.6`$. For $`T>T_m`$, when all the vortices diffuse randomly, the noise spectra again becomes white.

In simulations where the chain is commensurate with the bulk vortex lattice, such that $`a=a^{}`$, the chain and bulk vortex mobility are the same. Thus, the increased mobility of the vortices along the chain, as shown in Fig. 1(b), is due strictly to an incommensuration effect. When the vortex density along the chain is increased, lowering $`a^{}`$ and increasing the number of incommensurations, the thermal depinning temperature of the chain decreases. The depinning and phase locking described below persist as the number of incommensurations increases.

A further probe of the relative pinning strength of the incommensurations along the chain, compared to the

pinning of the bulk vortices, is the application of a dc drive. We consider the case with weak random pinning where, in the absence of a chain, there is a single well defined elastic depinning transition. In Fig. 3(a) the velocity-force curve for the incommensurate case shows two depinning transitions; the first is the depinning along the chain, and the second at $`f_d=1.0`$ is the bulk vortex depinning. A similar two stage depinning has been observed in artificial channels and periodic pinning arrays . Imaging experiments in the vortex chain state also find that the chain depins first under an external driving force . In Fig. 3(b) we plot the positions of the vortices in the chain as a function of time just above the depinning transition, showing that the motion is not continuous, but occurs in a pulse or stick-slip motion. Individual vortices remain stationary most of the time, and then jump one lattice constant in sequence. The disturbance or incommensuration moves much faster than the individual vortices. This is the same type of motion that occurs during the thermal depinning; however, in the driven case the incommensuration moves in only one direction. A similar type of 1D soliton motion has been

seen for vortex motion in 2D periodic pinning arrays in simulations  and experiments  where the vortices effectively channel along the pinning rows.

We have also investigated the effects of a combined dc and ac drive on the moving chain state. Since the motion along the chain occurs over an effective periodic substrate, interference effects should occur between the ac drive and the frequency induced by the dc drive. In this case, it is the incommensurations, rather than the individual vortices in the chain, that undergo phase locking. In Fig. 4(a) we show a series of velocity-force curves for the moving chain for different ac amplitudes and fixed ac frequency. Here clear phase locking effects can be seen as steps in the velocity curves. The widths of the critical depinning force and the higher order steps vary nonmonotonically with ac amplitude. If the phase locking steps are of the Shapiro step type, the $`n`$th step should vary as the Bessel function $`J_n(A)`$ . In Fig. 4(b) and Fig. 4(c), the widths of the depinning force $`\eta _0`$ and the first step $`\eta _1`$ fit to $`J_0`$ and $`J_1`$, respectively, showing that the phase locking is best described as a Shapiro step similar to those seen in sliding CDW’s  as well as 2D vortex systems with periodic pinning . The steps shown here occur below the bulk depinning transition, for $`f_d/f_p<1.0`$, where only the chain vortices are depinned.

In conclusion, we have investigated vortices confined in 1D embedded in an ordered 2D vortex lattice with weak pinning. The vortices in the chain are incommensurate with the periodic potential created by the bulk vortices. We confirm the conjecture in Ref.  that the smearing of the vortex chain state is due to an incommensuration effect, where for increasing temperature the incommensurations thermally depin and become mobile before the bulk vortices do. The thermal motion along the chain occurs through a soliton like pulse which moves randomly through the sample. The vortices hop by one lattice constant in the positive or negative direction as the pulse moves through. The onset of motion along the chain can also be observed as an increase in the voltage noise with a $`1/f^\alpha `$ characteristic spectra. Under application of an external dc drive, we observe a two stage depinning where the vortices in the chain depin first and form running solitons. When an ac drive is added to the dc drive, we observe phase locking of the moving incommensurations in the form of Shapiro steps in the velocity-force curves. Our predictions can be tested through noise measurements in the vortex chain state under application of external drives or temperature, and also apply to vortex motion in artificial channels.

Note: Just before submitting this work we became aware of recent experiments that find phase locking of dc and ac driven vortices in artificial channels .

This work was supported by the US Department of Energy under Contract No. W-7405-ENG-36. We thank P. Kes and A. Tonomura for useful discussions.