text:
# The Nature of Composite LINER/H II Galaxies, As Revealed from High-Resolution VLA Observations

## 1 Introduction

LINER galaxies are a class of active galaxies characterized by the presence of strong, nuclear, low-ionization emission lines (Heckman 1980). Like Seyfert nuclei, LINERs emit much stronger optical low-ionization forbidden lines compared to H II nuclei, whose line emission is powered by photoionization from young massive stars, but LINERs have a characteristically lower ionization state than Seyferts. In a recent optical survey of nearly 500 nearby galaxies, Ho, Filippenko, & Sargent (1995, 1997a, 1997b) find that $``$20% of all galaxies brighter than $`B_T`$ = 12.5 mag display LINER-type spectra. An additional 13% of the objects show spectra intermediate between those of “pure” LINERs and “pure” H II nuclei. Ho, Filippenko, & Sargent (1993; see also Ho 1996) dubbed this class as “transition objects,” and they hypothesized that these might be composite systems in which the optical signal of a weak active nucleus (the LINER component) has been spatially blended by circumnuclear star-forming regions (the H II component).

Models have shown that photoionization by hot stars are able to reproduce the optical spectra of LINERs (Terlevich & Melnick 1985), especially those of the transition-type variety (Shields 1992; Filippenko & Terlevich 1992). Other researchers have advocated shock-wave heating as a viable excitation mechanism for LINERs (Fosbury et al. 1978; Heckman 1980; Dopita & Sutherland 1995). At the same time, an increasing body of evidence suggests that a significant fraction of LINERs do contain irrefutable AGN characteristics (see reviews by Filippenko 1996 and Ho 1999b). Radio surveys, in particular, have shown that many LINERs exhibit weak nuclear radio emission (Heckman 1980; Sadler, Jenkins, & Kotanyi 1989; Wrobel & Heeschen 1991; Slee et al. 1994; Falcke, Wilson, & Ho 1997; Van Dyk & Ho 1998). As discussed in detail by Ho (1999a), weak nuclear radio emission in early-type (elliptical and S0) galaxies is most likely related to accretion-driven energy production in these objects, but at a very low level compared to powerful radio galaxies, which are also usually hosted by early-type systems.

Radio observations provide an expecially attractive tool to further our understanding of the LINER phenomenon. Since radio frequencies are not plagued by dust obscuration and photoelectric absorption, high-resolution radio observations potentially offer the cleanest and most efficient method to isolate an AGN core. Although radio emission also arises from supernova remnants and H II regions, AGN cores can be identified through their high degree of central concentration, small angular size, flat spectral indices, and high brightness temperatures.

Given that LINERs are so numerous, they may be the most common type of AGN known. Their large space densities could make a major impact on the faint end of the local AGN luminosity function, which in turn has ramifications for many astrophysical issues ranging from the cosmological evolution of AGN to the contribution of AGN to the cosmic X-ray background. We are therefore conducting high-resolution radio observations of selected LINER samples in order to investigate the interrelation between LINERs and classical AGN. For the present paper — the first in a series — we have concentrated on transition-type LINERs. We report on sensitive Very Large Array (VLA)111The VLA is a facility of the National Radio Astronomy Observatory (NRAO) which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. C-array observations of a sample of 25 such transition LINERs at 8.4 GHz. The resulting images, having angular resolution of about 2$``$ $``$ .5, allow us to assess the strength and spatial distribution of the nuclear radio emission. For comparison, we have also observed a small sample of 12 pure H II nuclei.

In the following discussion and in upcoming papers, we adopt a Hubble constant of H0 = 75 km s-1 Mpc-1.

## 2 Sample Selection

The sample of 25 transition LINER/H II galaxies and 12 pure H II nuclei was extracted from the original magnitude-limited Palomar survey of 486 bright northern galaxies (Ho et al. 1995), following the classification criteria outlined in Ho et al. (1997a). The Palomar sample contains 65 transition nuclei (Ho et al. 1997a), from which we chose a subset of 25 that had little or no arcsec-scale resolution radio continuum data published, and that fell within the observation window assigned to us by the VLA. In addition, we selected for comparison a small sample of 12 pure H II nuclei (out of 206 such objects in the full Palomar sample). These were chosen to mimic the Hubble type distribution of the sample of transition objects, again taking into consideration scheduling constraints. We are not aware of any strong biases introduced by our sample selection.

Table 1 lists the 37 sample objects. In addition to the optical positions of the galaxies (columns 3 and 4), this table lists their distance (column 5) and their Hubble type (column 6). The Hubble types and distances are taken from Ho et al. (1997a). The “reference” column (7), together with Table 2, supplies references to earlier radio studies; it is readily apparent that several of the sample objects are well-studied objects, while some are not. We will discuss our observational results in light of these existing data in Section 4. The sample galaxies span a wide range of Hubble types, but the majority are spirals.

– TABLEs 1/2: Sample of transition LINER galaxies and H II nuclei –

## 3 Observations and Data Reduction

The resulting sample of 25 composite LINER/H II and 12 pure H II galaxies was observed with the X-band (8.4 GHz) system of the VLA in two observing sessions, on 1997 August 21 and August 25. The array was in its C-configuration, yielding a typical resolution of 2$``$ $``$ .5. Two bands of 50 MHz each were combined, yielding a total bandwidth of 100 MHz. Snapshot observations combining two scans of 8–9 minutes each, at different hour angles, were obtained in order to improve the shape of the synthesized beam.

Secondary phase and amplitude calibrators were observed before and after each scan. The primary calibrator was 3C 286, with appropriate baseline constraints, and the adopted 8.4 GHz flux densities, as provided by the VLA staff, were 5.22 Jy and 5.20 Jy, at IF1 and IF2, respectively. The calibration uncertainty is dominated by the uncertainty in the absolute flux density of the primary calibrator, which is a few percent. The array performed well: judging from the calibration sources, the antenna phase and amplitude calibration appeared stable to within a few percent.

The radio data are of high quality, and there was no need for extensive flagging of discrepant points. Reduction of the data was performed using standard NRAO AIPS (version 15OCT98) image processing routines. AIPS task IMAGR was employed to Fourier transform the measured visibilities and obtain CLEAN (Högbom 1974; Clark 1980) maps of the sources. Full-resolution maps, having synthesized beams of roughly 2$``$ $``$ .5, and tapered maps, having beams between about 5$``$ $``$  and 12$``$ $``$ , were made. This proved useful to detect weak low surface brightness emission. As a rule of thumb, 1$``$ $``$  corresponds to 50–100 pc for the typical distances involved. Most images reach the theoretical noise level of $`0.04`$ mJy/beam (Perley, Schwab, & Bridle 1989) to within a factor of two. Self-calibration was employed in the analysis of two of the stronger sources (NGC 4552 and NGC 5354), leading to considerable improvements.

## 4 Results

### 4.1 Radio Maps

Noting that a non-detection should be taken to imply no correlated flux density in excess of 0.5 mJy on the arcsec scale, radio contour maps of the detected galaxies are presented in Figures 1–7, at the end of the paper. The radio maps are overlayed on optical images taken from the Digitized Sky Survey (DSS). In a number of cases only the high-resolution map (typical beamsize 2$``$ $``$ .5) or the low-resolution map (typical beamsize 7$``$ $``$  – 10$``$ $``$ ) are shown. Unrelated background radio sources were found in several fields; their number is entirely as expected given the size of the primary beam and the sensitivity of the observations. In cases of large distances from the field center, primary beam corrections were applied. Some of the background sources have been identified using NED (NASA/IPAC Extragalactic Database) and the Cambridge APM (Automatic Plate Measuring) facility. Where possible, we have included them by their name in the tables below. Previously uncatalogued background sources are designated by their sky orientation with respect to the field center for the relevant galaxy. All maps use contouring according to CLEV $`\times `$ (–3, 3, 6, …, MLEV) mJy/beam where MLEV are powers $`3\times 2^\mathrm{n}`$ up to a maximum of 96 and CLEV is the typical noise level on each radio map (see Table 3).

Table 3 lists the various map parameters, with the following column headings: the source name (1), the applied taper (2), the resulting restoring beam (3), the position angle of the beam (4), the image noise level (5), and the relevant figure number (6). The H II nuclei appear at the bottom of the table.

– TABLE 3: Radio map parameters –

Given that good phase stability was obtained, as judged from the VLA phase calibrators, the astrometric accuracy of our overlay procedure must be dominated by the DSS accuracy, which is known to be about 0.6$``$ $``$  (Véron-Cetty & Véron 1996). This is born out by the images of several unresolved and slightly resolved nuclei (NGC 4552, NGC 5354, NGC 5838, NGC 5846) which appear accurate to within a few tenths of an arcsec.

### 4.2 Comments on Individual Galaxies

With reference to the images in the Appendix, we proceed by discussing the results of our VLA imaging in the light of earlier studies of the sample objects. This discussion also deals with the detected background sources and includes important bibliographic information. We will frequently compare our measurements with sample galaxy flux densities obtained within the VLA 1.4 GHz NVSS and FIRST surveys, having typical resolution of 45$``$ $``$  and 5$``$ $``$ , respectively (Condon et al. 1998a; Becker, White, & Helfand 1995), as well as with the Green Bank 20 cm and 6 cm surveys, having typical resolution of $``$700$``$ $``$  and 3$``$ .5, respectively (White & Becker 1992, WB92 hereafter; Becker, White, & Edwards 1991, BWE91 hereafter).

IC 520 : No compact structure was seen in this object, not even on the tapered map. It was also not detected in the NVSS. However, Condon (1987) reports a detection at 1.49 GHz: at 0$``$ .8 resolution, he measured a 3.1 mJy radio source.

NGC 2541 : We did not detect compact structure, not even on the tapered map. This source was also not detected in the FIRST or in the NVSS survey. However, Condon (1987) reports a possible detection at 1.49 GHz, 0$``$ .8 resolution, of a 3.2 mJy source, slightly extended to the SW relative to the optical position of the galaxy.

NGC 2985 : We detect a 1.9 mJy source. At 2$``$ $``$ .5 resolution, the source is slightly extended in the N-S direction, along the orientation of the optical galaxy. The tapered map shows a 2.7 mJy core and about 1 mJy weak low surface brightness emission to the SW. The object was also detected in the NVSS: 44.1 mJy. At 1.49 GHz, 1$``$ .0 resolution, Condon (1987) measured a total flux of 61.9 mJy: the radio emission of NGC 2985 must be extended over tens of arcseconds.

NGC 3593 : On both the untapered and tapered maps we measure about 20 mJy of diffuse, 45$``$ $``$  E-W elongated emission, oriented along the major axis of the host of this H II nucleus. NVSS measured 87 mJy emission, WB92 measured 132 mJy (hence the radio source exceeds one arcminute in size), whereas BWE91 measured 53 mJy. The 1.4 GHz, 5$``$ $``$  resolution map of Condon et al. (1990) shows an E-W oriented source of 63.4 mJy: these and our VLA data imply a steep-spectrum radio source.

We also detect an unrelated 10.6 mJy source, 3$``$  to the N, extended towards the E. As judged from the APM survey, the POSS plates show an as yet unclassified nonstellar source at this position: background source NGC 3593N is most likely a weak distant radio galaxy.

NGC 3627, M 66 : This interacting galaxy belongs to the Leo triplet. The 8.4 GHz maps show a 2$``$  triple source aligned with the inner bar of the optical galaxy. The outer radio components are related to star formation in the disk of the galaxy (Urbanik, Gräve, & Klein 1985). We measure 3.9 mJy integrated emission in the compact central component, while $``$25 mJy are distributed in the NW and SE components. From comparison with low-resolution surveys (NVSS: 324.9 mJy; WB92: 434 mJy; BWE91: 141 mJy) we infer that most of the radio emission of this steep-spectrum radio source is resolved out by our observations. The triple radio structure was also observed by Saikia et al. (1994), at 5 GHz with 2$``$ $``$  resolution. Combining our data with the Saikia et al. (1994) and the Hummel et al. (1987) data, we conclude that the nuclear radio source in NGC3627 must be of variable, flat-spectrum nature. More recently, high resolution (0$``$ $``$ .15) 15 GHz observations detected an unresolved 1.4 mJy core (N. Nagar, private communication).

