text:
# A Point Source Excess in Abell 1185: Intergalactic Globular Clusters?Based on observations with the NASA/ESA Hubble Space Telescope obtained 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 Based in part 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 NASA, and was made possible by the generous financial support of the W. M. Keck Foundation.

## 1 Introduction

Although most of the stars residing in galaxy clusters are contained within distinct galaxies, it has been known since the study of Zwicky (1951) that clusters often have diffuse stellar components that are intergalactic in nature. The existence of this low surface brightness starlight has been firmly established by numerous subsequent studies ($`e.g.,`$ Thuan & Kormendy 1977; Melnick, White & Hoessel 1977; Vílchez-Gómez, Pelló & Sanahuja 1994; Bernstein $`etal.`$1995) but, due to the difficulties inherent in the analysis of extended, low surface brightness features, there remains a surprising level of disagreement on its basic properties including spatial distribution, luminosity and color (see the excellent review in Durrell $`etal.`$2002).

Our ability to characterize the properties of this intergalactic stellar component has improved dramatically in recent years, thanks in large part to advances in instrumentation and to the targetted analysis of tracer populations such as individual red giant branch stars (Ferguson, Tanvir & von Hippel 1998; Durrell $`etal.`$2002), planetary nebulae (Theuns & Warren 1997; Méndez $`etal.`$1997 Feldmeier, Ciardullo & Jacoby 1998; Feldmeier $`etal.`$2002; Okamura $`etal.`$2002), type Ia supernovae (Gal-Yam $`etal.`$2002) and even HII regions (Gerhard $`etal.`$2002). These studies typically find the intergalactic component to comprise 10 – 20% of the total cluster luminosity.

Globular clusters (GCs) represent an additional, and potentially powerful, tracer of this intergalactic component. Indeed, the possibility that galaxy clusters may contain a population of intergalactic GCs (IGCs) has been considered many times (van den Bergh 1958; Muzzio, Martínez & Rabolli 1984; White 1987; West $`etal.`$1995). Such an IGC population might either form in situ, as suggested by West (1993) and Cen (2001), or arise through the disruption or tidal stripping of cluster galaxies. Ample observational evidence shows that these latter processes must play a role in the evolution of galaxy clusters ($`e.g.,`$ Weil, Bland-Hawthorn, & Malin 1997; Gregg & West 1998), as predicted by numerical simulations of cluster evolution (see, $`e.g.`$, Dubinski 1999 and references therein).

If, as suggested by the above studies, 10 – 20% of the total stellar luminosity in clusters takes the form of a diffuse intergalactic component, then it is possible to estimate the total number of IGCs expected within a given cluster. Early-type galaxies in rich cluster environments have GC specific frequencies (Harris & van den Bergh 1981) of $`S_N4`$ (Harris 1991). A cluster such as Virgo, with a total luminosity of $`L_V1.5\times 10^{12}`$ $`L_{V,}`$ inside a radius of $`6^{}`$ (Sandage, Binggeli & Tammann 1985)111For an adopted Virgo distance of 16 Mpc (Ferrarese $`etal.`$1996) should then contain $`10^4`$ GCs associated with its intergalactic component.

However, the unambiguous detection of IGCs has proven difficult, with conflicting claims on whether or not such IGCs might already have been detected in nearby clusters ($`e.g.,`$ West $`etal.`$1995; Harris, Harris & McLaughlin 1998; Côté $`etal.`$2001; Hilker 2001; Marín-Franch & Aparicio 2002). Much of the confusion in this debate stems from the fact that, if IGCs follow the same distribution within a given cluster as the other baryonic components ($`i.e.,`$ cluster galaxies and X-ray emitting gas), they should be most apparent at, or near, the dynamical center of the cluster. Since most galaxy clusters have a giant elliptical galaxy in close projection to their center, the direct detection of IGCs has been hampered by the difficulties involved in disentangling the hypothesized population of intergalactic GCs from those intrinsic to the central galaxy.

The rich galaxy cluster Abell 1185 is peculiar in this respect as its brightest cluster galaxy is offset by $``$ 3, or about 120 kpc, from the centroid of its X-ray gas distribution (Mahdavi $`etal.`$1996). Assuming that the X-ray centroid marks the dynamical center of the cluster mass distribution, and is not the result of an ongoing merger, then the detection of IGCs in Abell 1185 promises to be more straightforward than in most clusters.

In this paper, we present deep HST observations of a field centered on the peak of the X-ray distribution of A1185. These observations were designed to test the IGC hypothesis by searching for the expected point source excess compared to a blank field. As shown below, the observations, which were designed to sample the brightest part of the GC luminosity function at the distance of A1185, do indeed contain a point source excess that starts at the expected magnitude for an old GC population.

## 2 Observations and Analysis

We used WFPC2 on HST to obtain F814W ($`I`$) imaging of A1185 as part of program GO-8164. The total exposure time was 13000 s, divided equally among ten, dithered images. The raw data were processed with the standard STScI pipeline using the best available calibration files. As a comparison field, we used a subset of the Hubble Deep Field North (HDF-N) observations (Williams $`etal.`$1996), chosen to have a total exposure time identical to that of the A1185 field. Given the enhanced noise and smaller field of view of the PC chip compared to the WF chips, we limit our analysis to the latter.

Cosmic ray masks and the shifts required to register the images were obtained using the DITHER package in IRAF.222IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. The images were combined by applying integer shifts and averaging to produce the final images. SExtractor (Bertin & Arnouts 1996) was used to detect and classify the sources present on the field, setting the detection threshold to 5 connected pixels having counts $`2\sigma `$ above the local background. The FWHM of point sources, required to properly classify the detections, was chosen on the basis of artificial star tests (see below) to be 0$`\stackrel{}{\mathrm{.}}`$16.

To carry out profile-fitting photometry, and to aid in setting the parameters for SExtractor, we constructed point spread functions (PSFs) for our frames as follows. We used DAOPHOT (Stetson 1987) to create artificial WFPC2 frames containing 121 stars positioned in a 11$`\times `$11 grid on the image. PSFs for the F814W filter were kindly provided by P.B. Stetson. Each artificial frame was then shifted, creating a series of images that had the same dither pattern as the actual images for both A1185 and HDF-N. Each set of images was combined into a final frame using the same procedures that were used for the actual observations. The final images created in this way thus contain point sources constructed from the single-frame PSF but have been subjected to the same DITHER process as the actual data. DAOPHOT was used to determine PSFs which vary quadratically with position in the frames for both the A1185 and HDF-N frames.

SExtractor assigns each detected object a classification parameter, CLASS, which ranges from $``$ 0 for galaxies to $``$ 1 for point sources. CLASS is obtained using a neural network trained on simulated images (for details see Bertin & Arnouts 1996). We consider point source candidates to have CLASS $`>0.8`$. An optimized value for the FWHM parameter was determined by carrying out artificial star tests on the A1185 and HDF-N images, adding 25 stars at a time (so as not to alter the crowding conditions of the frame) until a total of $`5000`$ stars were added to each chip. The FWHM quoted above was found to be the minimum value that ensures that the bulk of the artificial stars brighter than our detection limit are classified correctly as point sources.

Figure 1 shows the distribution of CLASS parameters for artificial stars with $`I_{F814W}26`$ mag recovered from the A1185 and HDF-N fields. The majority (85%) of these objects have CLASS $`0.8`$. The fraction of misclassified artificial stars depends sensitively on the magnitude cutoff, falling to $``$ 8% for $`I_{F814W}25.5`$ and $``$ 5% for $`I_{F814W}24.5`$. Based on these artificial star tests, we find the $`50\%`$ completeness level lies at $`I_{F814W}26`$ mag for A1185 and $`I_{F814W}26.2`$ mag for HDF-N. As Figure 1 makes clear, the distributions for the two fields are similar, as is necessary if we are to compare directly the counts for the two fields.

PSF-fitting photometry was then carried out with DAOPHOT, taking the sky from an annulus with inner and outer radii of 5 and 10 pixels. Instrumental magnitudes were transformed into the VEGAMAG system using zero-points taken from the HST WFPC2 Data Handbook (Baggett $`etal.`$2002). A correction for foreground extinction was performed using the reddening curves of Cardelli, Clayton & Mathis (1989), with the value of $`E(BV)`$ taken from the DIRBE maps of Schlegel, Finkbeiner & Davis (1998). A final correction of $`0.1`$ mag was applied to correct from a PSF radius of $`0\stackrel{}{\mathrm{.}}5`$ to one of infinite aperture (Holtzman $`etal.`$1995).