NGC 3628 : Like NGC 3627, this galaxy also belongs to the Leo triplet. Our 8.4 GHz full-resolution map shows some extended emission to the N and an eastern extension that is roughly aligned with the projected disk of the galaxy. The nuclear region is clearly resolved. The total flux density on our low-resolution map is about 70 mJy, 70$`\%`$ of which is in $``$15$``$ $``$  diffuse emission. The NVSS measured 292 mJy whereas WB92, at lower resolution, detected 402 mJy, indicating resolved large-scale emission. At 1.4 GHz Condon et al. (1990) measured 203 mJy and 205 mJy at 5$``$ $``$  and 1$``$ $``$ .5 resolution, respectively. These and our data imply a steep-spectrum radio source, in agreement with the results from the Effelsberg 100m telescope (Schlickeiser, Werner, & Wielebinski 1984). Carral, Turner, & Ho (1990) observed this object at 15 GHz using the VLA in A-array and measured a total flux of 23 mJy. These authors sampled the innermost part of the source and obtain a $`4`$$``$ $``$  string of a dozen components aligned with the major axis of the galaxy, that they suggest could be star-forming regions in the disk.

NGC 3675 : We detect about 1.5 mJy of low surface brightness emission on our tapered maps. NVSS detected a 48.9 mJy source, whereas FIRST measured 8.04 mJy: the source must be strongly resolved. 1.49 GHz VLA D-array observations (Condon 1987; Gioia & Fabbiano 1987) detected a $``$50 mJy source, oriented N-S, along the major axis of the optical galaxy. Condon, Frayer, & Broderick (1991) report a flux density of 27 mJy at 4.85 GHz and a steep spectral index between 4.8 and 1.4 GHz.

NGC 3681 : We did not detect a radio source on the untapered or tapered maps. In contrast, NVSS measured a weak source of about 4.2 mJy – at much lower resolution.

NGC 3684 : This H II nucleus had no radio source detected, either on the untapered or on the tapered maps, but NVSS detected a 15.9 mJy source. The radio emission of NGC 3684 must be diffuse.

We do, however, detect a partially resolved 18.1 mJy background source, 3$``$  to the SW, which APM identifies with a 19.5 mag stellar-like object. Further identification is still lacking.

NGC 4013 : We detect a moderately resolved 1.1 mJy core and a small extension to the NE, parallel to the projected disk of this edge-on galaxy. There is an additional 3 mJy of weak disk emission. On the basis of the widely different NVSS and FIRST detections (40.5 vs. 11.9 mJy), this disk emission must extend tens of arcseconds. Hummel, Beck, & Dettmar (1991) present 10$``$ $``$  resolution 5 GHz data of NGC 4013 (UGC 6963) which show both the prominent core and the extended disk emission. Recent high resolution (0$``$ $``$ .15) 15 GHz observations did not detect NGC 4013 above a 10$`\sigma `$ limit of 1 mJy (N. Nagar, private communication).

We also detect a slightly resolved 3.8 mJy background source, 2$``$  to the far SE. This source is also clearly seen in the map of Condon (1987). The APM facility reveals a 19.6 mag stellar-like object, which lacks further identification as yet.

NGC 4100 : We detect about 7 mJy of low surface brightness emission associated with this H II nucleus. The tapered map shows a slight extension to the SE along the major axis of the optical galaxy. NVSS detected a 50.3 mJy source and FIRST detected 17.8 mJy, indicating resolution effects.

NGC 4217 : We detect about 22 mJy, distributed along $``$1$``$ .5 of the projected galactic disk of this H II nucleus. From low-resolution observations (NVSS: 123 mJy; WB92: 139 mJy; BWE91: 40 mJy) we conclude that NGC 4217 harbours extended low surface brightness radio emission, which must have a steep radio spectrum.

NGC 4245 : No radio source associated with this H II nucleus was detected on the untapered or tapered maps. Neither NVSS nor FIRST has detected this object at 1.4 GHz.

NGC 4321, M 100 : We detect about 16 mJy of low surface brightness emission in the nuclear region of this well-known grand design spiral in the Virgo cluster. The tapered map shows emission slightly extended to the NW. NVSS detected 87 mJy, while WB92 measured 323 mJy and BWE91 87 mJy. At 1.49 GHz, 0$``$ .9 resolution Condon (1987) measured 180 mJy total flux in a $``$ 3$``$ $`\times `$ 2$``$ region. At 4.9 GHz, 1$``$ $``$ .5 resolution Collison et al. (1994) detect a ring-like structure coincident with our main feature. At 8.5 GHz, 0$``$ $``$ .2 resolution, these authors detect two unresolved sources, having peaks of 0.37 mJy/beam (E) and 0.22 mJy/beam (W), respectively. The eastern radio source is coincident with their 4.9 GHz main component while the one to the W is coincident with the optical nucleus and has a flat spectrum (Collison et al. 1994).

We also detect an unresolved source, 1.5$``$  to the SE, with 1.3 mJy integrated flux, apparently located near the middle of the southern spiral arm of the galaxy. It is clearly visible on the NGC 4321 map. We have identified this object with SN 1979C (Weiler et al. 1981). Collison et al. (1994) measured 2 mJy for this SN 1979C, at 4.9 GHz.

NGC 4369, Mrk 439 : We detect 3.8 mJy of low surface brightness disk emission associated with this H II nucleus. On the basis of widely different NVSS and FIRST detections (24.3 vs. 4.67 mJy), this disk emission must extend over tens of arcsec. The Condon et al. (1990) 1.49 GHz image, at 15$``$ $``$  resolution, shows a 18.3 mJy source with extended emission to the E.

NGC 4405, IC 788 : No compact emission was detected in the untapered or tapered maps. NVSS detected weak (4.5 mJy) emission for this H II galaxy.

We did, however, detect a partially resolved 10.1 mJy source, 2$``$  SW of the NGC 4405 target position. As judged from the APM, the POSS plates show a 19.2 mag stellar-like object at this position, which we subsequently identify as the $`z=1.929`$ QSO, LBQS 1223+1626 (e.g., Hewett, Foltz, & Chaffee 1995).

NGC 4414 : We detect $``$29 mJy of low surface brightness emission, distributed along $``$1$``$ .5, in a structure aligned with the galaxy’s major axis. NVSS detected a 242 mJy source while FIRST detected a double-peaked source with 44 and 64 mJy components. The Condon (1983), 1.465 GHz, 12$``$ $``$ .5 resolution map and the Condon (1987), 1.49 GHz, 1$``$ .0 resolution map show comparable N-S structure, measuring $``$0.22 Jy. Using this last value plus the 78 mJy measured at 4.85 GHz, 15$``$ $``$  resolution, Condon et al. (1991) obtain a steep spectral index.

NGC 4424 : On the tapered map we detect a resolved $``$1 mJy core plus 1.5 mJy weak low surface brightness emission to the E along the projected disk of this candidate merger galaxy (Kenney et al. 1996). NVSS detected a weak 4.5 mJy source. Our full-resolution data show the core of this H II nucleus to be somewhat resolved. This is borne out of our recent high resolution (0$``$ $``$ .25) 8.4 GHz observations which did not detect NGC 4424 above a 3$`\sigma `$ limit of 0.2 mJy.

We also detect an unrelated, slightly resolved 13.3 mJy source, 2$``$ .5 to the SE. The APM identifies it with a 20.2 mag stellar-like object. Further properties of this background source remain as yet unknown.

NGC 4470 : We detect $``$3 mJy of weak low surface brightness emission, extending 30$``$ $``$  along the galaxy major axis. NVSS detected 17.1 mJy associated with this H II nucleus.

We also detect an unrelated, resolved 16.0 mJy radio source, 1$``$ .5 to the NE of NGC 4470, which is not identified on the POSS plate (APM). We identify this object with the radio source TXS 1227+081 (Douglas et al. 1996), which has as yet no optical identification.

NGC 4552, M 89 : We detect about 77 mJy in a strong core plus possibly a jet-like structure to the NE. NVSS detected a 103 mJy source. Since BWE91 measured 64 mJy at 5 GHz, the radio source in NGC 4552 must have a relatively flat radio spectrum, which is in agreement with findings by Condon et al. (1991). The radio source appears to be variable (Wrobel & Heeschen 1984; Ekers, Fanti, & Miley 1983; Sramek 1975a, 1975b; Ekers & Ekers 1973). Our recent subarcsec resolution 8.4 GHz images have confirmed this object’s compactness.

NGC 4643 : No compact radio structure was detected on the untapered or the tapered maps. NVSS also did not detect this source. However, we have detected several very weak sources ($``$1 mJy) on the tapered maps, some of which may be due to weak diffuse disk emission.

NGC 4710 : We detect about 6 mJy of weak disk emission associated with this H II nucleus. NVSS measured a 19.3 mJy source. At 1.49 GHz, Condon et al. (1990) measure 17.2 mJy and 14.7 mJy at 15 and 5$``$ $``$  resolution, respectively. The data imply a resolved steep-spectrum radio source.

We also detect a 4.8 mJy resolved background source, 2.5$``$  to the E of NGC 4710. No indentification could be found on the POSS.

NGC 4713 : We detect $``$1.2 mJy weak low surface brightness disk emission on our tapered maps. NVSS measured a 46.9 mJy source: our observations must have resolved out much of the flux.

NGC 4800 : We detect about 1.2 mJy of weak low surface brightness emission, stretching $`30`$$``$ $``$  along the galaxy major axis. NVSS detected a 23.5 mJy source and FIRST 12.2 mJy. Therefore, our observations are resolving out much of the emission of this H II nucleus.

NGC 4826 , M 64: We detect about 21 mJy of low surface brightness emission. The central region on the full-resolution map (not presented here) shows a double-peaked component, which can be compared to the complex inner triple structure of the 15 and 5 GHz, 2$``$ $``$  resolution maps of Turner & Ho (1994). Combining the data indicates that this is a steep spectrum source. NVSS detected a 101 mJy source and BWE91 detected a 56 mJy source. The 1.49 GHz, 1$``$ .0 resolution data of Condon et al. (1998b) and Condon (1987) show a $``$100 mJy source, as do the data of Gioia & Fabbiano (1987) at 40$``$ $``$  resolution.

NGC 4845 : On our tapered maps we detect a 2 mJy core and some emission to the N, perpendicular to the galaxy major axis, plus 10 mJy of $`20`$$``$ $``$  elongated disk emission. NVSS detected a 43.8 mJy source associated with this H II nucleus. The data of Condon et al. (1990) imply resolution effects on the 10 arcsec scale.

NGC 5012 : We detect some very weak ($``$1 mJy) features near the target phase center. NVSS detected a 31.4 mJy source, but as FIRST did not detect it, we must be dealing with extended low surface brightness emission.

NGC 5354 : We detect an unresolved 11.7 mJy core. NVSS and FIRST detected a 8.4 and 8.0 mJy source, respectively, implying that the NGC 5354 radio source is unresolved at a resolution of 5$``$ $``$ . The 4.85 GHz, 15$``$ $``$  resolution maps of Condon et al. (1991) show 7 mJy emission, implying that the nuclear radio source has an inverted spectrum. Our recent subarcsec resolution 8.4 GHz images have confirmed this object’s compactness.

Paired galaxy NGC 5353 is also detected, 1$``$ .25 to the S, with a 26.7 mJy core. Condon et al. (1991) measure 29 mJy at 4.85 GHz, 15$``$ $``$  resolution, implying that NGC 5353 also hosts a flat spectrum radio source.

NGC 5656 : We detect $``$1 mJy of weak low surface brightness emission, stretching over $``$25$``$ $``$  along the major axis of the galaxy. NVSS detected a 22 mJy source.

NGC 5678 : We detect about 8.5 mJy of low surface brightness emission stretching over $``$1$``$ .2 along the major axis of the galaxy. The tapered map shows what could be a double-peaked source, similar to the 6$``$ $``$  resolution, 1.4 GHz map of Condon (1983). NVSS detected a 112 mJy source while BWE91 measured 68 mJy. The NVSS and BWE91 results, combined with Condon (1983) reporting 109 mJy at 1.4 GHz, imply a steep spectrum source.

NGC 5838 : We have detected a slightly resolved 2.2 mJy source. NVSS detected 3.0 mJy, while Wrobel & Heeschen (1991) report a 2 mJy source (5 GHz, 5$``$ $``$  resolution): the nuclear radio source in NGC 5838 must have a flat radio spectrum. Our recent subarcsec resolution 8.4 GHz images have confirmed this object’s compactness.

There is also an unrelated, unresolved 1.3 mJy source, about 1$``$ .5 to the S of NGC 5838. The POSS plates reveal an 18.4 mag star-like object at that position, for which as yet no redshift is available.

NGC 5846 : We have detected a partially resolved 7.1 mJy source. NVSS measured a 22 mJy source. 1.4 GHz VLA observations by Möllenhoff, Hummel, & Bender (1992) measured an unresolved core of 9 mJy plus 10 mJy of additional diffuse emission: NGC 5846 must possess a compact flat-spectrum core component. Our recent subarcsec resolution 8.4 GHz images have confirmed this object’s compactness.

NGC 5879 : No radio structure was detected on the untapered or the tapered maps. NVSS detected a 21 mJy source, which hence must arise in an extended low surface brightness region.

However, we did detect a 292.1 mJy, slightly resolved background source, about 2$``$  NE of our target. It shows up only on the POSS red plate as a noise-like source. We have identified this source as QSO 1508+5714, at a redshift of 4.301 (Hook et al. 1995). It is in the WENSS (Rengelink et al. 1997) and in BWE91 (282 mJy) as well as in WB92 (149 mJy). Patnaik et al. (1992) have observed the source with the VLA A-array at 8.5 GHz and measured 153 mJy: 1508+5714 must be a variable flat-spectrum quasar.

NGC 5921 : We have detected a very weak 0.5 mJy core and some additional weak extended emission. At 1.49 GHz and 0$``$ .9 resolution, Condon (1987) measured a 20.8 mJy source with an additional 2.8 mJy eastern component which we do not detect at 8.4 GHz. All emission must be of low surface brightness nature.

NGC 6384 : Following the NVSS non-detection, we also do not detect (compact) radio emission from NGC 6384. However, Condon (1987), at 1.49 GHz and 1$``$ .2 resolution, found the radio source, having $``$ 35 mJy, to be very extended.

NGC 6482 : We have not detected this source at 8.4 GHz and neither did NVSS.

We did, however, detect a slightly resolved 23.6 mJy background source, about 3$``$  to the S of NGC 6482. On the full-resolution map this source is slightly extended. The POSS plates show a 20.9 mag stellar object. We have identified this source with the radio source 1749+2302 included in GB6, with 26 mJy flux density.

NGC 6503 : We do not detect radio emission from this galaxy on the full-resolution or on the tapered maps. 1.4 GHz data (NVSS and Condon 1987) indicate a NW-SE extended $``$40 mJy source, which consequently must be diffuse.

We also detect a background source, about 6$``$  to the SW. Due to the large distance from the phase center, which implies a large and uncertain primary beam correction, we cannot assess the flux density of this source. The source appears as a 16.9 mag stellar-like object on the POSS plates. We have identified this object with the BL Lacertae source 1749+701 at redshift 0.77 (Hughes, Aller, & Aller 1992). This object has been extensively observed at other radio-frequencies and is included in the WENSS (Rengelink et al. 1997) as well as BWE91 and WB92. The only other 8.5 GHz measurement of this source comes from Patnaik et al. (1992). They measure 558 mJy with the VLA A-array.

Table 4 summarizes the properties of the background sources, including the previously identified, as well as the new, identifications. Tabulated 8.4 GHz fluxes have been corrected for primary beam attenuation.

– TABLE 4: Field/background sources –

### 4.3 Radio Source Parameters

For each galaxy with a 8.4 GHz detection, we list in Table 5 the radio source parameters, measured both from the full-resolution and the tapered maps. Peak and integrated values (columns 3 and 6), as well as well as sky positions of bidimensional gaussian fits (columns 4 and 5) are tabulated. These values were obtained using the AIPS task IMFIT. As we are often dealing with resolved and/or asymmetric emission, such measurements are often inaccurate and must represent lower limits. In such cases, we have estimated the integrated 8.4 GHz flux density from inspection of the map and the relevant clean component file and we have given it the prefix “$``$” (column 6). Also, given the fact that we have fitted single gaussian to the brightest source components, the accuracy of the radio peak positions is judged to be $`<`$ $``$1$``$ $``$  (see images). For NGC 3627, the quoted peak and integrated value are of the central core. In parentheses we include the total flux of the three components (see Section 4.2).

– TABLE 5: Radio parameters of detected sources –

From a quick comparison with the integrated NVSS 1.4 GHz flux densities, it is found that in many cases the high-resolution observations must have resolved a substantial fraction of the emission. This can be quantified by examination of the spectral index values. Table 6 lists spectral indices $`\alpha _{1.4}^{8.4}`$ ($`S_\nu \nu ^\alpha `$) (column 5), obtained by combining the NVSS flux densities (column 2) and the integrated 8.4 GHz flux densities (column 3) taken from Table 5. Due to the fact that the NVSS beam (45$``$ $``$ ) exceeds the 8.4 GHz beam (7$``$ $``$ –10$``$ $``$ ) by a substantial factor, these spectral index values should be considered as upper limits. It indeed appears that in most of these cases the radio emission was documented to be extended over tens of arcseconds. We also compile calculated NVSS radio luminosities (column 4) for the 37 galaxies under consideration. Again, for NGC 3627 quoted values refer to the compact core except when in parentheses, in which case the values refer to all three source components.

– TABLE 6: summary –

Most sample galaxies are seen to display steep spectral indices, as commonly found for star-forming late-type galaxies (Condon 1992). The following objects have flat ($`\alpha 0.6`$) spectra: NGC 4552, NGC 5354, NGC 5838 and NGC 5846 and the H II nucleus in NGC 4424. We will return to this issue in the next Section.

## 5 Discussion

Our VLA imaging observations have yielded X-band detections of 27 of the 37 sample galaxies. These 27 include NGC 4643 and NGC 5012, which were marginally detected. As for the non-detections, they are equally distributed over the transition LINERs and the H II nuclei. It appears that non-detection at X-band, at a resolution of 2$``$ $``$ .5–12$``$ $``$ , correlates with 1.4 GHz (NVSS) weakness and/or low surface brightness. Table 6 indicates that the NVSS radio luminosity distributions of the subsamples are comparable. The integrated 1.4 GHz radio luminosities imply that the sample sources span the usual luminosity range for nearby galaxies (1026 – 1030.5 erg s-1 Hz-1; see, e.g., Condon 1987, 1992). With the exception of NGC 4424 (see Section 4.2), none of the H II nuclei display compact nuclear emission and/or a flat radio spectrum. This is in strong contrast to the transition LINER galaxies.

Purely based on the radio morphology, the transition LINER objects can be divided into two categories. The first is made up of galaxies displaying extended, steep-spectrum, low surface brightness radio emission, usually tracing the optical isophotes of the host galaxy. The second category refers to objects that, in addition to extended emission, also show compact nuclear radio emission. The first category displays radio morphologies that are consistent with their being due to large-scale star formation; in several cases the radio morphology is seen to trace the H II regions in the host. However, we stress that the presence of a very weak nuclear component in these objects may still be masked by the more dominant radio emission from star-forming regions. The second category, which we will refer to in the following as AGN candidates, is comprised of NGC 3627, NGC 4552, NGC 5354, NGC 5838, and NGC 5846. The last four of these are characterized with a flat radio spectral index between 1.4 and 8.4 GHz (see end of Section 4.3). Judging from high-resolution observations which isolate the nuclear emission, the compact nucleus of NGC 3627 also displays a flat radio spectrum (see Section 4.2). These AGN candidates are on average more luminous than the non-AGN and, as we shall see below, their hosts are early rather than late-type galaxies. While at first sight the H II nucleus NGC 4424 also belongs in this AGN candidate class, the radio image (Fig. 3d) shows it to be resolved. Higher resolution observations recently carried out by us indeed resolve all the nuclear emission in NGC 4424 (see Section 4.2), in contrast to the AGN candidates from the transition LINER sample. As such there is a clear separation between the non-AGN and AGN candidate classes.

Apparently, galaxies with composite LINER and star-formation spectra separate out into objects with and without a clear signature of an AGN, that is to say, a compact, flat-spectrum radio source. The prime question to address, of course, is whether the optical spectra of these classes differ in any way. Deferring the full analysis of this issue to another paper, we here just note that the class of AGN candidates indeed displays somewhat stronger \[N II\] and \[S II\] lines in their optical spectra. As we will demonstrate in the forthcoming paper in this series, the behavior of the so called $`u`$-parameter, which measures the radio/far-IR ratio, also supports the weak AGN classification. However, like in the case of Seyfert galaxies showing a $``$30% incidence rate of compact radio cores (e.g., Norris et al. 1990), there is not a one-to-one correspondence between compact radio emission and nuclear activity. We conclude that the radio properties permit to isolate AGN candidates among the sample of transition LINERs, although the absence of radio cores cannot be used to argue against their being AGN. These conclusions strongly support the hypothesis which resulted from the VLA observations of transition LINER NGC 7331 by Cowan et al. (1994).

Inspection of Table 6 readily shows that these AGN candidates are hosted by rather early-type galaxies (E–Sb). The relevant radio core luminosities L$`_{8.4\mathrm{GHz}}`$ range from $``$ 1026.2 – 1028.5 erg s-1 Hz-1, which is several orders of magnitude less than the core luminosities in FR I or FR II type (Fanaroff & Riley 1974) radio galaxies and quasars, and in the range of the weakest radio cores in Seyfert galaxies (Giuricin et al. 1996). This, then, implies that weak LINER AGN preferentially occur in bulge-dominated hosts, not necessarily just in ellipticals, consistent with statistical results from other lines of evidence (Ho et al. 1997b; Ho 1999b). This is in agreement with the case of Seyfert galaxies, for which the radio core emission was also found to be stronger in early-type hosts (Giuricin et al. 1996). As such, these weak AGN differ from the more powerful FR I radio galaxies, which are commonly associated with elliptical galaxies, not seldom brightest cluster galaxies (Zirbel & Baum 1995). They exceed the very weak nuclear source in the nearby transition LINER NGC 7331 (Cowan et al. 1994) by a factor of $`15`$, and are comparable in strength to the well-known variable radio core of the LINER nucleus in M 81 (e.g., Ho et al. 1999). Hence, it must be concluded that at least some weak LINER AGN reveal themselves by low-power radio cores. This supports the analysis of Ho (1999a), who proposed that the weak nuclear radio sources (radio power 1026 – 1029 erg s-1 Hz-1) in nearby elliptical and S0 galaxies are the low-luminosity counterparts of more powerful AGN.

Finally, it is intriguing that the transition LINER in the elliptical galaxy NGC 6482 remained undetected in our observations (as well as in the NVSS). Given that NGC 6482 is the most distant object in our sample, our non-detection would still allow a $`10^{27}`$ erg s-1 Hz-1 compact radio source to reside in this object.

Full analysis of our sample, including optical emission-line and infrared luminosities, will be presented in forthcoming papers. To address the physical origin behind the findings presented here will be the challenge for future work.

## 6 Conclusions

A sample of composite LINER/H II galaxies and pure H II nuclei was studied using the VLA (C-array) at 8.4 GHz. On the basis of their radio morphological properties, the composite sources can be divided into objects with low surface brightness emission confined mainly to the plane of the galaxy and objects displaying compact nuclear components. The former objects are similar in radio morphology to the H II nuclei in our sample, whereas the latter show morphologies consistent with AGN. All five of the LINER AGN are hosted by bulge-dominated galaxies (Hubble type E–Sb), and four of them show flat spectral indices between 8.4 and 1.4 GHz. In terms of radio luminosity, the present LINER AGN populate the range of low core radio luminosities, which implies that classical AGN of low luminosity exist in a wide range of galaxy types.

###### Acknowledgements.