Those objects having a DAOPHOT $`\chi `$ statistic333$`\chi `$ is a robust estimate of the ratio of the observed to expected scatter about the model profile (Stetson & Harris 1988). greater than $`1.75`$ were discarded and any point sources brighter than $`3\sigma `$ above the expected turnover magnitude of the globular cluster luminosity function at the distance of A1185 were omitted from the analysis. The threshold on $`\chi `$ was determined on the basis of the artificial star tests. We assumed $`M_{I,TO}=8.4`$ mag and $`\sigma =1.4`$ mag for the absolute magnitude of the turnover and the dispersion of the Gaussian representation of the GC luminosity function, respectively (Harris 2002). The distance to A1185 is calculated directly from its redshift, $`z=0.0325`$, assuming $`H_0=72`$ km s-1 Mpc-1 (Freedman $`etal.`$2001). This procedure results in the rejection of all objects brighter than $`I_{F814W}=23.1`$ mag and eliminates $`9`$ objects in the A1185 field and $`10`$ in HDF-N. For comparison, $`67`$ galactic stars with $`I23`$ are predicted by the Bahcall-Soneira Galaxy model (Bahcall & Soneira 1980; Bahcall 1986) in these fields

As a final check on the point-source nature of the remaining detections, we examined the radial profile of the difference between the observed profile of each source with the PSF profile expected at that position in the field. While a trend between the median difference and radius might then be expected if the sources were extended no such trend was detected, supporting the SExtractor classifications.

## 3 The Nature of the Excess

Our final catalog consists of 99 point sources in the A1185 field. By contrast, we find 12 sources in the HDF-N field, which serves as our control. The luminosity functions for both fields are shown in Figure 2. For A1185 we also show the expected luminosity function for GCs at the distance of A1185, multiplied by the completeness function and normalized to the observed number of sources. Our counts for the HDF-N field are in good agreement with those of Elson, Santiago & Gilmore (1996) who find 14 point sources in the same magnitude interval. In what follows, we examine a number of possible explanations for this point source excess.

### 3.1 Galactic Field Stars

We first consider the possibility that the A1185 field contains a larger number Galactic field stars than does the HDF-N field. However, this explanation appears unlikely given the size of the excess ($`i.e.,`$ nearly an order of magnitude) and the roughly similar Galactic latitudes of the two fields: $`b`$ = 67$`\stackrel{}{\mathrm{.}}`$8 for A1185 and $`b`$ = 54$`\stackrel{}{\mathrm{.}}`$8 for the HDF-N field. Indeed, according to the Bahcall-Soneira Galaxy model (Bahcall & Soneira 1980; Bahcall 1986), both fields are predicted to contain only 4-6 stars in the magnitude range $`23I26`$. We conclude that Galactic field stars are unlikely to be the origin of the observed excess.

### 3.2 Background Galaxy Cluster

Could the excess be due to a distant background galaxy cluster which happens to fall in the A1185 field? The low surface density of high-redshift clusters (e.g. Postman $`etal.`$2002) suggests that this explanation is unlikely. More importantly, we expect almost all background giant ellipticals to be resolved on our WF frames: i.e., the angular diameter distance turns around at $`z1.5`$ in such a way that, for $`z1`$, 1′′ corresponds to $``$ 6 kpc. For instance, Lyman break galaxies at $`z3`$, though physically compact, are still clearly resolved in deep WFPC2 images, with half-light radii $`>`$ 0$`\stackrel{}{\mathrm{.}}`$2 (Steidel $`etal.`$1996). A definitive test should be possible with multi-color data (and particularly IR colors) for the point source population in A1185.

### 3.3 Nuclei of dE,N Galaxies

A third possibility is that our point source catalog includes a population of dE,N galaxies. At the distance of A1185, the bright nuclei of such galaxies would be unresolved. As Binggeli & Cameron (1991) have shown, the nuclei of dE,N galaxies mimic the bright end of the GC luminosity function in the Virgo cluster, so it is conceivable that our sample may include some dE,N galaxies, particularly since our WFPC2 field is located near the center of A1185 where the galaxy surface density is expected to be high. There are two separate issues that must be addressed to test the plausibility of this explanation: (1) the probability that dE,N galaxies would be classified incorrectly by SExtractor as point sources; and (2) overall number of dE,N galaxies expected in our A1185 field. This latter point will be examined in detail below.

To address the first issue, we carried out simulations in which artificial dE,N galaxies were added to our frames. These frames were then run through the same reduction procedures as the actual observations and, at the end of each run, the classifications returned by SExtractor were recorded. The dE,N galaxies were simulated using the data presented in Table 1 of Lotz $`etal.`$(2001) which gives total magnitudes, exponential length scales, nuclei $`V`$-band magnitudes and $`VI`$ colors for a sample of 27 dE,N galaxies in the Virgo and Fornax clusters that have been observed with HST. To simulate the dE,N galaxies, we first used DAOPHOT to add point sources with the tabulated $`I`$ magnitudes after shifting to the distance of A1185. The task MKOBJECTS in IRAF was then used to add an exponential disk having the appropriate exponential scale, ellipticity and total magnitude (minus the contribution of the nucleus). The colors of the exponential disks were taken to be $`(VI)=0.650.022M_V`$ and $`(BV)=0.410.022M_V`$; these relations were obtained by combining a \[Fe/H\]–$`M_V`$ relation for dwarf galaxies (Côté $`etal.`$2000) with the color–metallicity relations presented in Barmby $`etal.`$(2000). Each galaxy in the Lotz $`etal.`$(2001) sample was simulated $`50`$ times per WF chip, with the galaxies being simulated in groups of $`9`$ so as not to overcrowd the fields. These simulations revealed that all dE,N galaxies were classified correctly as non-stellar in all trials. We therefore conclude that the nuclei of dE,N galaxies cannot be the origin of the point source excess in A1185.

### 3.4 GCs associated with Dwarf Galaxies

A related possibility is that the point sources detected in the A1185 field represent bona fide GCs, but ones that are associated with dwarf galaxies that may be present in A1185. We examine this possibility by calculating the expected number of dE and dE,N galaxies in the A1185 field. Since the velocity dispersions of Virgo and A1185 are similar (Binggeli, Tammann & Sandage 1987; Mahdavi $`etal.`$1996), their dwarf populations can be compared directly. Based on the Virgo dE+dE,N luminosity function published by Sandage, Binggeli & Tammann (1985), expect $`1000`$ dE and dE,N galaxies in Virgo out to a distance of $`5^{}`$. We have adopted a cutoff at the faint end of the luminosity function of $`M_B=11.5`$ as there are no known GC systems belonging to dwarfs fainter than that. Given that A1185 is $`9`$ times more distant than Virgo, and assuming that the surface density of dwarfs in Virgo has the approximate form $`e^{r/\alpha }`$ with $`\alpha 1\stackrel{}{\mathrm{.}}7`$ (Ferguson & Sandage 1989) 444This value of $`\alpha `$ is appropriate for dE + dE,N fainter than $`M_B13`$. Brighter dwarfs exhibit in general larger $`\alpha `$ and this would drive down the estimated number of dwarfs in our field., the expected dE surface density at the center of A1185 is roughly 5500 deg-2. Thus, we expect our WFPC2 field to contain $`8`$ dE and dE,N galaxies. Given that these galaxies would typically contain only $`5`$ clusters each ($`e.g.,`$ Lotz $`etal.`$2001), and that our photometry samples only the brightest 16% of the GC luminosity function at the distance of A1185, we expect only $``$ 8 point sources to be GCs associated with dwarf galaxies in A1185.

### 3.5 Faint dE galaxies

We have also explored the possibility that the point source catalog in A1185 has been contaminated by faint dE galaxies with small exponential scale-lengths. This possibility has been investigated by adding simulated dE galaxies to the WFPC2 frames, with the galaxies modeled as pure exponential disks having scale-lengths of $`r_e=0.1`$ kpc. This choice of scale-length corresponds to the low-end of the values found for Local Group dwarfs (Mateo 1998). It is these highly concentrated dwarfs that are cause for the greatest concern, as they are the ones most likely to be misclassified at the distance of A1185. We added 150 artificial dE galaxies having absolute magnitudes $`M_I=11.7,11.2,10.7,10.2,9.7`$ and $`9.2`$ (1200 in total) to each WF chip and analyzed the images in the same way as the actual data. Only for $`M_I=10.7`$ and $`10.2`$ were any dwarfs misclassified, and then only in $``$1% of the cases. Combining this result with the fact that we expect $``$ 8 dE galaxies in this magnitude range (determined from an extrapolation of the fitted Schechter function for the dE population in Virgo; Sandage $`etal.`$1985), we find that the number of compact dE galaxies in our point source catalog is completely negligible.

### 3.6 A1185 Globular Clusters

The calculations above suggest that the number of dwarf galaxies and Galactic field stars expected in the A1185 field fall far short of that needed to explain the point source excess. We conclude that the point source excess is probably the result of bona fide GCs in A1185, but is it plausible that the sources are GC associated with nearby ellipiticals rather with the cluster as a whole? As Figure 3 shows, a few of the sources are close in projection to neighboring galaxies, but the bulk of the population is not obviously associated with galaxies. To examine this issue more quantitatively, we consider the GC systems of the six bright galaxies that are nearest to the WFPC2 field: NGC 3550 (BCG), NGC 3552, CGCG 155-081, 2MASSXi J1110398+284224, MCG+05-27-004 NED01 and MCG+05-27-004 NED02. These six objects are marked by crosses in Figure 3, and labelled according to the above order.