M. E. F. is supported by grant PRAXIS XXI/BD/15830/98 from the Fundação para a Ciência e Tecnologia, Ministério da Ciência e Tecnologia, Portugal. P. D. B. acknowledges a visitor’s grant from the Space Telescope Science Institute, where a large part of this paper could be written. L. C. H. is partly funded by NASA grant NAG 5-3556, and by NASA grants GO-06837.01-95A and AR-07527.02-96A from the Space Telescope Science Institute (operated by AURA, Inc., under NASA contract NAS5-26555). We want to thank Neil Nagar for providing us with some of his preliminary results. This research was supported in part by the European Commission TMR Programme, Research Network Contract ERBFMRXCT96-0034 “CERES.” We made extensive use of the APM (Automatic Plate Measuring) Facility, run by the Institute of Astronomy in Cambridge, the STScI DSS (Digitized Sky Survey), produced under US government grant NAGW – 2166, and NED (NASA/IPAC Extragalactic Database), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.


format: text

text:
# The Specific Heat of a Trapped Fermi Gas: an Analytical Approach

##

Trapped ultra-cold Fermi gases have received a lot of attention in the last few years. This was partly triggered by the first observations of Bose-Einstein condensation in 1995 . Since then, the field has vastly expanded. The trapping and cooling of fermions is in a much less developed state than that of bosons but it will probably not be long before quantum degeneracy is achieved . (For recent experiments see .)

On the theoretical side, Butts and Rokhsar  have calculated the specific heat numerically and provided analytical results for the spatial and momentum distribution of a non-self-interacting spin polarised Fermi gas at $`T=0`$ in a harmonic trap using the semi-classical (Thomas-Fermi) approximation. Their results are valid for large particle numbers. In addition, Schneider and Wallis  made a numerical study of a similar gas and focused on the effects of small particle numbers. The present paper complements that of Butts and Rokhsar in that we provide an analytical expression for the specific heat of a non-self-interacting spin polarised Fermi gas in the context of a semi-classical approximation. We then compare it to the exact result given by a numerical calculation. As we shall see, the semi-classical approximation produces extremely good results for $`kT\mathrm{}w_{x,y,z}`$, where $`w_{x,y,z}`$ are the frequencies of the trap in the three spatial directions, for particle numbers as low as $`N=1000`$. This high level of accuracy is perhaps unexpected. We will discuss this and an explanation will be given. The characteristic temperature for quantum degeneracy, the Fermi temperature, is given by $`kT_F=\mathrm{}w(6N)^{\frac{1}{3}}`$ for an isotropic trap. The condition we impose is therefore valid for a quantum degenerate regime. It is also worth noting that this condition holds in current experiments with trapped fermions .

The effects of interactions has been dealt with in several works  and in particular the possibility of a BCS transition has been studied in the case of 6Li . Considering the gas to be non-interacting is a very good approximation for a dilute neutral atomic gas. The only case where interactions could be important is trapped 6Li with two hyper-fine spin states . 6Li has an anomalously large s-wave scattering length. At least two different states are needed for s-wave scattering as it is forbidden for particles in the same state. However, the results of ref.  are for $`kT<<\mathrm{}w`$, lying outside the range of temperatures concerned in the present work.

Consider a gas of fermions in a harmonic potential. Neglecting interactions between the particles, the energy levels of each particle are

$`E_{n_1n_2n_3}`$ $`=`$ $`(n_1+{\displaystyle \frac{1}{2}})\mathrm{}w_1+(n_2+{\displaystyle \frac{1}{2}})\mathrm{}w_2+(n_3+{\displaystyle \frac{1}{2}})\mathrm{}w_3,`$ (1)
$`n_i`$ $`=`$ $`0,1,2,\mathrm{},`$ (2)
where $`w_{1,2,3}`$ are the frequencies of the trap.

We will use grand–canonical statistics throughout. This makes the calculations much easier and is justified as the difference between canonical and grand-canonical results is minute for the particle numbers we consider .

The number of particles in the system is

$$N=\underset{n}{}[e^{\beta (E_n\mu )}+1]^1.$$
(3)
The internal energy of the system is given by

$$U=\underset{n}{}E_n[e^{\beta (E_n\mu )}+1]^1.$$
(4)
$`\mu `$ is the chemical potential and the sum in (2) and (3) is over all particle states. Here it is convenient for simplicity of the expression, to make the substitutions

$`w_1`$ $`=`$ $`w,`$ (5)
$`w_2`$ $`=`$ $`\lambda w,`$ (6)
$`w_3`$ $`=`$ $`\lambda ^{}w,`$ (7)
$`x`$ $`=`$ $`\beta \mathrm{}w,`$ (8)
$`\mu `$ $`=`$ $`\mathrm{}w\left({\displaystyle \frac{1+\lambda +\lambda ^{}}{2}}ϵ\right),`$ (9)
where $`x`$, $`ϵ`$, $`\lambda `$ and $`\lambda ^{}`$ are newly defined dimensionless variables. The internal energy becomes

$$U=\frac{1}{2}\mathrm{}w(1+\lambda +\lambda ^{})N+\frac{1}{2}\mathrm{}wu,$$
(10)
where

$`u`$ $`=`$ $`2{\displaystyle \underset{n_1=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n_2=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n_3=0}{\overset{\mathrm{}}{}}}k[e^{x(k+ϵ)}+1]^1,`$ (11)
$`k`$ $`=`$ $`n_1+\lambda n_2+\lambda ^{}n_3.`$ (12)
Only the second term contributes to the specific heat since $`N`$ is held fixed. We have

$$C=\left(\frac{U}{T}\right)_{N,w}=\frac{1}{2}\mathrm{}w\left(\frac{u}{T}\right)_{N,w}.$$
(13)
Note that the partial derivative at constant particle number and trap frequency is the natural way of defining the specific heat for a gas in the situation concerned. After a change of variables from $`T`$ to $`x`$ this becomes

$$\frac{C}{k_B}=\frac{1}{2}x^2\left(\frac{u}{x}\right)_{N,w},$$
(14)
where $`k_B`$ is Boltzmann’s constant.

The simplest approximation consists in replacing the triple sums in (6) by a triple integral. Thus,

$$uv=2_0^+\mathrm{}_0^+\mathrm{}_0^+\mathrm{}k[e^{x(k+ϵ)}+1]^1𝑑y_1𝑑y_2𝑑y_3$$
(15)
and $`k`$ is now the real function $`k=y_1+\lambda y_2+\lambda ^{}y_3`$. Changing the integrating variable to $`k`$ we have

$$v=\frac{1}{\lambda \lambda ^{}}_0^+\mathrm{}k^3[e^{x(k+ϵ)}+1]^1𝑑k.$$
(16)
Note that for the isotropic case we could have transformed the triple sum in (6) into a single sum and only then convert it to an integral, yielding

$$v=\frac{1}{\lambda \lambda ^{}}_0^+\mathrm{}(k^3+3k^2+2k)[e^{x(k+ϵ)}+1]^1𝑑k,$$
(17)
which, at first sight, might be thought to be a more accurate expression. For the anisotropic case, though not trivial, we could likewise have a density of states expanded to three terms . We will coment on this shortly.

To solve the integral in (10) we Taylor expand the integrand in powers of $`e^{x(k+ϵ)}`$ and integrate each resulting term. Here we have to consider two different cases: $`ϵ>0`$ and $`ϵ<0`$. For $`ϵ>0`$ this yields

$$v=\frac{6}{\lambda \lambda ^{}x^4}\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n+1}\frac{e^{nxϵ}}{n^4}.$$
(18)
From (8) and (12) and retaining the fact that $`v`$ is an approximated version of $`u`$, we have

$$\frac{C}{k_B}\frac{12}{\lambda \lambda ^{}x^3}A_4+\frac{3}{\lambda \lambda ^{}x^2}\left(\frac{(ϵx)}{x}\right)_{N,w}A_3,$$
(19)
where $`A_i`$ denotes the sum

$$A_i=\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n+1}\frac{e^{nx|ϵ|}}{n^i}.$$
(20)
The use of the absolute value of $`ϵ`$ is redundant in the case we are considering ($`ϵ>0`$) but it is useful when we consider the $`ϵ<0`$ case.

To obtain an expression for the partial derivative in (13) we use the equality

$$\left(\frac{N}{x}\right)_{N,w}=0,$$
(21)
with $`N`$, the number of particles, being

$$N\frac{A_3}{\lambda \lambda ^{}x^3}.$$
(22)
Expression (16) is obtained using the same process as was used to obtain (13). We then have

$$\left(\frac{(ϵx)}{x}\right)_{N,w}\frac{3}{x}\frac{A_3}{A_2}$$
(23)
and

$$\frac{C}{k_B}\frac{3}{\lambda \lambda ^{}x^3}\left(4A_43\frac{A_3^2}{A_2}\right).$$
(24)
Finally we note that the sums $`A_i`$ can be put in terms of polylogarithms , yielding

$$\lambda \lambda ^{}\frac{C}{k_B}3x^3\left[4Li_4(e^{xϵ})\frac{1}{2}Li_4(e^{2xϵ})3\frac{[Li_3(e^{xϵ})\frac{1}{4}Li_3(e^{2xϵ})]^2}{Li_2(e^{xϵ})\frac{1}{2}Li_2(e^{2xϵ})}\right].$$
(25)
For $`ϵ<0`$, the procedure is similar but slightly more complicated. The integral of expression (10) must be divided in two parts, according to whether the exponential is greater or smaller than $`1`$, as the Taylor expansion of the integrand is different for the two cases. The final result for $`ϵ<0`$ is

$$\lambda \lambda ^{}\frac{C}{k_B}12A_4x^3+21\zeta (4)x^3+6\zeta (2)ϵ^2x^1+\frac{1}{2}ϵ^4x+\frac{(A_3\zeta (2)ϵx\frac{1}{6}ϵ^3x^3)^2}{A_2\zeta (2)\frac{1}{2}ϵ^2x^2},$$
(26)
where $`\zeta `$ is the Riemann zeta function. Putting the sums $`A_i`$ in terms of polylogarithms we have

$`\lambda \lambda ^{}{\displaystyle \frac{C}{k_B}}`$ $``$ $`[12Li_4(e^{xϵ})+{\displaystyle \frac{3}{2}}Li_4(e^{2xϵ})+21\zeta (4)]x^3+6\zeta (2)ϵ^2x^1+{\displaystyle \frac{1}{2}}ϵ^4x+`$ (28)
$`{\displaystyle \frac{[Li_3(e^{xϵ})\frac{1}{4}Li_3(e^{2xϵ})\zeta (2)ϵx\frac{1}{6}ϵ^3x^3]^2}{Li_2(e^{xϵ})\frac{1}{2}Li_2(e^{2xϵ})\zeta (2)\frac{1}{2}ϵ^2x^2}}.`$
Note that these analytical expressions are not an approximation over the Thomas-Fermi result. They are the exact Thomas-Fermi result, which of course, is in itself an approximation. We computed the specific heat numerically using an exact expression composed of several sums, which is easily obtained inserting (6) in (8) and using (15) to get an expression for $`((ϵx)/x)_{N,w}`$. This is then compared to the values for the specific heat given by (18) and (20). The accuracy of the Thomas-Fermi approximation is extremely good for $`x`$ of order 1 or less, as can be seen in figures 1 and 2. As an example, for $`N=1000`$ and $`x_1=x_2=x_3<2.5`$ the error in $`C`$ is always less than $`1\%`$ and if we take $`x_i<1`$, the error is always less than $`0.1\%`$. For larger values of $`x`$ the approximation rapidly deteriorates. The results for a higher number of particles are even more accurate, as expected from this kind of approximation .