To estimate the number of GCs that are associated with these galaxies and fall in the WFPC2 field, it is first necessary to assume a radial density profile for their GC systems. On theoretical grounds, the profiles of non-central galaxies are expected to be tidally truncated by the mean cluster field ($`e.g.,`$ Merritt 1984), but precise values for the truncation radii of cluster galaxies have not been determined observationally, nor has the process been probed adequately in numerical simulations. As a first estimate, we use the theoretical predictions for the tidal radius put forward in Merritt (1984), assuming $`\alpha =\beta =1.5`$ in his equation (6a). The cluster velocity dispersion is taken to be $`\sigma _v=736`$ km s-1 (Mahdavi $`etal.`$1996) and the velocity dispersion for each galaxy is estimated using the Faber-Jackson relation (Faber & Jackson 1976). For the radial density profiles inside the tidal radius, we use power laws of the form $`r^\alpha `$ with $`\alpha `$ given by

$$\alpha =0.29M_V8.0$$
(1)
(Kaisler $`etal.`$1996). Absolute magnitudes, $`M_V`$, were measured for the A1185 galaxies by running SExtractor on a $`30`$s $`V`$-band exposure taken with the Low-Resolution Imaging Spectrometer (Oke $`etal.`$1995) at the W.M Keck Observatory on 19 March, 1999. To estimate the total number of GCs belonging to each galaxy, we assume a constant GC specific frequency of $`S_N=4`$ (Harris 1991). Integrating the various GC density profiles over the WF chips, and taking into account the fact that we observe only the brightest $`16\%`$ of the GC luminosity function, we find a total contribution of $`15`$ GCs from the GC systems of these galaxies.555 The galaxy 2MASSXi J1110398+284224 lies partly on the WF4, so we neglect any GCs falling within 10′′ of its center since these GCs would go undetected in our reduction pipeline. Of the 15 expected GCs, $`13`$ would be on WF4. Although this detector is indeed observed to contain largest number of point sources (39), the factor of six discrepancy between the observed and predicted counts suggest that the point source excess is unlikely to be due to GCs associated with neighboring galaxies, if the galaxies are tidally truncated as described above.

However, this estimate depends rather sensitively on the assumed recipe for tidal limitation. We have therefore calculated the expected contribution when this effect is neglected entirely. Since we are now interested in the spatial distribution of GCs at very large projected radii, the power-law representations of the radial density is inadequate (Rhode & Zepf 2000). We instead assume that all galaxies have GC surface density profiles similar to the $`r^{1/4}`$ profile of M49. Beyond $`r=100`$ kpc, the profiles are set to zero, consistent with the results of Rhode & Zepf (2000). This is a rather conservative approach, as the fainter galaxies will certainly have less extended profiles (see Eq. 1). Following the same procedure as described above, we then find an expected contribution of $`62`$ GCs, still lower than the observed number of point sources, but now within a factor of $``$ 50% of the actual excess. We caution, however, that the previous estimate is probably more realistic, as the cluster galaxies have almost certainly been tidally truncated at some level. Indeed, dynamical modeling of the GC system of M87 — the central elliptical in the Virgo cluster — has shown the velocity dispersion profile of its GC system to rise beyond $`r20`$ kpc, as expected if the GCs are orbiting in the potential well defined by both the galaxy and its parent cluster (Côté $`etal.`$2001). At a projected distance of 100 kpc, the gravitational potential in which the GCs orbit is dominated not by the central galaxy, but by the cluster itself.

Could a population of IGCs be responsible for the observed point source excess? West $`etal.`$(1995) proposed a phenomenological model in which the high specific frequencies of some centrally dominant galaxies ($`e.g.,`$ Harris 1991) are explained by the presence of a population of IGCs. In this scenario, putative populations of IGCs provide an “excess” GC population that elevates the observed specific frequencies above a universal value of $`S_N=4`$. Following West $`etal.`$(1995) and Blakeslee $`etal.`$(1997), we assume that the surface density of IGCs is given by

$$\mathrm{\Sigma }_{\mathrm{IGC}}=\mathrm{\Sigma }_0T_X/(1+r^2/r_c^2)$$

where where $`r_c`$ and $`T_X`$ are the cluster’s core radius and X-ray temperature, respectively. To find the constant of proportionality in this relation, we performed a linear fit to the data presented in Figure 13 of Blakeslee, Tonry & Metzger (1997). Taking into account that the area of their detector was 22 arcmin2, the predicted IGC surface density is

$$\mathrm{\Sigma }_{\mathrm{IGC}}=33T_X/(1+r^2/r_c^2)\mathrm{arcmin}^2$$

Using a value of $`T_X=3.9`$ keV for A1185 (Jones & Forman 1999) and recalling that we are sensitive to only the brightest $`16\%`$ of the GC luminosity function, we expect our WFPC2 field to contain $`110`$ IGCs. The predicted number of IGCs is thus of the right order needed to explain the point source excess in A1185.

## 4 Conclusions

Having considered the expected contributions of many different sources, it seems clear that the observed point source excess is best explained by a population of GCs in A1185. Whether these GCs are intergalactic in nature, or whether they are associated with neighboring galaxies, is an open question. The difficulty in deciding between these scenarios lies in the unknown extent to which the cluster mean field has imposed a tidal limitation on the individual cluster galaxies. If tidal truncation has been less effective than assumed in the above analysis, then many, and perhaps most, of the observed point sources observed might be explained by GCs associated with individual galaxies. On the other hand, if the cluster galaxies have been tidally limited as prescribed in Merritt (1984), the IGC hypothesis emerges as the most plausible explanation for the observed excess.

Based on the evidence above, it seems very likely that the point sources observed in our WFCP2 field are bonafide GCs residing in A1185. Nevertheless, this assertion can be tested with deeper imaging of this field. Any population of GCs, whether or not it is intergalactic in nature, should be characterized a near-Gaussian luminosity function. At the distance of A1185, the GC luminosity function is expected to show a turnover at $`I_{TO}27`$, which is well within the imaging capabilities of HST and the Advanced Camera for Surveys. While the luminosity function alone would would not allow us to distinguish between GCs belonging to individuals galaxies and those associated with the cluster as a whole, it may be possible to discriminate between these alternatives by examining the spatial distribution of an expanded sample of sources. Additional clues to the nature of these GCs, including ages and metallicities, will be possible once optical and infrared colors are in hand.

We would like to thank P.B. Stetson for providing PSFs for the F814W filter and an anonymous referee for helpful comments. Support for program GO8164 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Additional support for this work was provided by the National Science Foundation from grant AST-0205960 and through a grant from the Association of Universities for Research in Astronomy, Inc., under NSF cooperative agreement AST-9613615 and by Fundación Andes under project No.C-13442. 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 work was based in part on observations obtained at the W. M. Keck Observatory, which is operated jointly by the California Institute of Technology and the University of California. We are grateful to the W. M. Keck Foundation for their vision and generosity.


format: text

text:
# One dimensional heat conductivity exponent from random collision model

## Abstract

We study numerically the thermal conductivity coefficient $`\kappa `$ as a function of system length $`L`$ for several different quasi one dimensional models: classical gases of hard spheres with both longitudinal and transverse degrees of freedom. We introduce a model that is ergodic and highly chaotic but also conserves energy and momentum, and is very useful because it shows scaling even at small system sizes. We find that $`\kappa L^\alpha `$ over more than two decades, with $`\alpha `$ very close to the analytical prediction of 1/3.

Since the surprising result obtained over thirty years ago that the heat current flowing across a one-dimensional chain of harmonic oscillators with a small temperature difference between the two ends is independent of the length $`L`$ of the chain Lieb , the conductivity of one dimensional systems has been studied analytically Exact ; Dhar ; LLP  and numerically LLP ; Numeric ; Grass ; Grass2 ; Casati  at great length. The standard approach to conductivity would predict that if the temperature gradient $`T`$ in a material is small, the heat current flowing through should be of the form $`j=\kappa T,`$ where $`\kappa `$ is a property of the material. On the other hand, the result for the harmonic oscillator chain, that a small temperature difference $`\mathrm{\Delta }T`$ results in a current $`j\mathrm{\Delta }T`$ that is independent of $`L`$, is equivalent to a conductivity $`\kappa L.`$ In the various models studied thereafter, singular conductivities $`\kappa L^\alpha `$ have been found, with a variety of possible values for $`\alpha .`$ On the other hand, for some one-dimensional models a conventional $`L`$-independent $`\kappa `$ has also been obtained.

It is now believed that the singular conductivity of these models has two possible causes. Firstly, if the model is integrable, as in the case of the harmonic oscillator chain, the system does not equilibrate thermally. The behavior of the conductivity then depends on the details of the system. In fact, it has been shown that by changing the coupling of the oscillator chain to the heat reservoirs at the ends, normally a benign procedure, one can tune the exponent $`\alpha `$ over a range Dhar . Secondly, even if the model is not integrable, if the internal interactions in the system conserve momentum, the conductivity is singular, due to advection of heat in long wavelength modes ONSR ; Grass ; Prozen . If a model is not integrable and does not conserve momentum, $`\kappa `$ should have a well defined limit as $`L`$ diverges dingaling . Analytical studies ONSR ; Grass  predict that for non-integrable momentum conserving systems, $`\alpha `$ should have a universal value of $`1/3.`$

Numerical results for momentum conserving systems have yielded values of $`\alpha `$ ranging from 0.25 to 0.5. Recent studies of one dimensional chains of hard point particles with alternating masses have shown unexpectedly large corrections to scaling even for systems of $`10^4`$ particles, with $`\alpha `$ estimated to be $`0.25`$ Casati  and $`0.33`$ Grass . Similar results have been obtained for chains of Fermi Pasta Ulam chains, with $`\alpha `$ estimated to be 0.37 Grass2 . For the system of hard point particles, the slow convergence to the asymptotic behavior has been justified Grass  by noticing that the system is always ‘close’ to an integrable model. Thus if one considers a chain of particles with equal masses, energy is carried ballistically, the system does not thermalize, and $`\alpha =1.`$ On the other hand, if the ratio of the masses of successive particles in the alternating chain is chosen to be very different from unity, the light particles are almost inconsequential, and the problem again reduces to one of equal mass particles. In Ref. Grass , a mass ratio of 2.62 was found to yield the longest power-law scaling range, from which $`\alpha `$ was estimated.

In this paper, we study how to eliminate residual effects of integrability by considering non-integrable and highly chaotic models. We start with Sinai and Chernov’s pencase model pencase  that was initially introduced by them to study one dimensional hydrodynamics. In this model, hard sphere particles are confined to a long narrow tube. We consider both periodic and hard wall boundary conditions in the transverse direction, and apply heat baths to the two ends. The extent in the transverse direction is taken to be slightly less than twice the diameter of the particles. This ensures that the particles cannot get past each other, but allows a large range of incidence angles at the collisions. Thus the transport of energy along the tube remains quasi one dimensional, with the transverse degree of freedom serving as an additional randomizing effect. This model (without heat reservoirs) has been proved to be ergodic in four dimensions and hyperbolic in three SSz .

The extra degree of freedom in our model limits the system sizes we simulate. For the range of system sizes that we consider, there is still insufficient universality to obtain the asymptotic value for $`\alpha `$ with confidence. Therefore, we have devised a novel 1d model that we call the “random collision” (RC) model, discussed in detail below, that is extremely chaotic, yet still satisfies energy and momentum conservation.

Probably due to the further randomization introduced by our new dynamics, the conductivity fits very well to the form $`\kappa L^\alpha ,`$ with $`\alpha `$ close to 1/3. If the masses of all the particles are taken to be equal, the estimated value of $`\alpha `$ is $`0.29\pm 0.01,`$ while when the particle masses alternate with a mass ratio of 2.62, one obtains $`\alpha `$ to be $`0.335\pm 0.01.`$ The small discrepancy from the theoretical prediction of $`\alpha =1/3,`$ although larger than the error bars, is within the range one might expect from corrections to scaling from the leading irrelevant operators.

We also show results for the autocorrelation function of the energy current, which is related to the conductivity through the Kubo formula Kubo . The autocorrelation function has a zero frequency limit that scales close to $`L^{1/3},`$ in accordance with previous predictions and our numerical results for the conductivity. However, the full correlation function $`C(\omega ,L)`$ shows quite complicated behavior, not easily described by any simple scaling form, possibly indicating the existence of multiple time scales.

One might be concerned whether the analytical derivation of $`\alpha =1/3,`$ which uses the hydrodynamic description of a one-dimensional normal fluid ONSR  is valid for the models we have considered here. Since the transverse direction is small, it should not affect long wavelength singularities in the dynamics. Of greater concern is, with periodic boundary conditions, the existence of the transverse momentum as another conserved quantity. Although this makes the hydrodynamic equations more complicated, and increases their number from three to four, we recall that the analytical calculation relies only on the fact that (without an applied temperature gradient) the system reaches thermal equilibrium, that it satisfies Galilean invariance, and that a finite sound velocity sets a cutoff to the dynamics for any finite system size. None of these conditions is violated by the introduction of the transverse momentum.

Figure 1 shows a log log plot of the conductivity as a function of the length the system $`L,`$ for the pencase model. The results for both periodic and hard wall boundary conditions are shown. All the particles were taken to have the same mass. The diameter of the particles was 0.6, while the cross section of the tube and the average longitudinal separation between centers of neighboring particles were both taken to be 1. The temperatures of the heat reservoirs at the two ends were taken to be $`1.0`$ and $`1.2`$ respectively; it was verified numerically that this temperature difference is sufficiently small for the system to be in the linear response regime. The heat reservoirs at the ends were implemented as follows: whenever an extremal particle collided with the reservoir adjoining it, its velocity was randomized, drawn from the distribution $`P(v_x,v_y)v_x\mathrm{exp}[m(v_x^2+v_y^2)/(2k_BT)],`$ where $`x`$ and $`y`$ are along the longitudinal and transverse direction respectively and $`T`$ is the temperature of the reservoir. (This is the velocity distribution for particles leaking out of a heat reservoir.) At each such collision, the energy exchanged by the system and the reservoir is kept track of, and used to calculate the time average of the energy current at both ends.

In the same figure, the conductivity as a function of system length is also shown for a different choice of model parameters: alternating particle masses, with a mass ratio of 2.62, a tube cross section of 1.14 and a longitudinal interparticle separation of 0.9. (Only hard wall transverse boundary conditions are shown for this case.)

The number of particles in the system ranged from 8 to 1024. This is substantially less than the largest system sizes used when simulating the one dimensional chain of hard point particles. However, introducing the transverse dimension should make the system no longer near integrability, and therefore allow the large $`L`$ limit to be reached quickly. Unfortunately, as seen in Figure 1, this is not the case. There is some curvature in all the plots; more importantly, there is substantial disagreement between the slopes obtained when the particles have equal or alternating masses. Note that the plots all curve downwards, from which one might be tempted to conclude that the asymptotic slope (i.e. the value of $`\alpha `$) would be smaller than obtained from the curves. However, prior experience with the one dimensional system Grass  indicates that it is possible for the curves to turn around at much larger system sizes. Thus one cannot obtain even an upper bound to $`\alpha `$ from the figure, and must conclude that the corrections to scaling are large for this model foot1 .

In order to randomize the dynamics further, enabling faster convergence to the asymptotic scaling form for large $`L,`$ we modify the model above. Firstly, the diameter of the particles is taken to be negligibly small, while keeping the cross-section of the tube as less than twice the diameter. Secondly, the particles are no longer disk shaped, but highly irregular. As a result, when two particles collide with each other, in the center of mass frame they can recoil in any direction, unrelated to the direction of impact. For any collision, we take the recoil angle in the center of mass frame to be a uniform random variable foot2 . The result of both these modifications together is that the transverse coordinate $`y`$ of the particles becomes redundant, and they effectively move along the $`x`$ axis. However, the transverse velocities $`v_y`$ are retained. In this RC model, each particle has both $`v_x`$ and $`v_y,`$ with the latter behaving as an auxiliary variable that is only important in collisions. In any interparticle collision, the total momentum in the $`x`$ and $`y`$ directions and the total energy are conserved. Collisions with the reservoirs at the two ends are still implemented as before. In the transverse direction, periodic boundary conditions correspond to $`v_y`$ for a particle remaining constant between collisions, whereas hard wall boundary conditions allow $`v_y`$ to be reversed. In the latter case, since in the zero cross-section limit any particle undergoes a huge number of collisions with the sidewalls between two collisions with its neighbors, one should change the sign of $`v_y`$ randomly between interparticle collisions.

Figure 2

shows a log log plot of the conductivity as a function of $`L`$ for the RC model. The particles are now points, the average interparticle separation is unity, and the temperatures of the reservoirs are 1 and 1.2 as before. In the linear response regime, the dependence of the energy current on the average interparticle separation and the reservoir temperatures can be found trivially, so the only parameters that can be varied (apart from the number of particles) are the masses of the particles. Only the results for periodic boundary conditions in the transverse direction are shown, as the results for the hard-wall case are very similar. As before, plots are shown for the case when all particles have the same mass and the case when the particle masses alternate with a mass ratio of 2.62. Both parameter choices yield plots that rapidly converge to a power law form, but with slightly different exponents. The estimated exponents are $`\alpha =0.29\pm 0.01`$ and $`\alpha =0.335\pm 0.01`$ respectively for the two cases. Both of these are very close to the analytical prediction of $`\alpha =1/3.`$ Although the difference between the two estimates for $`\alpha `$ is larger than the error bars, we expect that this is due to corrections to scaling from irrelevant operators in a renormalization group analysis; an $`O(1)`$ bare strength for irrelevant operators can produce effective values of $`\alpha `$ that differ from $`1/3`$ by the desired amount. Thus the numerical results of Figure 2 are a strong indication of the validity of the prediction of $`\alpha =1/3.`$