Given that for this approximation we have replaced triple sums by triple integrals, such a level of accuracy is quite surprising. Note that the expression inside the sums in (6) falls off exponentially. Therefore the first few terms are very important and expression (9) should differ significantly from (6). And indeed it does. The reason our results are so accurate is that the value of $`ϵ`$ used in (18) and (20) is not the same as the one used in (6). For our Thomas-Fermi approximation we took $`ϵ`$ from expression (16), ie, an approximate $`ϵ`$. So, the insertion of an approximate $`ϵ`$ in an approximate expression leads to a fortunate cancelation of errors. We have tested this by inserting the exact value of $`ϵ`$ into (18) and (20). This leads to much worse results. Also, we have tried to further improve our approximation by considering the first two terms of the cubic polynomial in (11) instead of only the first term and doing the same to the number of particles, obtaining an additional term in (16) as well as additional terms in (18) and (20). To our surprise, this only increases the error. We can only assume that the very high level of accuracy of the Thomas-Fermi approximation is a happy coincidence. In fact, if we take the exact value of $`ϵ`$ and input it in expressions (18) and (20) and in the equivalent expressions obtained using $`k^3+3k^2`$ in (11) instead of only $`k^3`$, then the latter yield the most accurate result, though still not a very good one.

###### Acknowledgements.

J.N. acknowledges financial support by the Portuguese Foundation for Science and Technology under grant PraxisXXI/ BD/5660/95.


text:
# A comparison between broad histogram and multicanonical methods

## I Introduction

Development of tools for optimization of computer simulations is a field of great interest and activity. Cluster updating algorithms , probability reweighting procedures  and, more recently, methods  that obtain directly the spectral degeneracy $`g(E)`$ are a few examples of very successful approaches (for reviews of these methods see, for instance,  and references therein). The Multicanonical (MUCA) and Broad Histogram (BHM) methods belong to the former category. The Entropic Sampling Method (ESM)  was proven to be an equivalent formulation of MUCA . From the knowledge of $`g(E)`$, these methods allow us to obtain any thermodynamical quantity of interest for the system under study, as the canonical average

$$Q_T=\frac{\underset{E}{}g(E)Q(E)\mathrm{exp}(E/T)}{_Eg(E)\mathrm{exp}(E/T)}$$
(1)
of some macroscopic quantity $`Q`$ (magnetization, density, correlations, etc). Both sums run over all allowed energies (for continuous spectra, they must be replaced by integrals), $`T`$ is the fixed temperature, and the Boltzmann constant was set to unity. The degeneracy function $`g(E)`$ simply counts the number of states with energy $`E`$ (which, must be interpreted as a density within a narrow window $`\mathrm{d}E`$, for continuous spectra). Also,

$$Q(E)=\frac{\underset{S[E]}{}Q_S}{g(E)}$$
(2)
is the microcanonical, fixed-$`E`$ average of $`Q`$. The sum runs uniformly over all states $`S`$ with energy $`E`$ (or, again, within the small window $`\mathrm{d}E`$, for continuous spectra). Note that neither $`g(E)`$ nor $`Q(E)`$ depend on the particular environment the system is actually interacting with, for instance the canonical heat bath represented by the exponential Boltzmann factors in eq. (1). Thus, these methods go far beyond the canonical ensemble: once $`g(E)`$ and $`Q(E)`$ were already determined for a given system, one can study its behavior under different environments, or different ensembles, using the same $`g(E)`$ and $`Q(E)`$. Accordingly, as an additional advantage, only one computer run is enough to evaluate the quantities of interest in a large range of temperatures and other parameters.

MUCA was introduced in 1991 by Berg and Neuhaus. The basic idea of the method is to sample microconfigurations performing a biased random walk (RW) in the configuration space leading to another unbiased random walk (i.e. uniform distribution) along the energy axis. Thus, the visiting probability for each energy level $`E`$ is inversely proportional to $`g(E)`$. By tuning the acceptance probability of movements in order to get a uniform distribution of visits along the energy axis, one is able to get $`g(E)`$ at the end. MUCA has been proven to be very useful and efficient to obtain results in many different problems, such as first order phase transitions, confinement/deconfinement phase transition in SU(3) gauge theory, relaxation paths , conformal studied of peptides, helix-coil transition and protein folding, evolutionary problems  and to study phase equilibrium in binary lipid bilayer  (for reviews of the method see ).

BHM was introduced three years ago by de Oliveira et al . It is based on an exact relation between the spectral degeneracy, $`g(E)`$, and the microcanonical averages of some special macroscopic quantities. A remarkable feature of BHM is its generality, since these macroscopic quantities can be averaged by different procedures . Because it is not restricted to a rule like a biased random walk, more adequate dynamics can be adopted for each different application. BHM has been applied to a variety of magnetic systems such as the 2D and 3D Ising Models (also with external fields, next nearest neighbor interactions), 2D and 3D XY and Heisenberg Models, Ising Spin Glass  with accurate results in a very efficient way. Theoretically, the method can be applied to any statistical system . Another distinguishing feature of BHM is that its numerical results do not rely on the number of visits $`H(E)`$ to each energy level, a quantity which is updated by one unit for each new visited state. On the contrary, BHM is based on microcanonical averages of macroscopic quantities. Each visited state contributes with a macroscopic upgrade for the measured quantities. Thus the numerical accuracy is much better than that obtained within all other methods, better yet for larger and larger systems, considering the same computer effort.

Besides the practical points described above, the most important feature of BHM is the following conceptual one. All reweighting methods  depend on the final distribution of visits $`H(E)`$ along the energy axis. Histogram methods  adopt a canonical dynamics, getting $`H_{T_0}(E)`$ for some fixed temperature $`T_0`$; a new distribution $`H_T(E)`$ is then analytically inferred for another (not simulated) temperature $`T`$. Following the same reasoning, one can also obtain $`g(E)`$ . Multicanonical approaches , on the other hand, tune appropriate dynamics in order to obtain a flat distribution $`H(E)`$. In both cases, the actually implemented transition probabilities from energy level $`E`$ to another value $`E^{}`$ are crucial. In other words, in both cases the results depend on the comparison of $`H(E)`$ with the neighboring $`H(E^{})`$. All those reweighting methods are, thus, extremely sensitive to the particular dynamic rule adopted during the computer run, i.e. to the prescribed transition probabilities from $`E`$ to $`E^{}`$.

BHM is not a reweighting method. It does not perform any reweighting on the distribution of visits $`H(E)`$. It needs only the knowledge of the microcanonical, fixed-$`E`$ averages of some particular macroscopic quantities. The possible transitions from the current energy $`E`$ to other values are exactly taken into account within these quantities (see section III), instead of performing a numerical measurement of the corresponding probabilities during the computer run. Thus, the only important role played by the actually implemented dynamic rule is to provide a good statistics within each energy level, separately: the relative weight of $`H(E)`$ as compared to $`H(E^{})`$, i.e. the relative visitation frequency for different energy levels, is completely irrelevant. One can even decide to sample more states inside a particularly important region of the energy axis (near the critical point, for instance) , instead of a flat distribution. In short, any dynamic rule can be adopted within BHM, the only constraint is to sample with uniform probability the various states belonging to the same energy level, not the relative probabilities concerning different energies.

In this paper we have used the formulation of MUCA given by Lee  called Entropic Sampling, which, from now on, we call ESM. We present a comparison between ESM/MUCA and BHM, focusing on both accuracy and the use of CPU time. We choose to start our study with the same example used in the original ESM paper by Lee : the $`4\times 4\times 4`$ simple cubic Ising model, for which the exact energy spectrum is known . Our results show that BHM gives more accurate results than ESM/MUCA with the same number of Monte Carlo steps. Despite the fact that one Monte Carlo step in BHM takes more CPU time, in measuring further macroscopic averages, the overall CPU time is smaller for the same accuracy, at least for this model. Also, BHM can be applied to larger lattices without the problems faced by ESM/MUCA, as we show in our simulations of the same model in a $`32\times 32`$ lattice.

This paper is structured as follows: in section II and III we review the implementation of ESM/MUCA and BHM (including a detailed description of the distinct dynamics adopted in this work). In section IV our numerical tests are presented and discussed. Conclusions are in section V.

## II The Multicanonical Method

The idea of the Multicanonical method is to obtain the spectral degeneracy of a given system using a biased RW in the configuration space . The transition probability between two states $`X_{\mathrm{old}}`$ and $`X_{\mathrm{new}}`$ is given by

$$\tau (X_{\mathrm{old}},X_{\mathrm{new}})=e^{[S(E(X_{\mathrm{new}}))S(E(X_{\mathrm{old}}))]}=\frac{g(E_{\mathrm{old}})}{g(E_{\mathrm{new}})}$$
(3)
where $`S(E(X))=\mathrm{ln}g(E)`$ is the entropy and $`E(X)`$ is the energy of state $`X`$. The transitional probability (3) satisfies a detailed balance equation and leads to a distribution of probabilities where a state is sampled with probability $`1/g(E)`$. The successive visitations along the energy axis follow a uniform distribution. However, $`g(E)`$ is not known, a priori. In order to obtain $`g(E)`$, Lee  proposes the following algorithm:

Step 1: Start with $`S(E)=0`$ for all states;

Step 2: Perform a few unbiased RW steps in the configuration space and store $`T(E)`$, the number of tossed movements to each energy $`E`$ (in this stage, $`T(E)=H(E)`$ because all movements are accepted);

Step 3: Update $`S(E)`$ according to