We have also measured the autocorrelation function for the energy current, $`C(t)=J(t+t^{})J(t^{}),`$ where the average is over $`t^{},`$ since this is related to the conductivity by the Kubo formula Kubo . This was done in the center of mass frame, as necessary for the Kubo formula bonetto . Figure 3

shows the Fourier transform of the correlation function, $`C(\omega ),`$ with periodic boundary conditions in the $`x`$ direction (without heat reservoirs), for the case when all the particles have identical mass. Based on the analytical prediction ONSR , we show $`L^{1/3}C(\omega L)`$ for a range of system sizes. The $`\omega 0`$ limit of all the curves appears to be the same, in accordance with the numerical results we have obtained for the conductivity, although the error bars are clearly larger due to the method of measurement. The collapse of the curves to a single plot at low frequencies is less clear; in fact, as shown in Figure 4,

a functional form of $`C(\omega ,L)\omega ^{1/4}`$ seems to work slightly better. The correlation function flattens out from this form at low frequencies, at a cutoff $`\omega _L`$ that appears (for $`L`$ ranging from 64 to 1024, not all shown in the figure) to decrease slightly faster than $`1/L`$ with $`L.`$ The existence of a well defined limit $`L^{1/3}C(\omega 0,L)`$ would require a $`1/L^{4/3}`$ decay for this cutoff scale. In addition to this cutoff, there is a clear resonance peak at $`\omega 1/L,`$ indicating that there are important dynamics on the $`tL`$ timescale.

It is not clear whether the correlation function converges to the simple scaling form for much larger system sizes, or whether there are indeed multiple time scales in the dynamics: $`t_1(L)L`$ and $`t_2(L)L^\beta `$ with $`\beta 4/3.`$ However, earlier work delosRios  on the one dimensional KPZ equation KPZ  and work in two dimensions halpin-healy , has found it difficult to obtain universal critical exponents from numerical simulations, even with large system sizes. The hydrodynamic description used ONSR  to obtain $`\alpha `$ is similar to the KPZ (Burgers) equation, but with three equations instead of one. In fact, $`\alpha `$ was correctly estimated earlier Grass  from the KPZ equation. Thus the slow convergence that we see for $`\alpha `$ for the pencase model — and, to a lesser extent, for the RC model — may be a similar phenomenon to that seen for the KPZ equation. Of course, even the hydrodynamic description for one-dimensional conduction ONSR  does involve two time scales: for propagation and decay of the sound modes.

In this paper, we have introduced a random collision model for studying dynamics of momentum conserving one dimensional systems, in order to obtain the scaling form of the thermal conductivity $`\kappa `$ as a function of system size $`L.`$ This has allowed us to probe the scaling regime better than has been possible with previous models that show a slower approach to the large $`L`$ limit and less robust behavior. Over a wide range of length scales, we find good agreement with the earlier analytical prediction of $`\kappa L^{1/3}.`$

We thank Abhishek Dhar and M.A. Moore for very useful discussions.


text:
# A Search for 6.7 GHz Methanol Masers in OH Megamaser Galaxies at 0.11<𝑧<0.27

## 1 Introduction

A survey for OH megamasers (OHMs) in luminous IR galaxies has recently been completed at the Arecibo Observatory111The Arecibo Observatory is part of the National Astronomy and Ionosphere Center (NAIC), which is operated by Cornell University under a cooperative agreement with the National Science Foundation. (Darling & Giovanelli, 2002a). The survey detected 50 previously unknown OHMs at $`z=0.1`$–0.25, more than doubling the sample of OHMs overall, and increasing the sample at $`z>0.1`$ sevenfold. OHMs are associated with starburst nuclei in major galaxy mergers. The primary factors responsible for producing OHMs appear to be the large-scale density enhancements of molecular gas produced exclusively in galaxy-scale tidal interactions and the strong far-IR radiation field which inverts the OH ground state, providing the maser pump. If this is the case, then the large-scale density enhancements of excited molecular gas combined with the intense IR radiation field can be expected to produce masing in other lines and other molecules. We focus here on the 6.7 GHz ($`5_16_0`$ A+) transition of methanol (CH3OH). This strong masing line is associated with compact HII  regions and OH masing (Menten, 1991). Galactic OH and CH3OH masing regions identify HII  regions, are likely to be commingled, and have similar lifetimes (Caswell et al., 1995; Caswell, Vaile, & Forster, 1995). Typically, 6.7 GHz methanol emission is an order of magnitude stronger than 1665 MHz OH and has a similar width but differs in the detailed line structure (Caswell et al., 1995). Hence, if the Galactic relationship between OH and CH3OH emission in compact HII  regions scales to megamaser luminosities, then one would expect the 6.7 GHz line be easily detectable at cosmological distances. If, however, OH megamasers represent a fundamentally different physical setting from Galactic OH masers then the scaling would not necessarily apply.

Detections of extragalactic methanol masers are limited to the Large Magellanic Cloud (Sinclair et al., 1992; Ellingsen et al., 1994b; Beasley et al., 1996), with isotropic integrated line luminosities similar to observed Galactic methanol masers (roughly 10-5–10-6 $`L_{}`$). There are two surveys for methanol megamasers described in the literature: both select OH and H2O megamaser hosts to search for 6.7 GHz methanol masers, and both obtain only upper limits on maser emission (Ellingsen et al., 1994a; Phillips et al., 1998). Phillips et al. (1998) also surveyed a number of IR-luminous starburst galaxies and made no detections. Does this paucity of detections indicate different mechanisms for OHMs from the OH masers associated with Galactic star formation, or are are methanol masers only associated with a subset of OHMs?

This paper presents the results of a sensitive survey for 6.7 GHz methanol masers in OH megamaser hosts conducted at the Arecibo Observatory. No detections were made, but the survey sets strong limits on the methanol emission for a wide variety of OH megamasers and explores a new region of the redshift-luminosity space available to these extremely rare objects. The sample selection, observation methods, and data reduction techniques are presented in §2. Section 3 presents the limits on 6.7 GHz methanol maser emission in the 25 OHMs observed, and §4 offers a brief interpretation of the implications of this survey for methanol megamasers and the nature of OHMs.

## 2 Data

### 2.1 Sample Selection

A sensitive survey for 6.7 GHz methanol masers associated with OH megamasers would optimally include all known OHMs within the Arecibo declination range ($`0^{}<\delta <37^{}`$) up to the highest redshift OHM known at $`z=0.2655`$ (Baan et al., 1992). Unfortunately, the survey sample must be restricted to $`z0.11`$ for the redshifted 6.7 GHz methanol line to fall within the frequency range of the C-band receiver (3.95–6.05 GHz). No 6–8 GHz receiver is currently available at Arecibo.

We select 28 OHMs at $`z0.11`$ spanning the range of properties found in this diverse population. Line widths range from unresolved to extremely broad (25–850 km s-1), peak flux densities span 2–40 mJy, and 1667:1665 MHz OH hyperfine ratios range from 2.5 to $`>15`$ (the LTE value is 1.8; see Table 1 and Darling & Giovanelli 2002b). The sample includes Seyferts type 1 and 2, LINERs, and pure starburst nuclei and spans the range of (ultra)luminous IR galaxies (\[U\]LIRGs) on the radio-far-IR relation, from radio loud to IR overluminous (Darling & Giovanelli, 2002a, 2003). The sample includes compact and extended OHMs as inferred from variability or quiescence of spectral features (Darling, in preparation). The sample also includes OH lines with smooth profiles, indicating a distributed masing gas, lines with high velocity wings characteristic of molecular flows, and multiply-peaked profiles, indicating a wide range of masing conditions and complicated kinematics. To maximize the chances of detection and to provide the strongest possible constraint on the existence of methanol megamasers, the methanol maser survey sample includes every variety of OHM and OHM host available. This sample is complementary to the Phillips et al. (1998) survey which includes OH and H2O megamasers and starburst galaxies.

### 2.2 Observations

Observations were conducted (remotely) at the Arecibo Observatory with the newly commissioned C-band (3.95–6.05 GHz) receiver on February 26–28, 2002. Each program source was observed in multiple 5-minute on-source, 5-minute off-source sessions with a noise diode calibration at the end of each. The radio frequency interference (RFI) environment at Arecibo generally requires fast sampling for RFI excision, but the effects of RFI at 4–6 GHz are relatively mild and sampling every 6 s was adequate. Spectra were produced using an autocorrelation spectrometer with 9-level sampling in two linear polarizations. We observed a 25 MHz band centered on the redshifted 6668.518 MHz $`5_16_0`$ A+ transition of CH3OH flanked by 25 MHz bands offset by $`\pm 20`$ MHz at the rest frequency and doppler tracked separately to detect any high velocity masing lines which might be present in the program sources. Each 25 MHz bandpass had 1024 channels for a spectral resolution of 1.3 km s-1 and a total velocity coverage of 3500 km s-1 in the rest frame of a galaxy at $`z=0.2`$.