$$S(E)=\{\begin{array}{cc}S(E)+\mathrm{ln}T(E)\hfill & \text{, if }T(E)0\hfill \\ S(E)\hfill & \text{, otherwise.}\hfill \end{array}$$
(4)
Step 4: Perform a much longer MC run using the transitional probability given by eq. (3), storing $`T(E)`$.

Step 5: Repeat 3 and 4. This is considered one iteration.

This implementation is known to be quite sensitive to the lengths of the MC runs in steps 2 and 4. In section IV we study, in two examples, how the accuracy depends on the total number and size of each iteration.

## III The Broad Histogram Method

BHM  enables us to directly calculate the energy spectrum $`g(E)`$, without any need for a particular choice of the dynamics to be used . Many distinct dynamic rules could be used, and indeed some were already tested .

Within BHM, the energy degeneracy is calculated through the following steps (alternatively, other quantities could replace $`E`$):

Step 1: Choice of a reversible protocol of allowed movements in the state space. Reversible means simply that for each allowed movement $`X_{\mathrm{old}}X_{\mathrm{new}}`$ the back movement $`X_{\mathrm{new}}X_{\mathrm{old}}`$ is also allowed. It is important to note that these movements are virtual, since they are not actually performed. In this work we take the flips of one single spin as the protocol of movements;

Step 2: For a configuration $`X`$, to compute $`N(X,\mathrm{\Delta }E)`$, the number of possible movements that change the energy $`E(X)`$ by a given amount $`\mathrm{\Delta }E`$. Therefore $`g(E)N(E,\mathrm{\Delta }E)`$ is the total number of movements between energy levels $`E`$ and $`E+\mathrm{\Delta }E`$, according to the definition (2) of microcanonical averages;

Step 3: Since the total number of possible movements from level $`E+\mathrm{\Delta }E`$ to level $`E`$ is equal to the total number of possible movements from level $`E`$ to level $`E+\mathrm{\Delta }E`$ (step 1, above), we can write down the equation

$$g(E)N(E,\mathrm{\Delta }E)=g(E+\mathrm{\Delta }E)N(E+\mathrm{\Delta }E,\mathrm{\Delta }E).$$
(5)
The relation above is exact for any statistical model and energy spectrum . It can be rewritten as

$$\mathrm{ln}g(E+\mathrm{\Delta }E)\mathrm{ln}g(E)=\mathrm{ln}\frac{N(E,\mathrm{\Delta }E)}{N(E+\mathrm{\Delta }E,\mathrm{\Delta }E)}$$
(6)
This equation can be easily solved for all values of $`E`$, after $`N(E,\mathrm{\Delta }E)`$ is obtained by any procedure, determining $`g(E)`$ along the whole energy axis. In cases where $`\mathrm{\Delta }E`$ can assume more than one value, eq. (5) becomes an overdetermined system of equations. However, the spectral degeneracy can be obtained without need of solving all equations simultaneously, since the spectral degeneracy is the same for all values of $`\mathrm{\Delta }E`$.

The exact Broad Histogram relation (5) is independent of the procedure by which $`N(E,\mathrm{\Delta }E)`$ is obtained . Therefore, virtually any procedure can be adopted in this task, for instance, an unbiased energy RW , a microcanonical simulation , or a mixture of both . Even the juxtaposition of histograms obtained through canonical simulations at different temperatures, a completely unphysical procedure, could be used, as in some of the results presented in , and explicitly used in  where BHM is re-formulated under a transition matrix  approach. Here, we are going to introduce an alternative procedure, referred as Entropic Sampling-based Dynamics for BHM (ESDYN, hereafter). First, one implements ESM, as described in section II, in order to perform the visitation in the configuration space. Additionally, for each visited state $`X`$, we store the values of $`N(X,\mathrm{\Delta }E)`$ cumulatively into $`E`$histograms. Therefore, at the end we have two choices for the determination of the spectral degeneracy, either by using the entropy accumulated through $`T(E)`$ (that is the traditional ESM/MUCA) or by using the accumulated $`N(E,\mathrm{\Delta }E)`$ and the BHM relation, eq.(5). Because of this special implementation, we can guarantee that exactly the same states are visited for both methods. Hence, the eventual difference in the performances reported in this work must be credited to the methods themselves and not to purely statistical factors.

The dynamic rule originally used in order to test BHM  prescribes an acceptance probability $`p=N(E+\mathrm{\Delta }E,\mathrm{\Delta }E)/N(E,\mathrm{\Delta }E)`$. Both the numerator and the denominator are read from the currently accumulated histograms, and thus $`p`$ varies during the simulation. Wang  has proposed a new approach: instead of using the dynamically updated values of $`N(E,\mathrm{\Delta }E)`$ as in , the transitional probabilities follow a previously obtained (from a canonical simulation, for example) distribution $`N_{\mathrm{fixed}}(E,\mathrm{\Delta }E)`$, kept fixed during the simulation. An alternative and simpler derivation of Wang’s dynamics can be done by using the BHM relation (5) itself. ¿From this, we readily obtain that Wang’s dynamics is the same as using the transitional probability $`p=g_{\mathrm{fixed}}(E)/g_{\mathrm{fixed}}(E+\mathrm{\Delta }E)`$ with approximated values $`g_{\mathrm{fixed}}(E)`$ kept fixed during the simulation. We refer to this dynamics as approximated Wang’s dynamics, since $`N_{\mathrm{fixed}}(E,\mathrm{\Delta }E)`$ (or $`g_{\mathrm{fixed}}(E)`$) is actually only an approximation of the real $`N(E,\mathrm{\Delta }E)`$ (or $`g(E)`$). We will use the results obtained by ESM or BHM with ESDYN as inputs to the approximate Wang’s dynamics. These dynamics will be called AWANG1 and AWANG2, respectively.

For comparison, we also implemented a dynamics that uses the ESM probabilities taken from the exact values of $`g(E)`$. It is worth noticing that, as pointed out in the previous paragraph, this is equivalent to Wang’s proposal with exact values for $`N_{\mathrm{fixed}}(E,\mathrm{\Delta }E)`$. We refer to this dynamics as WANG. Its purpose is only to test the accuracies of the other approaches, once one does not know the exact $`g(E)`$ a priori, in real implementations.

## IV Numerical Tests

We start our comparison with the smallest system, since it is also present in the ESM original paper by Lee . The partition function for the $`4\times 4\times 4`$ simple cubic Ising model is exactly known . It is given by the polynomial function

$$Z(\beta )=\underset{n=0}{\overset{96}{}}C(n)u^n$$
(7)
where $`u=\mathrm{exp}(4\beta )`$, $`\beta =1/T`$, $`T`$ is the temperature. The energy spectrum is written in terms of the coefficients $`C(n)`$ as $`g(E)=g(2n)=C(n)`$ for $`n=0`$ to $`n=96`$. For this model, only the first $`49`$ coefficients are necessary by symmetry. The other $`48`$ coefficients are mirror images of the first ones. Our results will be expressed in terms of $`S(E)=\mathrm{ln}g(E)`$.

In order to compare the entropies as obtained by both ESM and BHM, we normalize the entropy such that $`g(96)=1`$, i.e. $`S(96)=0`$. This point corresponds to the center of the energy spectrum, or, alternatively, to infinite temperature. Of course, the error relative to the exact value vanishes for $`E=96`$. In fig. (1), we compare the normalized (with respect to its exact value) entropy as function of $`E`$ obtained by ESM and BHM (with four different dynamics). In AWANG1 and AWANG2 dynamics we use as $`N_{\mathrm{fixed}}(E,\mathrm{\Delta }E)`$ the results obtained by ESM and BHM with ESDYN, respectively. BHM with the entropic sampling dynamics or using the exact relation proposed by Wang present errors within the same order of magnitude, while pure ESM gives the worst results, as clearly seen in the inset. It is also clear that AWANG1 and AWANG2 results are worse than BHM with ESDYN.

The methods can be better compared by the ratio of their relative errors rather than their absolute values. We define the relative errors in the entropy, for a given energy $`E`$, as

$$ϵ(E)=\left|\frac{S(E)S(E)_{\mathrm{exact}}}{S(E)_{\mathrm{exact}}}\right|$$
(8)
In fig. (2) we show the ratio between the BHM relative errors (obtained by ESDYN, AWANG1 and AWANG2 dynamics) and the ESM ones. The inset shows the ratio between the errors from BHM with WANG and the ESM ones. In its worst performance, the error obtained by BHM with ESDYN is roughly one third of the one obtained by ESM. BHM with Wang’s exact dynamics gives slightly better results once it uses the exact $`g(E)`$ as input in order to get $`g(E)`$ as output. BHM with ESDYN gives on average results $`11`$ times more accurate than ESM, for this number of iterations. However, BHM with the AWANG1 and AWANG2 dynamics give relatively poorer results. Therefore, we have shown that the ES dynamics is a powerful approach (among other possibilities ) to obtain with great accuracy the microcanonical averages $`N(E,\mathrm{\Delta }E)`$ needed by BHM. Let us stress that the results for $`g(E)`$ obtained with AWANG2 are worse than the input they use, namely the values of $`g(E)`$ obtained as output of BHM with ESDYN. In order to obtain good results with Wang’s dynamics, we need a pretty good estimative of $`g(E)`$ as input, and not only some crude estimation as claimed in .

In fig. (3), we show the time evolution of the mean error as a function of the number of Monte Carlo steps (MCS). One iteration in ESM corresponds to perform many RW steps, according to eq. (3), storing the number of visits at each energy and, after this fixed number of RW steps, to update the entropy, according to eq. (4). As we pointed out before, the ESM performance is quite sensitive to the choice of the number of RW steps before each entropy update. For a fixed number of MCS ($`10^6`$), we plot the time evolution for both ESM and BHM, but with different number of iterations. As one can see in fig.(3), the best results for ESM correspond to the smaller number of iterations since, in this case, more RW steps are performed and, consequently, we have better statistics for the determination of the entropy.

It is worth noticing that the ESM/MUCA errors seem to stabilize after a few iterations. We believe that this effect occurs because the number of visits to each energy is not sufficient to provide a good statistics in determining the entropy. The only way to improve the accuracy in the spectral degeneracy is to perform more RW steps between each entropy update (it does not increase the CPU time if the total number of MCS is kept constant). Conversely, for BHM the error decreases monotonically because macroscopic quantities are stored, leading to a good statistics even if some states are not frequently visited. After the very first steps, the errors decay as $`t^{1/2}`$, as expected. It is also remarkable the accuracy when using BHM: we reach the same accuracy of the best performance of ESM/MUCA, by performing roughly 30 times less MCS.

Of course, accuracy is not the only important factor concerning the efficiency of a computational method. BHM with ES dynamics has one additional step compared to the traditional implementation of ESM, namely the storage of the macroscopic quantities $`N(X,\mathrm{\Delta }E)`$. So, we need to know the cost of this additional step. In table I, we show the CPU time (in seconds) spent in both implementations on a 433 MHz DEC Alpha, in order to obtain the results shown in fig.(1). For direct comparison, we also present the CPU time relative to the ESM CPU time. All dynamics tested within BHM takes roughly the same CPU time. As one can see, BHM uses twice more CPU time than ESM. However, for the same number of steps, and the best strategy for ESM, BHM is at least 10 times more accurate. If we consider accuracy and CPU use, we argue that BHM is more efficient than ESM.

Up to now, we have tested both approaches in a very small lattice. Nevertheless, we can show that BHM is even more efficient for larger systems, as expected due to the macroscopic character of the quantities $`N(X,\mathrm{\Delta }E)`$. For $`L^d`$ Ising spins on a lattice, for instance, even restricting the allowed movements only to single-spin flips, the total number of movements starting from $`X`$ is just $`L^d`$. Thus, being a finite fraction of them, $`N(X,\mathrm{\Delta }E)`$ is a macroscopic quantity (this is true along the whole energy axis, except at the ground state where $`g(E)`$ presents a macroscopic jump relative to the neighboring energy levels). In fig. (4) we present results for the time evolution of the mean error for a $`32\times 32`$ square lattice Ising model (a lattice that is 16 times larger than the one in the previous results). The exact solution for this system is also known . Here, the accuracy of BHM is two orders of magnitude higher than that of ESM. Again the errors decrease as $`t^{1/2}`$ for BHM with ESDYN, while the ESM mean errors seem to stabilize. In summary, BHM with ESDYN can obtain accurate results for lattices much larger than the ones considered as limit for ESM.

## V Conclusions

The Multicanonical  and the Broad Histogram  methods are completely distinct from canonical Monte Carlo methods, once they focus on the determination of the energy spectrum degeneracy $`g(E)`$. This quantity is independent of thermodynamic concepts and depends only on the particular system under study. It does not depend on the interactions of the system with the environment. Thus, once one has determined $`g(E)`$, the effects of different environments can be studied using always the same data for $`g(E)`$. Different temperatures, for instance, can be studied without need of a new computer run for each $`T`$.

The goal of this paper is to discuss the conceptual differences between Multicanonical and Broad Histogram frameworks, and to compare both methods concerning accuracy and speed. We obtained the energy spectrum of the Ising model in $`4\times 4\times 4`$ and $`32\times 32`$ lattices, for which the exact results are known and, therefore, provide a good basis for comparison. Our findings show that a combination of the Broad Histogram method and the Entropic Sampling random walk dynamics (BHM with ESDYN) gives very accurate results and, in addition, it needs much less Monte Carlo steps to obtain the same accuracy as the pure Entropic Sampling method. This advantage of the Broad Histogram method grows with the system size, and it does not present the limitations of the Multicanonical or Entropic Sampling methods concerning large systems.

The reason for the better performance is that the BHM  uses the microcanonical averages $`N(E,\mathrm{\Delta }E)`$  of the macroscopic quantity $`N(X,\mathrm{\Delta }E)`$ \- the number of potential movements which could be done starting from the current state $`X`$, leading to an energy variation of $`\mathrm{\Delta }E`$. In this way, each new visited state contributes with a macroscopic value for the averages one measures during the computer simulation. Being macroscopic quantities, the larger the system, the more accurate are the results for these averages. Conversely, Histogram  and Multicanonical  approaches rely exclusively on $`H(E)`$, the number of visits to each energy. Therefore, each new averaging state contributes with only one more count to the averages being measured, i.e., $`H(E)H(E)+1`$, independent of the system size.

Under a conceptual point of view, BHM is also completely distinct from the other methods which are based on the final distribution of visits $`H(E)`$. Alternatively, it is based on the determination of microcanonical, fixed-$`E`$ averages $`N(E,\mathrm{\Delta }E)`$ , concerning each energy level separately. Thus, the relative frequency of visitation between distinct energy levels, which is sensitive to the particular dynamic rule one adopts, i.e. the comparison between $`H(E)`$ and $`H(E^{})`$, does not matter. The only requirement for the dynamics is to provide a uniform sampling probability for the states belonging to the same energy level. The transition probabilities from one level to the others are irrelevant.

## VI Acknowledgments

This work has been partially supported by Brazilian agencies CAPES, CNPq and FAPERJ. The authors acknowledge D. C. Marcucci and J. S. Sá Martins for suggestions, discussions and critical readings of the manuscript.

FIGURE CAPTIONS

* Entropies (normalized by its exact values) for the $`4\times 4\times 4`$ Ising ferromagnet, obtained by ESM/MUCA and BHM for $`100`$ iterations of $`10^6`$ Monte Carlo steps, each. The inset shows a detailed view of the first fourth of the whole spectrum. ESM/MUCA gives the largest errors. AWANG1 and AWANG2 dynamics also give worse results than BHM with both ES or WANG dynamics.
* Ratio between the relative errors for BHM and the ESM/MUCA ones. The horizontal lines are the mean relative errors. BHM is on average more than ten times more accurate than ESM/MUCA, for this number of iterations. The exact Wang’s dynamics is, as expected, slightly more accurate than ESDYN (as shown in the inset) once it uses the exact $`g(E)`$ as input in order to get $`g(E)`$ as output. AWANG1 and AWANG2 give the worst results among the four distinct dynamics presently used in order to test BHM (see text).
* Time evolution of the mean error for BHM with ESDYN (a) and ESM/MUCA (b). In both cases, exactly the same averaging states are visited: thus, the differences are due to the methods themselves, not to statistics. As described in the text, the entropic sampling dynamics is quite sensitive to the number of Monte Carlo steps between each iteration, that correspond to an update of the entropy. Here, we present the results for $`N=10,100,1000`$ iterations, that means $`10^5,10^4`$ and $`10^3`$ Monte Carlo steps, respectively, between each iteration. The more Monte Carlo steps, the smaller is the error for the entropic sampling. Conversely BHM does not seem to depend on the computational strategy adopted. Moreover, the errors seem to stabilize after some steps, in the ESM case. On the contrary, for BHM, to get a better accuracy is a simple matter of increasing the computer time, once the errors decay as $`t^{1/2}`$. The results are averaged over ten realizations.
* Time evolution of the mean error for BHM with ESDYN (a) and ESM/MUCA (b) for the $`32\times 32`$ square lattice Ising Model. Now we consider $`10^7`$ MC steps. Again we present the results for $`N=100,1000`$ iterations, that means $`10^5`$ and $`10^4`$ Monte Carlo steps, respectively, between each iteration. The results are averaged over ten realizations. The ratio between the accuracies of BHM and ESM/MUCA is even higher for large systems.


text:
# Opening up extra dimensions at ultra-large scales

## Abstract

The standard picture of viable higher-dimensional theories is that direct manifestations of extra dimensions occur at short distances only, whereas long-distance physics is described by effective four-dimensional theories. We show that this is not necessarily true in models with infinite extra dimensions. As an example, we consider a five-dimensional scenario with three 3-branes in which gravity is five-dimensional both at short and very long distance scales, with conventional four-dimensional gravity operating at intermediate length scales. A phenomenologically acceptable range of validity of four-dimensional gravity extending from microscopic to cosmological scales is obtained without strong fine-tuning of parameters.

preprint: DTP/00/13, hep-th/0002072

Our world is often thought to have more than four fundamental space-time dimensions. This point of view is strongly supported by string/M-theory, and higher dimensional theories are currently being developed in various directions. The standard lore is that in phenomenologically viable models, extra dimensions open up at short distances only, whereas above a certain length scale, and all the way up to infinite distances, physics is described by effective four-dimensional theories.

In this letter we show that the latter picture is not universal: there are models in which extra dimensions open up both at short and very long distances. Namely, gravity may become higher-dimensional in both of these extremes. At intermediate length scales, physics is described by conventional four-dimensional laws, so the phenomenology of these models is still acceptable.

The starting point for our discussion is the observation  that extra dimensions may be infinite, with the usual matter residing on a 3-brane embedded in higher-dimensional space. In the original Randall–Sundrum (RS) model , gravity is effectively four-dimensional at large scales due to the existence of a graviton bound state localized near the 3-brane in five dimensions. We point out that in other five-dimensional models with an infinite extra dimension (in particular, in the model introduced in Ref.), localization of the graviton may be incomplete, and it is this property that leads to the restoration of the five-dimensional form of gravity at very long distances. The latter feature is manifest both in the case of static sources, where the four-dimensional Newton’s gravity law changes to its five-dimensional counterpart above a certain length scale, and in dynamics of gravitational waves which disappear into the fifth dimension thus violating four-dimensional energy conservation, after travelling long enough distance. It is likely that these phenomena are not peculiar to five-dimensional theories and occur in a number of models with more than one extra dimension.

It has been found recently  that exotic large-distance effects may also appear in models with compact extra dimensions, due to the possible presence of very light Kaluza–Klein states. This interesting scenario is, however, considerably different from ours: in compact models, extra dimensions show up at large distances somewhat indirectly, through the spectrum of the corresponding four-dimensional effective theory: four-dimensional energy is conserved, and so on. In our case, physics at ultra-large distances in intrinsically five-dimensional.

As our concrete example, let us consider a five-dimensional model of Ref. . The model contains one brane with tension $`\sigma >0`$ and two branes with equal tensions $`\sigma /2`$ placed at equal distances to the right and to the left of the positive tension brane in the fifth direction. The two negative tension branes are introduced for simplicity, to have $`Z_2`$ symmetry, $`zz`$, in analogy to RS model (hereafter $`z`$ denotes the fifth coordinate). We assume that conventional matter resides on the positive tension brane, and in what follows we will be interested in gravitational interactions of this matter. The $`Z_2`$ symmetry enables us to consider explicitly only the region to the right of the positive tension brane.

The physical set-up (for $`z>0`$) is as follows: The bulk cosmological constant between the branes, $`\mathrm{\Lambda }`$, is negative as in the RS model, however, in contrast to that model, is zero to the right of the negative tension brane. With appropriately tuned $`\mathrm{\Lambda }`$, there exists a solution to the five-dimensional Einstein equations for which both positive and negative tension branes are at rest at $`z=0`$ and $`z=z_c`$ respectively, $`z_c`$ being an arbitrary constant. The metric of this solution is

$$ds^2=a^2(z)\eta _{\mu \nu }dx^\mu dx^\nu dz^2$$
(1)
where

$$a(z)=\{\begin{array}{cc}e^{kz}\hfill & 0<z<z_c\hfill \\ e^{kz_c}a_{}\hfill & z>z_c\hfill \end{array}$$
(2)
The constant $`k`$ is related to $`\sigma `$ and $`\mathrm{\Lambda }`$ as follows: $`\sigma =\frac{3k}{4\pi G_5}`$, $`\mathrm{\Lambda }=\sigma k`$, where $`G_5`$ is the five-dimensional Newton constant. The four-dimensional hypersurfaces $`z=const.`$ are flat, the five-dimensional space-time is flat to the right of the negative-tension brane and anti-de Sitter between the branes. The spacetime to the left of the positive tension brane is of course a mirror image of this set-up.

This background has two length scales, $`k^1`$ and $`\zeta _ck^1e^{kz_c}`$. We will consider the case of large enough $`z_c`$, in which the two scales are well separated, $`\zeta _ck^1`$. We will see that gravity in this model is effectively four-dimensional at distances $`r`$ belonging to the interval $`k^1r\zeta _c(k\zeta _c)^2`$, and is five-dimensional both at short distances, $`rk^1`$ (this situation is exactly the same as in RS model), and at long distances, $`r\zeta _c(k\zeta _c)^2`$. In the latter régime of very long distances the five-dimensional gravitational constant gets effectively renormalized and no longer coincides with $`G_5`$.

To find the gravity law experienced by matter residing on the positive tension brane, let us study gravitational perturbations about the background metric (1). We will work in the Gaussian Normal (GN) gauge, $`g_{zz}=1`$, $`g_{z\mu }=0`$. The linearized theory is described by the metric

$$ds^2=a^2(z)\eta _{\mu \nu }dx^\mu dx^\nu +h_{\mu \nu }(x,z)dx^\mu dx^\nu dz^2$$
(3)
There are two types of linearized excitations in this model. One of them is a four-dimensional scalar — the radion — that corresponds to the relative motion of the branes (see, e.g., Refs. ) The wave function of the radion is localized between the branes, and its interactions with matter on the positive tension brane result in a scalar force of the Brans–Dicke type. If the distance between the branes is not stabilized by one or another mechanism, the radion has zero four-dimensional mass and gives rise to the usual $`1/r`$ potential in an effective four-dimensional theory. In the particular model under discussion, the radion excitation has been studied in Ref. . It is conceivable that the distance between the branes may be stabilized (cf. Ref. ); in which case the interactions due to the radion switch off at large distances.

In this letter we are interested in other types of excitation that leave the branes at rest, namely the five-dimensional gravitons. When the radion is disregarded, there exists a frame which is GN with respect to both branes simultaneously. In this frame, the transverse-traceless gauge can be chosen, $`h_\mu ^\mu =0`$, $`h_{\mu ,\nu }^\nu =0`$, and the linearized Einstein equations take one and the same simple form for all components of $`h_{\mu \nu }`$,

$$\{\begin{array}{cc}h^{\prime \prime }4k^2h\frac{1}{a^2}\mathrm{}^{(4)}h=0\hfill & 0<z<z_c\hfill \\ h^{\prime \prime }\frac{1}{a_{}^2}\mathrm{}^{(4)}h=0\hfill & z>z_c\hfill \end{array}$$
(4)
The Israel junction conditions on the branes are

$$\{\begin{array}{cc}h^{}+2kh=0\hfill & \text{at }z=0\hfill \\ \left[h^{}\right]2kh=0\hfill & \text{at }z=z_c\hfill \end{array}$$
(5)
where $`\left[h^{}\right]`$ is the discontinuity of the $`z`$-derivative of the metric perturbation at $`z_c`$, and four-dimensional indices are omitted. A general perturbation is a superposition of modes, $`h=\psi (z)e^{ip_\mu x^\mu }`$ with $`p^2=m^2`$, where $`\psi `$ obeys the following set of equations in the bulk,

$$\{\begin{array}{cc}\psi ^{\prime \prime }4k^2\psi +\frac{m^2}{a^2}\psi =0\hfill & 0<z<z_c\hfill \\ \psi ^{\prime \prime }+\frac{m^2}{a_{}^2}\psi =0\hfill & z>z_c\hfill \end{array}$$
(6)
with the junction conditions (5) (replacing $`h`$ by $`\psi `$). It is straightforward to check that there are no negative modes, i.e., normalizable solutions to these equations with $`m^2<0`$. There are no normalizable solutions with $`m^20`$ either, so the spectrum is continuous, beginning at $`m^2=0`$. To write the modes explicitly, it is convenient to introduce a new coordinate between the branes, $`\zeta =\frac{1}{k}e^{kz}`$, in terms of which the background metric is conformally flat. Then the modes have the following form,

$$\psi _m=\{\begin{array}{cc}C_m\left[N_1\left(\frac{m}{k}\right)J_2(m\zeta )J_1\left(\frac{m}{k}\right)N_2(m\zeta )\right]\hfill & 0<z<z_c\hfill \\ A_m\mathrm{cos}\left(\frac{m}{a_{}}(zz_c)\right)+B_m\mathrm{sin}\left(\frac{m}{a_{}}(zz_c)\right)\hfill & z>z_c\hfill \end{array}$$
(7)
where $`N`$ and $`J`$ are the Bessel functions. The constants $`A_m,B_m`$ and $`C_m`$ obey two relations due to the junction conditions at the negative tension brane. Explicitly,

$`A_m`$ $`=`$ $`C_m\left[N_1\left({\displaystyle \frac{m}{k}}\right)J_2(m\zeta _c)J_1\left({\displaystyle \frac{m}{k}}\right)N_2(m\zeta _c)\right]`$ (9)
$`B_m`$ $`=`$ $`C_m\left[N_1\left({\displaystyle \frac{m}{k}}\right)J_1(m\zeta _c)J_1\left({\displaystyle \frac{m}{k}}\right)N_1(m\zeta _c)\right]`$ (10)
The remaining overall constant $`C_m`$ is obtained from the normalization condition. The latter is determined by the explicit form of Eq. (6) and reads

$$\psi _m^{}(z)\psi _m^{}(z)\frac{dz}{a^2(z)}=\delta (mm^{})$$
(11)
One makes use of the asymptotic behaviour of $`\psi _m`$ at $`z\mathrm{}`$ and finds

$$\frac{\pi }{a_{}}(|A_m|^2+|B_m|^2)=1$$
(12)
which fixes $`C_m`$ from (9) and (10).

It is instructive to consider two limiting cases. At $`m\zeta _c1`$ we obtain by making use of the asymptotics of the Bessel functions,

$$C_m^2=\frac{m}{2k}\left[J_1^2\left(\frac{m}{k}\right)+N_1^2\left(\frac{m}{k}\right)\right]^1$$
(13)
which coincides, as one might expect, with the normalization factor for the massive modes in RS model. In the opposite case $`m\zeta _c1`$ (notice that this automatically implies $`m/k1`$), the expansion of the Bessel functions in Eqs. (9) and (10) yields

$$C_m^2=\frac{\pi }{(k\zeta _c)^3}\left(1+\frac{4}{(m\zeta _c)^2(k\zeta _c)^4}\right)^1$$
(14)
It is now straightforward to calculate the static gravitational potential between two unit masses placed on the positive-tension brane at a distance $`r`$ from each other. This potential is generated by the exchange of the massive modes (cf. Refs.  – recall that the zero mode here is not localized)

$$V(r)=G_5_0^{\mathrm{}}𝑑m\frac{e^{mr}}{r}\psi _m^2(z=0)$$
(15)
It is convenient to divide this integral into two parts,

$$V(r)=G_5_0^{\zeta _c^1}𝑑m\frac{e^{mr}}{r}\psi _m^2(0)+G_5_{\zeta _c^1}^{\mathrm{}}𝑑m\frac{e^{mr}}{r}\psi _m^2(0)$$
(16)
At $`rk^1`$, the second term in Eq. (16) is small and it is similar to the contribution of the continuum modes to the gravitational potential in RS model. It gives short distance corrections to Newton’s law,

$$\mathrm{\Delta }V_{short}(r)\frac{G_5}{kr^3}=\frac{G_N}{r}\frac{1}{k^2r^2}$$
(17)
where $`G_N=G_5k`$ is the four-dimensional Newton constant.

Of greater interest is the first term in Eq. (16) which dominates at $`rk^1`$. Substituting the normalization factor (14) into this term, we find

$$V(r)=\frac{G_5}{r}_0^{\zeta _c^1}𝑑m\frac{\pi }{(k\zeta _c)^3}\left(1+\frac{4}{(m\zeta _c)^2(k\zeta _c)^4}\right)^1\frac{4k^2}{\pi ^2m^2}e^{mr}$$
(18)
This integral is always saturated at $`m\stackrel{<}{_{}}r_c^1\zeta _c^1`$, where

$$r_c=\zeta _c(k\zeta _c)^2k^1e^{3kz_c}$$
(19)
Therefore, we can extend the integration to infinity and obtain

$`V(r)`$ $`=`$ $`{\displaystyle \frac{G_N}{r}}{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle \frac{e^{\frac{2r}{r_c}x}}{x^2+1}}`$ (20)
$`=`$ $`{\displaystyle \frac{2G_N}{\pi r}}\left[\text{ci}(2r/r_c)\mathrm{sin}(2r/r_c)\text{si}(2r/r_c)\mathrm{cos}(2r/r_c)\right]`$ (21)
where $`x=mr_c/2`$, and $`\text{ci/si}(t)=_t^{\mathrm{}}\frac{\mathrm{cos}/\mathrm{sin}(u)}{u}𝑑u`$ are the sine and cosine integrals. We see that $`V(r)`$ behaves in a peculiar way. At $`rr_c`$, the exponential factor in Eq. (20) can be set equal to one and the four-dimensional Newton law is restored, $`V(r)=G_N/r`$. Hence, at intermediate distances, $`k^1rr_c`$, the collection of continuous modes with $`mr_c^1`$ has the same effect as the graviton bound state in RS model. However, in the opposite case, $`rr_c`$, we find

$$V(r)=\frac{G_Nr_c}{\pi r^2}$$
(22)
which has the form of “Newton’s law” of five-dimensional gravity with a renormalized gravitational constant.

It is clear from Eq. (20) that at intermediate distances, $`k^1rr_c`$, the four-dimensional Newtonian potential obtains not only short distance corrections, Eq. (17), but also long distance ones, $`V(r)=G_N/r+\mathrm{\Delta }V_{short}(r)+\mathrm{\Delta }V_{long}(r)`$. The long distance corrections are suppressed by $`r/r_c`$, the leading term being

$$\mathrm{\Delta }V_{long}(r)=\frac{G_N}{r}\frac{r}{r_c}\frac{4}{\pi }\left(\mathrm{ln}\frac{2r}{r_c}+𝐂1\right)$$
(23)
where $`𝐂`$ is the Euler constant. The two types of corrections, Eqs. (17) and (23), are comparable at roughly $`r\zeta _c`$. At larger $`r`$, deviations from the four-dimensional Newton law are predominantly due to the long-distance effects.

In our scenario, the approximate four-dimensional gravity law is valid over a finite range of distances. Without strong fine-tuning however, this range is large, as required by phenomenology. Indeed, the exponential factor in Eq. (19) leads to a very large $`r_c`$ even for microscopic separations, $`z_c`$, between the branes. As an example, for $`kM_{Pl}`$ we only require $`z_c50l_{Pl}`$ to have $`r_c10^{28}`$ cm, the present horizon size of the Universe, i.e., with mild assumptions about $`z_c`$, the four-dimensional description of gravity is valid from the Planck to cosmological scales (in this example, long distance corrections to Newton’s gravity law dominate over short distance ones at $`r\stackrel{>}{_{}}\zeta _c10^{13}`$ cm).

An interesting consequence of the incomplete localization of graviton is the dissipation into the fifth dimension of gravitational waves propagating to large distances. Let us consider gravitational waves generated by a periodic pointlike source on the brane, $`T(x,z)=T(𝐱)e^{i\omega t}\delta (z)`$, where the four-dimensional indices are again omitted. The gravitational field on the brane is given by the convolution of the source, $`T(𝐱)`$, with the Green’s function

$$G(𝐱𝐱^{};\omega )=8\pi G_5d(t^{}t)G(x,x^{};z=z^{}=0)e^{i\omega (tt^{})}$$
(24)
Here the five-dimensional gravitational constant is included for convenience of comparison with four-dimensional formulae, and $`G(x,x^{};z,z^{})`$ is the retarded Green’s function of the linearized Einstein equations. The latter is constructed from the full set of eigenmodes in the usual way:

$$G(x,x^{};z,z^{})=_0^{\mathrm{}}𝑑m\psi _m(z)\psi _m(z^{})\frac{1}{(2\pi )^4}d^4p\frac{e^{ipx}}{m^2p^2iϵp^0}$$
(25)
After substitution of (25) into (24) and simplifications we get

$$G(𝐱𝐱^{};\omega )=\frac{2G_5}{r}_0^{\mathrm{}}𝑑m\psi _m^2(0)e^{ip_\omega r}$$
(26)
where $`r=|𝐱𝐱^{}|`$, $`p_\omega =\sqrt{\omega ^2m^2}`$ when $`m<\omega `$ and $`p_\omega =i\sqrt{m^2\omega ^2}`$ when $`m>\omega `$. We see that the gravitational field on the brane has the form of a superposition of massive four-dimensional modes. Only modes with $`m<\omega `$ are actually radiated, the other ones exponentially fall off from the source. Thus, as long as we are interested in gravitational waves, we can integrate in (26) only up to $`m=\omega `$.

Let us study the case $`r_c^1\omega k`$. Then, the main contribution to (26) is given by modes with $`mr_c^1`$, whereas modes with larger masses give rise to corrections suppressed by $`\omega /k`$ (cf.(16), (17)). In this region of $`m`$, the eigenfunctions $`\psi _m`$ are given by the explicit expressions (7), (14) and $`p_\omega `$ is approximated by $`\omega \frac{m^2}{2\omega }`$. In this way we get (cf.(20))

$$G(𝐱𝐱^{};\omega )=\frac{2G_N}{r}e^{i\omega r}\frac{2}{\pi }_0^{\mathrm{}}𝑑x\frac{e^{i\frac{2r}{\omega r_c^2}x^2}}{1+x^2}$$
(27)
where again $`x=mr_c/2`$ and we extended the integration to infinity. This integral is expressed in terms of the error function:

$$G(𝐱𝐱^{};\omega )=\frac{2G_N}{r}e^{i\omega r}\left(1\mathrm{erf}(\beta )\right)e^{\beta ^2}$$
(28)
where $`\beta =e^{i\frac{\pi }{4}}\sqrt{\frac{2r}{\omega r_c^2}}`$. When $`\beta 1`$ (that is $`r\omega r_c^2`$) one has $`\mathrm{erf}(\beta )0`$, and we obtain the usual $`1/r`$-dependence of the gravity wave amplitude on the distance to the source. In the opposite case $`\beta 1`$ ($`r\omega r_c^2`$), we make use of $`\mathrm{erf}(\beta )=1\frac{1}{\sqrt{\pi }\beta }e^{\beta ^2}`$ to obtain

$$G(𝐱𝐱^{};\omega )=\frac{2G_N}{r^{3/2}}\sqrt{\frac{\omega r_c^2}{2\pi }}e^{i\omega ri\frac{\pi }{4}}$$
(29)
that is, the amplitude is proportional to $`r^{3/2}`$. Thus, the gravitational waves dissipate into the fifth dimension and from the point of view of a four-dimensional observer on the brane, energy of gravity waves is not conserved.

This effect becomes considerable after the wave travels the distance of order $`r\omega r_c^2`$ from the source. Note that at $`\omega r_c^1`$ this distance is much larger than the distance $`r_c`$ at which the violation of four-dimensional Newton’s law is appreciable. This difference between the two distance scales can in fact be seen to be a relativistic effect. The collection of five-dimensional graviton states with $`mr_c^1`$ may be viewed as an RS bound state which becomes metastable in our model. The width of this metastable state is of order $`\mathrm{\Gamma }=\mathrm{\Delta }mmr_c^1`$. In its own reference frame, the graviton disappears into the fifth dimension with time scale $`\tau \mathrm{\Gamma }^1r_c`$. This time scale determines the distance at which four-dimensional gravity of static sources gets modified. On the other hand, when the graviton moves in four dimensional space-time with momentum $`p\omega `$, there is an additional gamma-factor of order $`\gamma \omega /m\omega r_c`$. The graviton therefore remains effectively four-dimensional during a time of order $`\gamma \tau \omega r_c^2`$ due to the relativistic time delay.

Because of this property, the leaking of the gravity waves into the extra dimension is negligible at relatively short wavelengths. However, the corresponding time scale becomes of order $`r_c`$ for wavelengths of the same order. This is another manifestation of our observation that extra dimensions open up at the length scale $`r_c`$.

Finally, we point out that in more complicated higher-dimensional models, the long distance properties of gravitational interactions may be even more intriguing. Indeed, let us modify the model discussed throughout this paper by compactifying the fifth dimension to a very large radius, $`z_{}z_c`$. Then the mass spectrum of Kaluza–Klein gravitons becomes discrete, and the graviton zero mode reappears. However, the spacing between the masses may be tiny, depending on $`z_{}`$. With appropriately chosen $`z_{}`$, the five-dimensional gravity law, Eq. (22), will itself be valid in a finite interval of distances, and the four-dimensional gravity will again be restored at largest scales, well above $`r_c`$.

We conclude that higher-dimensional theories provide, somewhat unexpectedly, valid alternatives to four-dimensional gravity at large distances. With hindsight, it is perhaps not surprising that this happens, since at the very large scale, the anti-de Sitter sandwich becomes very slim, and spacetime is nearly flat, however, from the perspective of our putative universe – the four-dimensional central wall – this conclusion is not so transparent. It would be worthwhile to explore such models in the context of cosmology and astrophysics. The long-distance phenomena in our Universe may become a window to microscopic extra dimensions!

We would like to thank Victor Berezin, Sergei Dubovsky, Dmitry Gorbunov, Maxim Libanov and Sergei Troitsky for useful discussions. We would also like to thank Ian Kogan for drawing our attention to Ref. . R.G. and V.R. acknowledge the hospitality of the Isaac Newton Institute for Mathematical Sciences, where this work has begun. R.G. was supported in part by the Royal Society, and V.R. and S.S by the Russian Foundation for Basic Research, grant 990218410.