The primary reflector was readjusted to $`2`$ mm rms prior to these C-band observations, which improved the telescope forward gain and the quality of the main beam. The beam size at 5 GHz is 54″$`\times `$63″ (azimuth$`\times `$zenith angle) and the pointing has a total rms error of about 5″. The gain is about 8.5 K Jy-1 at low zenith angles at 5 GHz, but is a strong function of zenith angle, falling to about 6.4 K Jy-1 at a zenith angle of 15 (the maximum zenith angle at Arecibo is 19.7). The gain is also a strong function of frequency, falling by roughly 20$`\%`$ from 5 GHz to 6 GHz. The system temperature is nearly flat at 34 K for zenith angles up to 15 at 5 GHz, but increases by $`15\%`$ from 5 GHz to 6 GHz222See www.naic.edu for current receiver and system characterization information.. The methanol line for the program sources falls in the spectral range 5270–6005 MHz, with the majority of the sources in the range 5470–5955 MHz.

### 2.3 Data Reduction

All data reduction was performed in AIPS++333The AIPS++ (Astronomical Information Processing System) is a product of the AIPS++ Consortium. AIPS++ is freely available for use under the Gnu Public License. Further information may be obtained from http://aips2.nrao.edu. with the NAIC single dish packages. No RFI excision was done. In the three objects where RFI was significant, IRAS 06487+2208, IRAS 12162+1047, and IRAS 14586+1432, it was overwhelming and persistent and could not be removed. This RFI filled the range 5808–5840 MHz (42531–44417 km s-1). On-off pairs and polarizations were averaged to obtain the final spectra. These spectra were Hanning smoothed, and linear baselines were subtracted. The bandpasses are generally excellent: they show no more than a linear slope and are consistent across the 65 MHz of the three combined bandpasses.

All spectra including the flanking fields were inspected after Hanning smoothing, a 6-channel boxcar smoothing to 146 kHz, and gaussian smoothing matched to the OH line width. The boxcar smoothing nearly matches the 156 kHz boxcar-smoothed rms noise quoted by Phillips et al. (1998). This is the noise level listed in column (9) of Table 1. The observations are sensitive to methanol lines from a few to $`2000`$ km s-1 wide and would also identify any high velocity lines up to 1750 km s-1 from the systemic redshift.

## 3 Results

No credible emission or absorption lines were identified in the 25 RFI-free program sources. Table 1 lists the OHMs observed in the redshifted 6668.518 MHz methanol line, including coordinates, optical heliocentric redshifts, OH line peak flux density, full width at half maximum, and isotropic line luminosity. We assume a cosmology with $`H_{}=75`$ km s-1, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ to compute luminosities. To provide the most conservative possible constraint on methanol emission relative to OH emission, we compute the upper limit on the isotropic methanol line luminosity from a 1.5$`\sigma `$ boxcar having width matching the OH line width $`\mathrm{\Delta }`$vOH, where $`\sigma `$ is the rms noise of the 6-channel boxcar smoothed spectrum. We compute the same 6.7 GHz line luminosity limit for the Phillips et al. (1998) OHM sample for confirmed OHMs with published widths (17 objects), enabling direct comparison of the two samples. This OH width-dependent methanol luminosity limit provides a direct measure of the limit on the OH-to-methanol line luminosity under the assumption that the two emission lines scale and are spatially coincident (this is true in Galactic HII  regions; Caswell, Vaile, & Forster 1995). The scaling of methanol masers with OH from Galactic HII  regions up to OH megamaser scales is a straw man hypothesis which has already been discounted by previous methanol surveys (Phillips et al., 1998; Ellingsen et al., 1994a), but this method for constraining methanol emission is a natural way to incorporate the available information about molecular masing regions in the sample galaxies.

The limits on the 6.7 GHz methanol emission in the survey sample are listed in Table 1: $`t_{on}`$ is the on-source integration time in minutes, rms$`_{CH_3OH}`$ is the rms noise in a 6-channel (146 kHz) boxcar averaged spectrum, and $`\mathrm{log}L_{CH_3OH}`$ is the limit on the methanol luminosity matched to the OH line width.

## 4 Discussion

Figure 1 illustrates the methanol-OH megamaser parameter space excluded by this work and the Phillips et al. (1998) survey. Although the Arecibo methanol survey is more sensitive than the Phillips et al. (1998) survey by $`1.5`$ dex, there is significant overlap in the 6.7 GHz line luminosity limits because the Arecibo survey is on average a factor of 3 higher in redshift. The long-dashed line in panels (b) and (d) of Figure 1 indicates the typical order of magnitude methanol to OH intensity ratio found in Galactic compact HII  regions (Caswell et al., 1995). Note that this relationship is derived from the 1665 MHz OH line whereas in OHMs the 1667 MHz line dominates. All of the upper limits on the 6.7 GHz methanol emission fall below this relationship. The discrepancy is most significant for the largest $`L_{OH}`$ objects, indicating a significant difference between the physical conditions producing OH megamasers and those responsible for Galactic OH masers.

There are other indications that OHMs represent a fundamentally different phenomenon from OH masing associated with compact HII  regions. Foremost of these is the 1667:1665 MHz hyperfine ratio: in OHMs the 1667 MHz line dominates while in compact HII  regions the 1665 MHz line dominates. This suggests a lower density environment for OHMs (Pavlakis & Kylafis, 1996). Second, a Galactic star formation rate scaled up to 100–1000 $`M_{}`$ yr-1 does not produce the observed OHM line luminosities found in (U)LIRGs and cannot explain the OHM fraction in these galaxies or the dependence of this fraction on the far-IR luminosity (Darling & Giovanelli, 2002a).

In light of the apparent differences between Galactic OH masers and OHMs, it may not be surprising that no methanol megamasers have been detected in several surveys. But most of the conditions required for methanol masers are still present in OHMs: adequate molecular columns with coherent velocity paths, a strong far-IR radiation bath, and copious continuum photons to amplify. The intensity of 100 and 60 $`\mu `$m radiation is more than adequate to pump methanol, even with low efficiency (van der Walt, Gaylard, & MacLeod, 1995). What component required for methanol masing is missing in OHM galaxies? As suggested by Phillips et al. (1998), it may be the large-scale presence of methanol itself. Methanol in giant molecular clouds or circumnuclear molecular tori may be depleted onto dust grains and restricted to star formation regions where it is liberated by grain heating. In a study of ultracompact HII  regions, Walsh et al. (1997) find that the number of 6.7 GHz methanol maser spots increases with peak maser luminosity, which is interpreted to indicate that maser emission depends on the abundance of methanol rather than on the far-IR radiation field. If methanol megamasers are abundance-limited rather than pump-limited, then we suggest that in a dense, massive starburst environment it is possible for the output from the occasional “normal” methanol maser to be amplified by another cloud with inverted methanol level populations. The required velocity coincidence along with geometrical constraints probably conspire to make this a rare occurrence, but one that could occasionally be found in the nuclei of OHM and luminous IR galaxies. Although no detections were made in this survey, a survey sensitive to $`1`$ km s-1 lines of a larger sample of IR galaxies may detect a few of these rare events.

Numerical modelling can address the absence of 6668 MHz methanol megamasers in OHM galaxies by extending the detailed combined methanol-OH maser models of massive star forming regions by Cragg, Sobolev, & Godfrey (2002) to physical conditions more appropriate for OHMs. Kylafis & Pavlakis (1998), Burdyuzha & Vikulov (1990), and Randell et al. (1995) suggest (with significant disagreement) $`n_{H_2}=10^3`$$`10^5`$ cm-3, $`n_{OH}/n_{H_2}=10^6`$$`10^8`$, $`T_{dust}60`$ K, $`T_{kin}=40`$–60 K, $`L=1`$–10 pc, and $`\mathrm{\Delta }\text{v}=1`$–100 km s-1, where $`L`$ is the masing path length and $`\mathrm{\Delta }\text{v}`$ is the turbulent velocity width in the masing path. The methanol density is a free parameter. Extrapolation of the Cragg et al. (2002) simulations to likely OHM settings suggests that methanol would not be detected by our survey if $`n_{CH_3OH}n_{OH}`$ or $`T_{dust}<110`$ K. While the scaling of $`L`$, $`n_{H_2}`$, and $`\mathrm{\Delta }\text{v}`$ from the values assumed by Cragg et al. ($`L=0.03`$ pc, $`n_{H_2}>10^4`$ cm-3, and $`\mathrm{\Delta }\text{v}=1`$ km s-1) to OHM conditions nearly cancel in the computation of the OH and methanol specific column density ($`N/\mathrm{\Delta }\text{v}`$), $`T_{dust}<110`$ K severely quenches methanol masing. In fact, extrapolation of the Cragg et al. models to $`T_{dust}<100`$ K indicates that no enhancement of methanol abundance or path length can counteract the sensitivity of masing methanol lines to the temperature of the IR radiation field. Global dust temperatures in (U)LIRGs are generally 30–40 K (Dunne et al., 2000), but they may be significantly higher in small-scale star forming regions or in the vicinity of enshrouded AGN. Although low dust temperature is an attractive explanation for the lack of methanol megamasers, confirmation would require measurements of $`T_{dust}`$ in the regions actually producing OH megamasers rather than a global temperature measurement. Note that the 1665 MHz line is always dominant over the 1667 MHz line in the Cragg et al. (2002) models whereas the reverse is observed in OHMs.

## 5 Conclusions

A sample of OH megamasers at $`z0.11`$ spanning a wide range of host galaxy types and OH emission properties has been observed for 6.7 GHz methanol megamaser emission. No lines were identified in the 25 RFI-free OHM galaxies within $`\pm 1750`$ km s-1 of the systemic redshift, and the observations were sensitive to lines from a few to $`2000`$ km s-1 in width. We set limits on methanol emission well below the Galactic methanol-OH intensity ratio of $`10`$, although it should be noted that the Galactic ratio uses the 1665 MHz OH line whereas in OHMs the 1667 MHz line dominates. Combined with 17 OHMs observed by Phillips et al. (1998), this study supports the notion that OH megamasing is a fundamentally different process from the OH masing associated with Galactic compact HII  regions and methanol masers. Maser models suggest that 6668 MHz methanol masers are strongly quenched in OH megamaser galaxies in which $`T_{dust}<100`$ K or $`n_{CH_3OH}n_{OH}`$.

The authors wish to thank the excellent staff of NAIC for observing assistance and support, especially Karen O’Neil, and the AIPS++ consortium for software support, especially Joe McMullin. We thank the anonymous referee for an excellent suggestion and thoughtful comments. This research was supported by NSF grant AST 00-98526 and 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. The National Astronomy and Ionosphere Center is operated by Cornell University under a cooperative agreement with the National Science Foundation.


text:
# Dynamic Nonlinear Effect on Lasing in a Random Medium

## Abstract

We have studied the dynamic effect of nonlinearity on lasing in disordered medium. The third-order nonlinearity not only changes the frequency and size of lasing modes, but also modifies the laser emission intensity and laser pulse width. When the nonlinear response time is comparable to or longer than the lifetime of the lasing mode, the nonlinearity changes the laser output through modifying the size of the lasing mode. When the nonlinear response is much faster than the buildup of the lasing mode, positive nonlinearity always extracts more laser emission from the random medium due to the enhancement of single particle scattering.

PACS numbers: 71.55.Jv,05.45.-a,42.55.-f,42.25.Dd

Over the past decade, there have been considerable interest in nonlinear optical processes in disordered media, e.g. second harmonic generation , phase conjugation . It has been shown that optical nonlinearity can modify light transport in disordered medium.  The coherent backscattering peak exhibits additional features in the presence of third-order nonlinearity.  Temporal fluctuation of scattered waves in nonlinear disordered medium leads to instability of speckle pattern.

Recent studies illustrate that adding gain to a disordered medium leads to lasing in the long-lived eigenmodes.  Despite the modes with long lifetime are preferably amplified, their wavefunction are not modified by the presence of gain. However, nonlinearity can change the eigenmodes of a disordered system. In random lasers nonlinearity is large due to high intensity and resonance enhancement. The laser light intensity is high, especially in the case of pulsed pumping. The nonlinearity at the lasing frequency, which is in resonance of a transition of the gain material, is resonantly enhanced.

In this paper, we study the effect of nonlinearity on lasing in disordered medium. The experimental evidence of the change of lasing modes is presented first. It is followed by numerical simulation of a model system. The algorithm is based on the finite-difference time-domain (FDTD) solution to Maxwell equations coupled with the rate equations of electronic populations. Third-order nonlinearity with finite response time is included. Our results demonstrate that nonlinearity not only modifies the frequency and size of the eigenmodes of a disordered system, but also changes laser emission intensity and laser pulse width of a random laser. The nonlinear response time plays a crucial role in the dynamic nonlinear effect on random lasing.

We used two types of random media in our experiments. One is poly(methyl methacrylate) (PMMA) sheets containing rhodamine 640 (Rh640) perchlorate dye and titanium dioxide (TiO2) microparticles, the other is ZnO nanoparticles with average diameter 100 nm. The dye concentration in PMMA is 10 mM. The TiO2 particles have an average diameter of 0.4 $`\mu `$m. The density of TiO2 particles is $``$ 1.4 $`\times `$ 1012 cm-3. The PMMA sheets are optically excited by the second harmonics of a mode-locked Nd:YAG laser. The pump spot is $``$ 100 $`\mu `$m in diameter at the sample surface. The spectrum of emission from the sample is measured by a 0.5-meter spectrometer with a liquid nitrogen cooled CCD array detector. When pump intensity exceeds the lasing threshold, discrete lasing modes emerge in the emission spectrum. As the pump intensity increases, the gain spectrum of dye molecules shifts towards longer wavelength. As shown in Fig. 1, additional lasing modes appear on the long wavelength side, while some lasing modes on the short wavelength side disappear. However, individual lasing modes shift towards shorter wavelength, as marked in Fig. 1. Typically optical gain pulls the frequencies of lasing modes towards the peak of gain spectrum (frequency pulling effect). In our experiment, the lasing modes move in the opposite direction as the gain spectrum. This suggests the spectral shift of lasing modes is caused by nonlinearity. Nonlinear Kerr effect has been observed in Rh640 solution containing TiO2 particles.  The refractive index $`n`$ changes linearly with light intensity $`I`$: $`n=n_0+n_2I`$, where $`n_0`$ is linear refractive index, $`n_2`$ is the nonlinear Kerr coefficient. The nonlinearity is contributed mainly by dye molecules. The value of $`n_2`$ in our samples is estimated to be $`10^{14}`$ to $`10^{15}`$ m2/W. Above the lasing threshold, the excitation intensity is on the order of $`10^{11}`$ W/m2. Thus the change of refractive index $`\mathrm{\Delta }n`$ ranges from $`10^3`$ to $`10^4`$. This can explain the decrease of lasing wavelengths as the pump intensity increases. The wavelength shift is in the range 0.1 – 1.0 nm, which is on the same order of the experimental shift.

In addition to time-integrated spectra, we measured time-resolved lasing spectra to track spectral shift of lasing modes in time. Figure 2 shows a spectral-temporal image of laser emission from ZnO powder. The ZnO powder is optically excited by 20 ps pulses from a frequency-tripled mode-locked Nd:YAG laser. The pump beam is focused to a 20 $`\mu `$m spot on the sample surface. Laser emission from the powder is dispersed by a 0.3-meter monochromator, then directed into a Hamamatsu streak camera. Individual lasing modes are separated by the monochromator, and their temporal evolutions are recorded by the streak camera. The lasing modes at shorter wavelength are red-shifted with time, while the lasing modes at longer wavelength are blue-shifted with time. The modes in the center do not shift in wavelength. For example, the mode labeled A shifts from 387.76 nm at t = 40 ps to 387.96 nm at t = 130 ps. The mode labeled B shifts from 390.69 nm at t= 40 ps to 390.58 nm at t= 105 ps. The wavelength of mode C remains nearly constant in time.

Temporal shift of lasing frequencies must result from dynamic change of the refractive index of ZnO. It is known that ZnO has a large third-order nonlinearity near its band edge.  The values of its third-order nonlinear coefficient $`\chi ^{(3)}`$ range from $`10^8`$ to $`10^6`$ esu at room temperature. The nonlinear response of ZnO is also fast: its nonlinear response time is measured to be 160 fs near band edge.

To understand the nonlinear effect on random lasing, we performed numerical simulation based on the finite-difference time-domain (FDTD) solution of Maxwell equations.  The model system is one-dimensional simplification of the real experiment. The binary layers are made of dielectric materials with dielectric constants of $`ϵ_1=ϵ_0`$ and $`ϵ_2=9ϵ_0`$, respectively. The thickness of the first layer, which simulates the gain medium, is a random variable: $`a=a_0(1+w_a\alpha )`$, where $`a_0=400`$ nm, $`w_a`$ describes the strength of randomness, $`\alpha `$ is random value in the range \[-0.5,0.5\]. The thickness of the second layer, which simulates the nonlinear medium, is $`b=b_0(1+w_b\beta )`$, where $`b_0=100`$ nm, $`w_b`$ is the strength of randomness, $`\beta `$ is random value in the range \[-0.5,0.5\]. For simplicity, the gain component and the nonlinear component are separated in our model. This separation does not change qualitatively the nonlinear effects that we study. The dynamics of electronic populations in the gain layers is described by four-level rate equations, with parameters similar to those in Ref.  In the nonlinear layers, the third-order nonlinearity is described by the nonlinear polarization

$$P_{NL}(x,t)=ϵ_0\chi ^{(3)}E(x,t)_{\mathrm{}}^{\mathrm{}}g(t\tau )|E(x,\tau )|^2𝑑\tau $$
(1)
$`\chi ^{(3)}`$ is the third-order nonlinear coefficient. The casual response function $`g(t\tau )=(1/\tau _0)\mathrm{exp}[(t\tau )/\tau _0]`$ for $`t\tau `$, and $`g(t\tau )=0`$ for $`t<\tau `$. $`\tau _0`$ is the nonlinear response time. To incorporate the nonlinearity with finite response time into the FDTD algorithm, we introduce a new function

$$G(x,t)_{\mathrm{}}^{\mathrm{}}g(t\tau )|E(x,\tau )|^2𝑑\tau =\frac{1}{\tau _0}_0^te^{(t\tau )/\tau _0}|E(x,\tau )|^2𝑑\tau .$$
(2)
The differential equation for $`G(x,t)`$ can be derived as

$$\frac{dG(x,t)}{dt}=\frac{G(x,t)}{\tau _0}+\frac{|E(x,t)|^2}{\tau _0}$$
(3)
$`E`$ is related to $`D`$ by accounting for both linear and nonlinear polarizations:

$$E(x,t)=\frac{D(x,t)P_L(x,t)P_{NL}(x,t)}{ϵ_0ϵ_{\mathrm{}}}$$
(4)
According to Eq.(4), The effective dielectric constant is

$$ϵ(x,t)=ϵ_{\mathrm{}}+\chi ^{(1)}+\chi ^{(3)}G(x,t).$$
(5)
When light frequency is far from any resonant frequencies, $`\tau _00`$, the nonlinearity becomes instantaneous. $`g(t\tau )=\delta (t\tau )`$, $`G(x,t)=|E(x,t)|^2`$. The nonlinear dielectric constant $`ϵ=ϵ_L+\chi ^{(3)}|E|^2`$, where $`ϵ_L=ϵ_{\mathrm{}}+\chi ^{(1)}`$ is the linear dielectric constant.

To simulate pulsed pumping, the external pumping term in the rate equations represents a Gaussian pulse of width 20 ps: $`P_r(t)=2\times 10^8\mathrm{exp}[\mathrm{ln2}(tt_0)^2/t_0^2]`$ 1/s, where $`t_0`$= 10 ps. We first turn off the nonlinearity by setting $`\chi ^{(3)}`$ = 0 to simulate lasing in linear random medium. The random structure consists of 50 layers with $`w_a=w_b=0.9`$. The radiative transition of the four-level system is centered at $`\lambda _a`$ = 389 nm with 7 nm spectral width. We choose the pumping rate so that only one mode lases. The frequency of the lasing mode does not exhibit any noticeable change during the lasing process. Next we include the nonlinearity. We choose the value of $`\chi ^{(3)}`$ to be $`1.2\times 10^{16}`$ m2/V2, close to the measured value for ZnO at room temperature.  Figure 3 plots temporal evolution of lasing spectra obtained by Fourier transform of the electric field in 2.6 ps intervals. For negative (positive) nonlinearity $`\chi ^{(3)}<0`$ ($`\chi ^{(3)}>0`$), the lasing mode first shifts towards shorter (longer) wavelength as the laser intensity increases, then shifts to longer (shorter) wavelength as the laser intensity decreases. These results can explain qualitatively the experimental observation of lasing frequency shift with time. As shown in Fig. 2, the laser pulses in ZnO powder rise quickly then decay slowly. The rising edge of the laser pulses is too short to detect temporal shift of the lasing frequency. However, the falling edge of the laser pulses is long enough to observe noticeable shift of lasing frequency with time. In the presence of population inversion, $`\chi ^{(3)}`$ is positive (negative) for photon energy below (above) the band gap. For the lasing modes whose frequencies are below (above) ZnO band edge, since $`\chi ^{(3)}>0`$ ($`\chi ^{(3)}<0`$), they shift towards shorter (longer) wavelength in the falling edge as in Fig. 3b (3a). For the lasing modes at the band gap, $`\chi ^{(3)}0`$, thus their frequencies remain nearly constant with time. The wavelength shift in our calculation is $``$ 1–2 $`\AA `$, close to the experimental value.

The spatial size of lasing modes also changes in the presence of third-order nonlinearity. When only one mode lases, the size of the lasing mode can be characterized by the inverse partition ratio $`r(t)=(|E(x,t)|^2𝑑x)^2/|E(x,t)|^4𝑑x`$. Figure 4 plots $`r(t)`$ of one lasing mode (the same one as in Fig. 3) in both linear and nonlinear cases. For this mode, when $`\chi ^{(3)}>0`$ ($`\chi ^{(3)}<0`$), $`r(t)`$ increases (decreases) during the lasing period. We repeat the calculation for lasing modes in many random configurations. However, we do not find any definite relation between the sign of the nonlinearity and the size change of the lasing modes. Nevertheless, in all simulations, we do observe a noticeable alteration in the mode size.

Optical nonlinearity not only changes the frequency and size of lasing modes, but also modifies the intensity and duration of laser pulses. The change of laser output is sensitive to the nonlinear response time $`\tau _0`$. We calculate the laser pulses by fixing the value of $`\chi ^{(3)}`$ and varying $`\tau _0`$. Figure 5 plots the total laser emission energy $`U(t)=(1/2)ϵ_0ϵ|E(x,t)|^2𝑑x`$ for $`\chi ^{(3)}=\pm 1.2\times 10^{16}`$ m2/V2. We first set $`\tau _0`$ equal to the experimental value (160 fs) of passive ZnO near band edge . For positive (negative) $`\chi ^{(3)}`$, the laser pulse becomes weaker (stronger) than that in the linear case. Next we repeat the calculation with smaller $`\tau _0`$ since the nonlinear response time is shortened in the presence of stimulated emission. As shown in Fig. 5, when $`\tau _0`$ = 6.5 fs, the behavior of laser pulse becomes totally different. For positive (negative) $`\chi ^{(3)}`$, the laser pulse becomes stronger (weaker) and longer (shorter) than that in the linear case.

From the calculation of many random lasing modes, we conclude that the effect of nonlinearity on laser pulse intensity and width depends on the relative magnitude of two time scales. One is the nonlinear response time $`\tau _0`$, the other is the lifetime $`\tau _c`$ of the lasing mode in the passive medium. $`\tau _c`$ can be obtained by launching a pulse at the frequency of the lasing mode in the passive medium and observing the temporal decay of the electric field.  When $`\tau _0`$ is comparable to or longer than $`\tau _c`$, the change of laser pulse intensity and width is related to the size change of the lasing mode. If the size of the lasing mode decreases (increases) in the presence of nonlinearity, light confinement gets better (worse). The decrease (increase) of light leakage is equivalent to an increase (decrease) of the quality factor of the random cavity. Hence, lasing lasts longer (shorter), and laser emission is stronger (weaker). When $`\tau _0`$ is much shorter than $`\tau _c`$, the change of laser output depends only on the sign of $`\chi ^{(3)}`$, i.e., positive (negative) nonlinearity always extract more (less) laser emission from the random medium at the same pumping rate. When the nonlinear response is much faster than the buildup of the lasing mode, the lasing mode cannot response fast enough to the nonlinear refractive index change. The phase of scattered light changes quickly due to rapid change of refractive index with intensity. The absence of constant phase relations among light waves scattered by different particles undermines the interference effect. Hence, the effect of single particle scattering becomes dominant over the collective effect of many particle scattering. For $`\chi ^{(3)}>0`$ ($`\chi ^{(3)}<0`$), the refractive index contrast of the binary layers increases (decreases) as the laser intensity increases. Light scattering of a single particle becomes stronger (weaker). The increase (decrease) of scattering strength results in more (less) efficient lasing, i.e., higher (lower) laser intensity and longer (shorter) lasing period.

In conclusion, we have studied the effect of nonlinearity on lasing in disordered medium. The third-order nonlinearity not only changes the frequency and size of lasing modes, but also modifies the laser emission intensity and laser pulse duration. When the nonlinear response time is comparable to or longer than the lifetime of the lasing mode, the nonlinearity changes the laser output through modifying the size of the lasing mode or equivalently the quality factor of the random cavity. When the nonlinear response is much faster than the buildup of the lasing mode, positive nonlinearity always extracts more laser emission from the random medium due to the enhancement of single particle scattering.

We acknowledge Dr. A.L. Burin and S.-H. Chang for stimulating discussions. This work is supported partially by the National Science Foundation under Grant No. DMR-0093949 and the David and Lucille Packard Foundation.