text:
# 1 Introduction

## 1 Introduction

At ambient pressure cesium metal crystallises in a bcc structure with experimental lattice parameter $`a_0=6.048\AA `$ at 0K, obtained from the high-temperature values by extrapolation. It is a simple $`s^1`$ metal with nearly free-electron character. Under pressure the $`s`$ valence electrons are transferred to more localised $`d`$-like states, and cesium exhibits an interesting sequence of phase transitions\[3-8\]. At the pressure 2.3GPa a transition from bcc (Cs-I) to fcc (Cs-II) occurs with a small reduction of volume. Near 4.2GPa, the fcc Cs-II phase undergoes an isostructural transition to a collapsed fcc phase Cs-III with a large volume reduction (9%). Then at 4.4GPa there follows a transition to the tetragonal Cs-IV phase with a decrease of the coordination number from 12 in Cs-III to 8 in Cs-IV (strictly, 4 nearest neighbours at 3.349Åand 4 second-nearest neighbours at 3.542Åat 8 GPa). At higher pressures more complicated phases appear.

Many detailed investigations have been devoted to the study of the phase transitions in cesium, both experimentally (see Ref. and literature cited therein) and theoretically (see, e.g., Refs.\[7-12\]). The pressure-induced electronic $`sd`$ transition was believed to be the driving force for destabilising the highly symmetric low-pressure structures (bcc or fcc) with respect to lower symmetry structures. The transition III$``$IV was attributed to the accelerated progress of the $`sd`$ transition and the unusual decrease of the coordination number from Cs-III to Cs-IV has been interpreted in terms of directional bonding induced by the $`sd`$ transition.

In general, the $`sd`$ electronic transition in cesium affects not only the static lattice properties by changing the electronic band structure but also the lattice vibrational properties by modifying the effective interatomic interactions. The structural behaviour of cesium may therefore be reflected in the phonon spectrum and it is therefore of interest to investigate the phononic behaviour with applied pressure. Glötzel and McMahan suggested that an anomaly in the lattice vibrational contribution to the pressure is a possible mechanism for the isostructural II$``$III transition. Recently, Christensen et al. calculated phonon dispersions of fcc-Cs within the quasiharmonic approximation to study the thermal expansion coefficient of cesium and the isostructural transition. They found that below $`V/V_0=0.375`$ ($`V_0`$ is the experimental volume of bcc-Cs at ambient conditions) fcc-Cs becomes unstable due to softening of a transverse  phonon branch, which is characterised by a negative shear elastic constant $`C^{}`$. Very recently, Xie et al. investigated the phonon instabilities of cesium using density-functional perturbation theory. Their calculations showed that the instability of the fcc cesium occurs at a volume of about 0.46$`V_0`$, which is also driven by a negative $`C^{}`$. However, it is not clear what happened with the phonons before the $`C^{}`$ goes soft. Using the density functional linear-response LMTO method, we calculate the accurate phonon spectrum of fcc-Cs at the volumes near the III-IV transition to study the lattice vibrations in fcc-Cs and to seek a possible dynamical instability of fcc-Cs before $`C^{}`$ becomes negative. According to our calculations, fcc-Cs becomes unstable at a volume slightly smaller than $`V/V_0=0.41`$. Soft $`T_{\left[1\overline{1}0\right]}`$-$`\left[\xi \xi 0\right]`$ phonons with $`\xi \frac{1}{3}`$ are found to be responsible for the dynamical instability.

This paper is organised as follows. Some computational details concerning linear-response LMTO calculations are described in the next section. Section 3 presents the results obtained and some discussions. Firstly, the phonon dispersion curves for fcc-Cs at the volume $`V/V_0=0.44`$ are analysed in terms of force constants between atoms; then the phonon frequencies at various volumes are discussed in relation to the instability of the fcc structure. Several concluding remarks are given in the last section.

## 2 Computational Details

The density functional linear-response LMTO method is used in the present study to calculate the phonon spectrum of fcc-Cs at the volumes around the Cs-III - Cs-IV transition. Phonon dispersion curves for a large number of simple and transition metals and a few compounds have been calculated by the linear-response LMTO method, and excellent agreement with experimental data were obtained. In the following we summarise some details of the calculations.

The dynamical matrix of fcc-Cs is calculated for a set of irreducible q points in a (8,8,8)-reciprocal lattice grid (29 points per 1/48th part of the Brillouin Zone (BZ)). The (I,J,K) reciprocal-lattice grid is defined in the usual manner: $`𝐪_{ijk}=\left(i/I\right)𝐆_1+\left(j/J\right)𝐆_2+\left(k/K\right)𝐆_3`$, where $`𝐆_i`$ are the primitive translations in reciprocal space. We chose the exchange-correlation potential of Vosko-Wilk-Nusair plus the non-local generalised-gradient-approximation (GGA-96) correction. Because this gave a better prediction of the equilibrium volume at ambient pressure and 0K which is only 4% smaller than the experimental zero-pressure volume at 0K. On the other hand, the local density approximation (LDA) gave a nearly 20% overbinding at ambient pressure. The calculated bcc to fcc transition pressure of 2.0GPa is in good agreement with the experimental value of 2.3GPa.

In the calculations a $`3\kappa spd`$ LMTO basis set is used with the one-center expansions inside the non-overlapping muffin-tin (MT) spheres performed up to $`l_{max}=6`$. In the interstitial region, the $`s`$-, $`p`$\- and $`d`$-basis functions are expanded in plane waves. The $`5s`$ semicore state is treated as valence state in a separate energy window. The induced charge densities and screened potentials are represented inside the MT spheres by spherical harmonics up to $`l_{max}=6`$ and by plane waves with a 148.4 Ry energy cutoff (9984 plane waves) in the interstitial region.

The $`𝐤`$-space integration is performed over a (16,16,16) grid (145 irreducible point) by means of the improved tetrahedron method, but the integration weights for these $`𝐤`$ points are calculated from a denser (32,32,32) grid (897 points in the irreducible BZ). This results in a more accurate representation of the Fermi surface with a small number of $`𝐤`$ points.

The phonons along the high-symmetry lines presented in the next section were calculated in a denser $`𝐪`$ mesh which fit to the (16,16,16) $`𝐤`$-grid.

## 3 Results and Discussions

### 3.1 Phonon dispersions at $`V/V_0=0.44`$

Fig. 1 shows the calculated phonon dispersion curves along some high-symmetry directions for fcc-Cs at the volume 0.44$`V_0`$, which is in the experimental volume range of the cubic to tetragonal phase transition. To the right the calculated phonon density of states (DOS) is plotted. The dash lines represent the fitted results by the Born-von Kármán model which is described in the next section. The calculated phonon frequencies at the high-symmetry zone-boundary (ZB) points $`L,X,K`$ and $`W`$ are listed in Table 1.

As observed in Fig. 1, the phonon frequencies of the $`L\left[\xi \xi 0\right]`$, $`T\left[00\xi \right]`$ and $`T\left[\xi \xi \xi \right]`$ branches exhibit rather normal dispersion while the $`T\left[\xi \xi 0\right]`$, $`L\left[00\xi \right]`$ and $`L\left[\xi \xi \xi \right]`$ dispersion curves show abnormal behaviour for intermediate and short wavelengths, which indicates long-range interactions between atoms. For the $`L\left[\xi \xi \xi \right]`$ phonon branch, small downwards anomaly is seen around $`\xi =\frac{1}{4}`$. Near the zone boundary, the $`L\left[\xi \xi \xi \right]`$ phonons have very little dispersion and the frequencies of the $`L\left[00\xi \right]`$ phonons are rather small compared to other fcc metals. The $`T\left[\xi \xi 0\right]`$ phonons have flat regions around $`\xi =\frac{1}{3}`$ and $`\frac{3}{4}`$ for the polarisation $`\left[1\overline{1}0\right]`$ ($`T_1`$) and $`\left[001\right]`$ ($`T_2`$), respectively. Since the phonon DOS in Fig. 1 was calculated at a coarse (8,8,8) grid, these abnormal behaviours are not clearly seen in the plotted DOS. Nevertheless, the sudden increase of DOS at the frequency just above 0.3 THz and the small sharp peak at the frequency near 1.6 THz indicate the presence of these anomalies.

In metals at equilibrium the long-range interatomic interactions are generally of little importance due to screening effect of the electrons. The lattice dynamical properties of metals can therefore be described well using a short-range force model. This, however, may not hold for cesium under high pressure. The pressure-induced electronic $`sd`$ transition causes the valence charge density of cesium to deviate from spherical symmetry with a consequent reduction in electronic screening and increased interactions between the atoms. To check this, we calculate the phonon dispersion curves of fcc-Cs at the volume $`V/V_0=0.44`$ using a short-range force model within the quasiharmonic approximation. The dispersion relations are derived as in the usual harmonic approximation by considering non-central forces between atoms up to second neighbours and the four interatomic force constants involved in the dynamical matrices are obtained from the three calculated elastic constants and the $`T\left[100\right]`$ zone boundary phonon frequency by expanding the elements of the dynamical matrix and finding its long-wavelength limit. The obtained dispersion curves are plotted in Fig. 1 by solid lines and the phonon frequencies at the zone boundaries are listed in Table 1. Note that the quasiharmonic results presented here are derived from a short-range force model, which does not exactly consider the electronic contribution to the dynamical matrix. The difference to the results from linear-response calculations may therefore reveal the effects of long-range interatomic interactions. It may be seen in Fig. 1 that the short-range non-central force model describes the $`T\left[00\xi \right]`$, $`T\left[\xi \xi \xi \right]`$ and $`L\left[\xi \xi 0\right]`$ phonon dispersions quite well. However, the previously mentioned anomalies in the $`L\left[00\xi \right]`$, $`L\left[\xi \xi \xi \right]`$ and $`T\left[\xi \xi 0\right]`$ dispersion curves calculated by the linear-response method are not reproduced and they may therefore be a consequence of long-range interactions.

In general, an abnormal behaviour of the phonon dispersion may be attributed to special properties of interatomic force constants. In the following subsection, we shall analyse the dispersion curves by using the Born-von Kármán model to find these special properties of the interactions between the atoms in fcc-Cs.

### 3.2 Born-von Kármán analysis of force constants

Since one of the motivations of the present work is to interpret the anomalies in the calculated phonon dispersion curves for fcc-Cs, it is necessary to use a model which can reproduce these anomalies. This could be achieved by a force-constant model based on the Born-von Kármán theory.

In the Born-von Kármán theory the frequencies $`\nu \left(𝐪\right)`$ of the normal vibration modes for an fcc lattice with one atom per unit cell are given as the solutions of a $`3\times 3`$ determinantal equation

$$\left|4\pi ^2M\nu ^2\left(𝐪\right)\delta _{\alpha \beta }C_{\alpha \beta }\left(𝐪\right)\right|=0,$$
(1)
where $`M`$ is the mass of the atom and

$$C_{\alpha \beta }\left(𝐪\right)=\underset{l^{}}{}\mathrm{\Phi }_{\alpha \beta }(l,l^{})exp\left[i𝐪𝐑(l,l^{})\right],$$
(2)
where $`\mathrm{\Phi }_{\alpha \beta }(l,l^{})`$ is the $`\alpha \beta `$ component of the force constant matrix between the atoms in the $`l`$th and $`l^{}`$th cell. For the phonon modes along the three high-symmetry directions $`\left[\xi \xi 0\right],\left[00\xi \right]`$ and $`\left[\xi \xi \xi \right]`$, the $`C_{\alpha \beta }`$ matrices are diagonal and the solution of Eq.(1) leads to equations in $`\nu ^2`$ of the form

$$4\pi ^2M\nu ^2=\underset{n}{}\mathrm{\Phi }_n\left[1cos\left(n\pi q/q_{max}\right)\right],$$
(3)
where $`q_{max}`$ is half the distance to the nearest reciprocal lattice point in the direction of $`𝐪`$, $`q=\left|𝐪\right|`$ and $`\mathrm{\Phi }_n`$ is a sum of $`\mathrm{\Phi }_{\alpha \beta }(l,l^{})`$ for which the phase $`𝐪𝐑(l,l^{})`$ is a constant. $`\mathrm{\Phi }_n`$ effectively represents a force between a plane of atoms and the planes of atoms normal to $`𝐪`$ and $`n`$ planes away. The summation includes $`N`$ terms so that $`\mathrm{\Phi }_n=0`$ for $`n>N`$. Thus, a Fourier series analysis of the squares of the phonon frequencies will yield the interplanar force constants.

Accordingly, a least-squares Fourier fit are performed using Eq.(3) for the phonon dispersion curves of fcc-Cs at the volume 0.44$`V_0`$ calculated from linear-response calculations. The fitted dispersion curves are plotted in Fig. 1 by the dashed lines. Within the accuracy of the calculated phonon frequencies, it is found that a satisfied fit can only be obtained by including at least four planes in the $`\left[00\xi \right]`$ branches, five planes in the $`T\left[\xi \xi \xi \right]`$ branch and even six planes in the $`\left[\xi \xi 0\right]`$ and $`L\left[\xi \xi \xi \right]`$ branches. In Table 2 we list the fitted interplanar force constants $`\mathrm{\Phi }_n`$. Since four planes along the $`\left[00\xi \right]`$ direction correspond to interactions out to at least eighth neighbour atoms and five (six) planes along the $`\left[\xi \xi \xi \right]`$ ($`\left[\xi \xi 0\right]`$) direction out to farther neighbours, the fitted results confirm the presence of long-range interactions between atoms in fcc-Cs. We plot the fitted $`\mathrm{\Phi }_n`$ for the $`T_{\left[1\overline{1}0\right]}`$\[$`\xi \xi `$0\] branch in Fig. 2a as a function of $`n`$, which effectively is the distance between planes of atoms. The $`\mathrm{\Phi }_n`$ are observed to oscillate with the distance between planes, especially with a prominently negative $`\mathrm{\Phi }_{n=2}`$. The values of $`\nu ^2`$ for the same phonon branch and the fitted curve with $`N=8`$ are shown in Fig. 2b together with the individual contributions from the first five Fourier components. It is seen that the $`n=2`$ planes produce a largely negative contribution to the frequencies of phonons, which partly cancels the dominant contribution from the $`n=1`$ planes. Consequently, the soft intermediate-wavelength $`T_{\left[1\overline{1}0\right]}\left[\xi \xi 0\right]`$ phonons result and a flat region of the dispersion curve is formed. A similar oscillating behaviour of the fitted $`\mathrm{\Phi }_n`$ is seen in some of the other branches although they are less pronounced. The anomalies observed in the $`T_{\left[001\right]}\left[\xi \xi 0\right]`$ and $`L\left[\xi \xi \xi \right]`$ dispersion curves can be associated with these oscillations. For the $`L\left[00\xi \right]`$ branch no oscillations in $`\mathrm{\Phi }_n`$ is found. The very low phonon frequencies for long-wavelengths compared to the short-range model calculation may be attributed to the small value of $`\mathrm{\Phi }_{n=1}`$.

The interplanar force constants can be used to derive interatomic force constants. With a general force model, a linear least-square fitting analysis for the $`\mathrm{\Phi }_n`$ in Table 2 only allows us to work out the interatomic force constants in fcc-Cs extending to the fourth neighbour. The truncation of the interactions with farther neighbours introduces large errors in the fitting procedure and makes the fitted interatomic force constants less reliable. Nevertheless, in Table 3 we list the best fit values of the interatomic constants in fcc cesium.

### 3.3 Soft mode and lattice instability

The volume dependence of the high-symmetry zone-boundary $`L,X,K`$, and $`W`$ phonons for fcc-Cs are presented in Fig. 3. From this the mode-Grüneisen parameter $`\gamma _i=\left(ln\nu _i/lnV\right)`$ was derived and listed in Table 1. It may be seen that all the high-symmetry phonons except the transverse $`X`$ and $`L`$ phonons decrease with increasing pressure and consequently have negative mode-Grüneisen parameters, which may be an indication that the fcc structure is unstable under pressure.

In Fig. 4 we show the $`T_{\left[1\overline{1}0\right]}`$\[$`\xi \xi `$0\] phonon-dispersion curves of fcc-Cs for the volumes $`V/V_0=0.37,0.40,0.41`$, and $`0.44`$. At volumes between $`V/V_0=0.41`$ and 0.40 the phonon frequencies around $`\xi =\frac{1}{3}`$ become imaginary. However, the shear elastic constant $`C^{}`$ does not become negative before the volume is reduced to $`V/V_0=0.37`$ where the entire phonon branch becomes unstable. Christensen et al. and Xie et al. inferred from their calculations that it is the pressure-induced negativity of $`C^{}`$ which is the cause of the instability of the fcc phase. They however did not trace the development of the instability in details and it is clear from Fig. 4 that it is a $`𝐪0`$ $`T_{\left[1\overline{1}0\right]}`$\[$`\xi \xi `$0\] soft mode which drives the dynamical instability of the fcc phase close to $`V/V_0=0.41`$. As a consequence of the soft mode, superstructures should exist as a precursor to the tetragonal Cs-IV phase transition. It is interesting to note that the calculated transition volume $`V/V_0=0.410.40`$ is close to the volume where the fcc to tetragonal transition is observed. In Fig. 5 we illustrate the displacement pattern of atoms in fcc-Cs corresponding a $`T_{\left[1\overline{1}0\right]}`$\[$`\frac{1}{4}\frac{1}{4}`$0\] phonon mode and the resulting Cs superstructure. As seen from Fig. 5, a structural element of the displacement pattern is triangular prisms, which are also present in the Cs-IV structure, even though other possible soft modes as well as the coupling between phonons could be important.

Usually a soft phonon mode leads to a second-order or nearly first-order phase transition. Early experimental results indicated that the Cs-III-Cs-IV transition is first-order in character with a 4.3% volume reduction. Therefore soft phonon modes should not be considered the unique driving force in the cesium III-IV phase transformation even though the soft mode leads to an instability in the fcc structure. The changes of the electronic properties under pressure, such as $`sd`$ transition, must play an important role.

## 4 Conclusions

Using a density functional linear-response LMTO method, we have studied the phonon spectra of fcc cesium at volumes near the III-IV transition. Several anomalies in the dispersion curves have been found.

Within the Born-von Kármán theory of lattice dynamics, the observed anomalies in the calculated phonon dispersion curves for fcc-Cs at the volume 0.44$`V_0`$ are connected with oscillating behaviours of the interplanar force constants against plane spacing $`n`$. In particular, a negative (110)-interplanar force-constant $`\mathrm{\Phi }_{n=2}`$ is found to be responsible for soft intermediate-wavelength $`T_{\left[1\overline{1}0\right]}\left[\xi \xi 0\right]`$ phonons.

The frequencies of high-symmetry $`K`$ and $`W`$ and longitudinal $`X`$ and $`L`$ phonons decrease with volume compression and the consequent negative mode-Grüneisen parameters indicate an instability of fcc lattice under pressure. When the volume is below $`0.41V_0`$, the frequencies of the $`T_{\left[1\overline{1}0\right]}\left[\xi \xi 0\right]`$ phonons with $`\xi `$ around $`\frac{1}{3}`$ becomes imaginary prior to the sign change of the shear elastic constant $`C^{}`$ and the fcc lattice is thus dynamically unstable. It is suggested that superstructures corresponding to the $`𝐪0`$ $`T_{\left[1\overline{1}0\right]}\left[\xi \xi 0\right]`$ soft phonon should be present as a precursor of tetragonal Cs-IV structure. It was pointed out in Ref. that the fcc structure of Cs-III could not be confirmed and that Cs-III appears to be a rather complicated structure. The present argument may stimulate further experimental study for reexamining the structure of Cs at the pressure near the Cs-II-Cs-III-Cs-IV transitions.

## Acknowledgements

Y. Kong gratefully acknowledge helpful discussions and advice from Prof. O. K. Andersen, and we are grateful to Prof. K. Syassen for a critical reading of the manuscript.

## Tables

## Figure Captions

Figure 1: Calculated phonon-dispersion curves along three symmetry directions for fcc-Cs at the volume $`V/V_0=0.44`$. The squares are the calculated results using the density-functional linear-response method and the solid lines represent the results from a quasiharmonic calculation (see text). The dash lines are the fitted results by the Born-von Kármán model. Also plotted in the right panel is the phonon density of states derived from the linear-response calculations.

Figure 2: (a) The interplanar force constants (up to $`n=8`$) for the $`T_{\left[1\overline{1}0\right]}`$\[$`\xi \xi `$0\] branch phonons of fcc-Cs at the volume $`V/V_0=0.44`$. A clear oscillating behaviour against the number of the plane (interplanar spacing) is seen. (b) Values of $`\nu ^2`$ for the same branch phonons plotted as a function of the reduced wave-vector. The fitted curve (solid line) with 8 planes are shown together with the contributions from the first five Fourier components.

Figure 3: Calculated phonon frequencies $`\nu _i`$ (THz) at the high-symmetry L, X, K and W points for fcc Cs as a function of the volume. (a) Longitudinal and (b) transverse branches.

Figure 4: Calculated $`T_{\left[1\overline{1}0\right]}`$\[$`\xi \xi `$0\] phonon-dispersion curve for fcc Cs at the volumes $`V/V_0=0.37,0.40,0.41`$, and $`0.44`$. The points are the calculated values and the lines result from the interpolation between points and are only a guide to the eye.

Figure 5: Schematic illustration of displacement pattern (left) of atoms in fcc-Cs corresponding a $`T_{\left[1\overline{1}0\right]}`$\[$`\frac{1}{4}\frac{1}{4}`$0\] phonon and the resulting Cs superstructure (right). The lengths of solid arrows are proportional to the displacement of atoms and the dash arrows indicate the ,  and  directions of fcc lattice.


format: text

text:
# Dynamic scaling theory of the Kosterlitz-Thouless-Berezinskii transition: ubiquitous finite size effects

## I Introduction

The universality class of the two-dimensional (2D) XY model has been fascinating from the start with its phase transition being driven by topological excitations, the unusual, exponential temperature dependence of the correlation length, and its apparent presence in systems ranging from 2D superfluids and superconductors to liquid crystals and magnetic systems. Despite its uniqueness, the detection of the Kosterlitz-Thouless-Berezinskii (KTB) transition, characterized by an unbinding of vortex pairs, has been challenging, but several techniques have been developed to probe its critical behavior. In two-dimensional superconductors and Josephson Junction arrays (JJA’s), an examination of the dc transport properties in zero magnetic-field has been the primary method, but the magneto-resistance and the ac kinetic inductance have also been used. More recently, magnetic flux noise measurements have been used to study the KTB transition. And, as we discuss shortly, numerical methods and scaling approaches have also been used. The advantage of the scaling methods is that they can reveal information on the dynamic critical exponent $`z`$, one of the key exponents in the KTB critical behavior, which relates the relaxation time $`\tau `$ to the vortex correlation length $`\xi `$ via $`\tau \xi ^z`$.

While there are many methods for studying the KTB transition, the most common is the “conventional” method which involves an analysis of the dc voltage-current ($`V`$-$`I`$) isotherms and the temperature dependent resistance $`R(T)`$. In this approach, the transition temperature $`T_c`$ is determined via the $`V`$-$`I`$ exponent $`\alpha (T,I)`$ ($`VI^\alpha `$), which is predicted to have a value of 3 at $`T_c`$ (assuming $`z=2`$,) and whose value should decrease rapidly for weak currents near $`T_c`$. The Ohmic resistance is then fitted by the Kosterlitz or Minnhagen form for the resistance, which are respectively, $`R\mathrm{exp}[bz/\sqrt{TT_c}]`$ and $`R\mathrm{exp}[bz\sqrt{(T_{c0}T_c)/(TT_c)}]`$, where $`b`$ is a system-specific parameter and $`T_{c0}`$ is the mean-field transition temperature.

There are many weaknesses to the conventional approach. For example, it is valid only in the weak-current limit where few experiments can probe effectively and where finite-size effects dominate. Moreover, the exponent $`\alpha `$ has a current dependence and so determining its value uniquely for a given temperature is not possible. In fact, researchers can use this flexibility to their advantage to show a rapid variation of $`\alpha `$ with temperature near $`T_c`$. This is done by determining $`\alpha `$ in different current ranges for each temperature so as to achieve the desired behavior. (The most accurate method would be to determine $`\alpha `$ for the $`V`$-$`I`$ isotherms for the same current so as to represent a single length scale.) Another problem with the traditional approach is the number of fitting parameters in the resistance formulas, which allow excessive freedom for interpreting results. For these reasons, it is important to examine other approaches to verifying and studying KTB behavior.

The Fisher-Fisher-Huse (FFH) dynamic scaling, which says that $`V\xi ^z/I`$ is a function of $`I\xi /T`$, has its advantages over the conventional approach. The most prominent is that it is valid over an extended current range and therefore incorporates the current dependences of $`\alpha (T,I)`$. Moreover, FFH scaling can be used to determine $`z`$. Recently, it has been used to study experimental, dc $`V`$-$`I`$ data but yields a large value of the dynamic critical exponent $`z5.6`$ for a wide variety of samples and systems. It has recently been reiterated that the FFH scaling is only valid in the absence of finite size effects and suggested that finite size could result in large values of $`z`$.

Finite size effects, which include a finite size system and the inherent magnetic penetration depth $`\lambda _{}`$ in superconductors, have long been known to complicate the verification of KTB behavior because they make free vortices possible for any temperature range. Kadin et al. and Lee and Garland were among the first to study this experimentally while Repaci et al. and Herbert et al. have recently extended this work. Theoretically, Simkin and Kosterlitz have investigated numerically and analytically the finite size effects in these systems, while Artemov has examined the effect of “Pearl” vortices on the magnetization in superconducting thin films. And more recently, Pierson and Valls have completed a renormalization group study of the two-dimensional Yukawa Gas, which incorporates the inherent finite size $`\lambda _{}`$ into the two-dimensional Coulomb Gas (2DGG).

Finite size (FS) scaling, which incorporates the effect of a finite size into the FFH scaling, is an established method in the study of numerical data in KTB behavior. Lee and Teitel were among the first to use this approach to show that the dynamic critical exponent $`z`$ has a value $`2`$ for the two-dimensional Lattice Coulomb Gas (2DLCG). Medvedyeva et al. have also used finite size scaling to reinforce the result $`z=2`$ for the 2D resistively-shunted-junction model. (Numerical studies in fact have provided many insights into KTB behavior. A review of those results is beyond the scope of this paper but we refer the reader to many excellent papers in the literature.)

Like the conventional approach and FFH scaling, FS scaling has its strengths and weaknesses. For example, because FS scaling requires the ability to carefully tune the system size, it is not conducive to studying experimental data where the finite length scales are more illusive. It is therefore important to keep in mind that all approaches provide a great deal of information and often complement one another.

The correlation length is an important quantity in FFH scaling but questions about its form below $`T_c`$ remain. In the original paper of Kosterlitz, the correlation length was taken to be infinite below the transition temperature owing to the divergence of the susceptibility there. In subsequent papers by Friesen, Pierson, and Simkin and Kosterlitz, the vortex correlation length is taken to be finite for this temperature regime (but still diverging at $`T_c`$ from below) and interpreted roughly as the size of the largest vortex pairs. Other authors maintain that the correlation length is infinite for $`T<T_c`$.

In this paper, we do a numerical study of the 2DLCG using Monte Carlo simulations to understand the difference in the value of $`z`$ obtained from finite size scaling and FFH scaling. We vary the system size to study the effect of finite sizes and verify that $`z=2`$ using FFH scaling, thereby making FFH scaling consistent with FS scaling. We also analyze the transport characteristics and find evidence that the vortex correlation length is finite below the phase transition, in contrast to the claims of others. The numerical approach that we use is that of Lee and Teitel extended to larger system sizes ($`L=200`$).

The paper is organized as follows: Section II describes the model that we use along with the numerical and scaling methods. In Section III, we present our results and analysis of the transport characteristics. Implications and a summary of our work is presented in Section IV.

## II Model and Methods

### A 2D Lattice Coulomb Gas with an applied field

To study the KTB critical behavior and its dynamic critical exponent, we use the 2DLCG model, whose relationship to the KTB transition and mapping to the 2DXY model are well known. This model is described in Refs.  and  and so will be only briefly described here.

The Hamiltonian for the 2DLCG is given by

$$H=\frac{1}{2}\underset{i,j}{}q_iV(𝐫_{ij})q_j,$$
(1)
where the sum is over all pairs of sites in the lattice, $`q_i`$ is the (integer) charge at site $`i`$, $`𝐫_{ij}=𝐫_i𝐫_j`$, and $`𝐫_i`$ is the location of the $`i`$th charge. $`V`$, which solves the discrete form of Poisson’s equation in two dimensions, is the Coulomb potential (and is therefore logarithmic at large distances). Using Fourier tranforms, $`V`$ in Eq. (1) may be replaced by

$$V(𝐫)=\frac{1}{N}\underset{k}{}V_k(e^{i𝐤𝐫}1)$$
(2)
where, for the square lattice considered here,

$$V_k=\frac{\pi }{2\mathrm{cos}(𝐤\widehat{𝐱})\mathrm{cos}(𝐤\widehat{𝐲})}.$$
(3)
The wave vectors are defined in the usual way:

$$𝐤=\frac{2\pi m_1}{L}\widehat{𝐱}+\frac{2\pi m_2}{L}\widehat{𝐲},$$
(4)
where $`m_1,m_2=0,1,2,\mathrm{},L1`$. In addition, since we wish to apply our results to 2D superconductors in zero-magnetic field, we use the constraint that the total charge in the system must be zero, which deals with the infrared divergence as $`𝐤0`$.

In the presence of an electric field, $`𝐄`$, an electric potential energy term of the form $`_iq_i𝐫_i𝐄`$ is added to the Hamiltonian. The new term results in a net movement of charge along the direction of the field $`𝐄`$ that we measure in terms the charge current density, $`𝐉`$, using the approach of Ref. .

The applied field $`E`$ and the resulting vortex current density $`J`$ in the 2DLCG model are analogous to the current $`I`$ and voltage $`V`$ in a superconducting thin film, but in an inverted way due to the mapping between the two systems. In the superconducting thin film, the current is injected into the sample, which causes vortex movement in a direction perpendicular to the current via the magnus force. This results in resistance and therefore a voltage $`V`$. For the 2DLCG, it is the voltage that is applied and that results in the vortex “charge” current in a direction parallel to the electric field. Following the notation of Lee and Teitel, the following correspondences between the two systems apply: $`EI`$ and $`JV`$. In this paper, we will examine the $`E`$-$`J`$ curves of the 2DLCG, which correspond to $`V`$-$`I`$ characteristics in actual superconducting thin films. Unless otherwise noted, the notation that we use will be that of the 2DLCG. Along these lines, the $`E`$-$`J`$ curves will be plotted as $`E`$ versus $`J`$ in our figures, so that they mirror the analogous curves in the 2D superconcuctors that are plotted as $`V`$ versus $`I`$. We will therefore use the notation $`V`$-$`I`$ for the current-voltage curves instead of the more standard $`I`$-$`V`$.

### B Monte Carlo Simulations

The Metropolis  Monte Carlo technique is used to calculate the $`J`$-$`E`$ data based on the Hamiltonian discussion in Section IIA. Briefly, the technique, as it applies to the 2DLCG model with a square ($`N=L\times L`$) lattice and integer charges, is as follows: (See Ref.  for more details.)

1. Create a trial system from the original system by choosing a site, $`i`$, and one of its four neighbors, $`j`$, at random. Add one unit of charge to site $`i`$ and subtract one unit from site $`j`$ so that the net charge in the system is always zero.

2. Compute $`\mathrm{\Delta }H=H_{trial}H_{original}`$. Since only two sites have been altered, most terms in $`H_{trial}`$ and $`H_{original}`$ are the same, simplifying the calculation of $`\mathrm{\Delta }H`$.

3. If $`\mathrm{\Delta }H<0`$ or if $`e^{\mathrm{\Delta }H/T}`$ is greater than a random number uniformly distributed on the interval $`[0,1)`$, accept the change and replace the the original system with the trial system for subsequent trials. Otherwise, leave the system unchanged.

When evaluating the Hamiltonian for a trial system, the minimum image convention and toroidal (periodic) boundary conditions are used. In addition, a “stability” term of the form $`_i(q_i^4q_i^2)`$ is added to the Hamiltonian to minimize the probability of creating multiply charged sites in the system.

For many combinations of $`E`$ and $`T`$, particularly when the charge density is small, a low rate of acceptance of trials is expected. When this occurs, it is convenient to calculate $`\mathrm{\Delta }H`$ without summing over the array. The method employed here, used and described by Lee and Teitel, is to store information on each site in the system in the form of the potential, which is updated for each site after each accepted trial. This technique can improve run times on the order of the reciprocal of the acceptance ratio, often factors of $`O(10^3)`$ or more.

A sequence of $`N`$ trials is called one Monte Carlo Step (MCS). All of our computer runs consist of an initial $`10,000`$ MCS to equilibrate the system during which no data are saved, followed by at least $`500,000`$ MCS, for which the average current is computed. The only exceptions were some large systems ($`L=200`$) at $`E>0.6`$ where a smaller number of MCS was used due to very long run times. At high electric field, the current is large and the shorter runs provide sufficient accuracy. At each $`E`$ and $`T`$, at least five runs were performed. The error bars on the points in Figure 1, Figure 4, and Figure 5 are the standard deviations of the average current from these runs. In some cases at low $`E`$, many additional runs were performed in an attempt to reduce the error bars.

### C Fisher-Fisher-Huse and Finite-Size Scaling Methods

Given the difficulties of the conventional approach to analyzing $`V`$-$`I`$ characteristics discussed in Section I, it is important to use other methods including the powerful scaling methods. In this paper, we will use Fisher-Fisher-Huse scaling and finite size scaling.

In the FFH theory, for 2D superconductors, the $`E`$-$`J`$ (in the notation of the 2DLCG) curves should scale as

$$J=E\xi ^z\chi _\pm (E\xi /T),$$
(5)
where $`\chi _{+()}(x)`$ is the scaling function for temperatures above (below) $`T_{KTB}`$. The two important asymptotic behaviors of $`\chi (x)`$ are $`lim_{x0}\chi _+(x)=\text{const.}`$ (Ohmic limit), and $`lim_x\mathrm{}\chi _\pm (x)x^z`$ (critical isotherm). As established in previous papers, it is advantageous to rewrite Eq. (5) as

$$\frac{E}{T}\left(\frac{E}{J}\right)^{1/z}=\epsilon _\pm (E\xi /T)$$
(6)
(where $`\epsilon _\pm (x)x/\chi _\pm ^{1/z}(x)`$) because only the $`x`$-scale is stretched in Eq. (6). In Eq. (5), both the $`x`$-scale and $`y`$-scale are stretched making it harder to judge a collapse of the scaled data.

The FFH scaling has three fitting parameters $`z`$, $`T_c`$, and $`b`$, which can be determined in the three scaling steps. The first, which yields estimates of $`z`$ and $`T_c`$, is to determine the critical isotherm. In the absence of finite size effects, the critical isotherm is straight on a log-log plot since $`JE^{z+1}`$ for all $`E`$ at $`T=T_c`$. The second step is to estimate $`b`$ in the exponent of the correlation length, which, for experimental data, one can do by examining the resistance since $`R(T)\xi ^z.`$ For numerical data, one assumes $`\xi \mathrm{exp}[C/|TT_c|^{1/2}]`$, which is generally expected to hold for $`T>T_c`$. The final step is to achieve a scaling collapse for Eq. (6), adjusting the three parameters within the error bars established in steps 1 and 2. One should also bear in mind that all these steps can be repeated to optimize the scaling collapse.

In the finite size scaling form, valid for $`\xi >L`$, the system size $`L`$ is the dominant length scale and therefore substitutes for $`\xi (T)`$ in Eq. (5) to arrive at

$$JL^z/E=\chi _\pm (EL).$$
(7)
Strictly, this equation works only for $`T=T_c`$, but is also valid for $`TT_c`$ since leading order corrections are small: $`O(L/\xi )`$. In this powerful approach, there is only one fitting parameter, $`z`$, making its results more definitive.

To achieve the best scaling collapse for the FFH and FS scaling, the principle method that we use is to sort the data according to the value of the scaling variable $`x`$ (whether it be $`E\xi /T`$ or $`EL`$) and then to minimize the difference of the logarithm of the values of two adjacent, scaled points using a least squares method. As a second means of optimization, we tried approximating the $`T>T_c`$ curve to a function and then minimizing the difference between the scaled points and the curve. However, because of the large number of fitting parameters, definitive results could not be achieved with this approach.

For various reasons, not all of the data are used in the optimization. At large values of the applied electric field ($`E0.8`$), the isotherms become ohmic due to a saturation of the vortex density. (See discussion in Section IIIA.) Because this behavior is not due to critical behavior, the data should not scale (collapse) with the other data and so is not included in our optimization. However, we do plot those points in the scaling plots (Figures 2 and 6), even though those points do not collapse. (In Refs.  and , the analogous, large $`I`$ data was not plotted with the scaled data.) In addition, because the dynamic critical exponent $`z`$ has a temperature dependence below $`T_c`$ due to the line of fixed points in the KTB critical behavior, the $`T<T_c`$ data are not used in the optimization of the FFH scaling collapse. This is because the temperature dependence of $`z`$ is not the same as the scaling variable $`x=E\xi /T`$ and so the FFH scaling is not expected to be as effective in this regime. Nonetheless, we do scale the $`T<T_c`$ data assuming a symmetric correlation length. Finally, data are also excluded from the optimization if the error bars are greater than $`50\%`$.

## III Results

The principal objectives of our work are to study the value of the dynamic critical exponent $`z`$ for the 2DCG in the context of Fisher-Fisher-Huse scaling and the influence of finite size effects, and to study the behavior of the correlation length $`\xi _{}(T)`$ below the phase transition.

### A Dynamic Critical Exponent $`z`$

The two principal system sizes used in our simulations were $`L=24`$ and $`L=100`$. The $`J`$-$`E`$ isotherms for the $`L=24`$ system are shown in Figure 1. As one can see the isothermal $`J`$-$`E`$ curves mimic closely the $`V`$-$`I`$ of thin-film superconductors. Reflecting the absence of vortex pairs, the higher temperature curves are primarily Ohmic (i.e., have slope 1 on a log-log plot). The lowest temperature curves are primarily non-Ohmic, indicating the presence of vortex pairs. The intermediate-temperature curves, on the other hand, are Ohmic for small $`E`$, become non-Ohmic at larger values of $`E`$, and abruptly turn Ohmic for the largest values of $`E`$. This of course reflects the fact that with each different value of $`E`$, one is probing a different length scale, which varies as the inverse of $`E`$. (See discussion in last paragraph of Section IIIB.) Consequently, at small $`E`$, long length scales are being probed where free vortices exist and Ohmic behavior is expected for $`T>T_c`$ (or when $`\lambda _{}`$ is finite). For large values of $`E`$, smaller length scales are probed where vortices are still bound in pairs and result in non-Ohmic behavior. The other salient feature of this graph is that the error bars become large at the small-$`E`$ end of the isotherms due to small numbers of vortex pairs.

The temperature of the phase transition in the $`L=24`$ data is not obvious in Figure 1. Using the FFH dynamic scaling criteria for the critical isotherm, which says that $`JE^{z+1}`$ for all values of $`E`$ at $`T=T_c`$, results in $`T_c0.14`$ and $`z6`$. The subsequent FFH scaling for the $`T>T_c`$ data yields an excellent collapse for $`z=6.01`$, $`T_c=0.127`$, and $`C=0.7475`$, as seen in Figure 2. (As noted in Section II, the crossover to Ohmic behavior at large $`E`$ data seen Figure 1 is not due to critical behavior and so should not be used to judge the collapse. It is this large-$`E`$ data that does not collapse with the other data in Figure 2.)

The values of the parameters in this collapse, optimized using our method described in Section II, are consistent with the large values of $`z`$ and low values of $`T_c`$ found previously. However, an examination of the inverse dielectric constant $`ϵ^1`$ and the vortex density reveals no critical behavior at that temperature with or without a current. An inspection of the lattice configurations also reveals no clues of a phase transition at $`T0.13`$. In Figure 3, characteristic vortex configurations in the lattice are show for various temperatures for a relatively large value of the electric field, $`E=0.3`$. Even at the lowest temperature shown here ($`T=0.17`$), there are few vortices and they are clearly bound to one other. (One can also see in that figure that the vortices seem to become unbound above the expected $`T_c`$ of 0.218.) Hence, the apparent FFH scaling collapse at the large value of $`z`$ and the small value of $`T_c`$ does not seem to correspond with a real phase transition. We now investigate if the FFH scaling collapse for these values of $`T_c`$ and $`z`$ could be due to finite size effects, as recently suggested by Medvedyeva et al.

To further explore the value of $`z`$ and the transition temperature, we employed finite-size scaling methods for $`T=0.210`$, a temperature just below the expected value $`T_c=0.218`$. The $`J`$-$`E`$ curves for values of $`L`$ ranging from $`12`$ to $`200`$ are shown in Figure 4. As one can see, the curves become straighter (on a log-log scale) as expected near the critical isotherm when the system size is increased, indicating that the low-$`E`$ Ohmic behavior for the $`L=24`$ data is attributable to finite size effects and that the true critical isotherm occurs at a larger temperature. The finite size scaling, shown in the inset to Figure 4, reinforces this idea. With $`z`$ being the only fitting parameter, the optimal value is 2.15, reinforcing the original results of Lee and Teitel (and the subsequent work of others) that $`z=2`$ and that $`T_c`$ occurs at a much larger value than obtained from the FFH scaling. The fact that the value of $`z`$ here is slightly larger than 2 can be attributed to the fact that the temperature is just below that of the expected value of the transition temperature and that $`z`$ increases in value for $`T<T_c`$.

The value of $`z`$ obtained in the finite size scaling reinforces the idea that the large value of $`z`$ found for the $`L=24`$ data in the FFH scaling is due to finite size effects. We have therefore calculated $`J`$-$`E`$ curves for $`L=100`$. (See Figure 5.) The $`J`$-$`E`$ curves are similar to those of $`L=24`$ both in terms of qualitative behavior and in terms of the error bars. The primary difference is that the intermediate-temperature curves are straighter than those of $`L=24`$ and the critical isotherm now occurs at a temperature $`T0.210`$ for $`z2`$. (Of course, if we could extend the $`J`$-$`E`$ data to smaller $`J`$, the $`J`$-$`E`$ curves would not be straight for $`T<T_c`$ due to finite size effects.) Using these values as a guide, the best collapse for this data with the FFH scaling, shown in Figure 6, was obtained using $`z=1.90`$, $`T_c=0.212`$, and $`C=0.65`$. As one can see, the collapse is quite good for the $`T>T_c`$ data. (Note that four points for $`T=0.22`$ were excluded from this figure due to the fact that the error bars were greater than $`50\%`$.) Even the $`T<T_c`$ data scale well despite the fact that $`z`$ in this regime has a temperature dependence different than that of the scaling parameter $`x`$. The $`T=0.090`$ and $`T=0.12`$ are the exception to the collapse but is most likely due to the fact that that data is well outside the critical region. The fact that the rest of the $`T<T_c`$ data scales reasonably well may be a coincidence or evidence for a finite correlation length. (We present evidence for a finite correlation length for $`T<T_c`$ in the next subsection.) The value of the $`T_c`$ achieved through our methods is slightly smaller than the expected value, which we attribute to finite size effects. We add that an acceptable collapse is also achieved for $`T_c=0.218`$ and $`z=2`$.

We believe that this work is a definitive verification of $`z=2.0\pm 0.2`$ found using FFH scaling and therefore puts FFH scaling in accord with finite size scaling. We conclude that the larger value of $`z`$ (5.6) found for the $`L=24`$ data and elsewhere is due to finite size effects. (See discussion in Section IV.) This work also vindicates the use of FFH scaling for systems free of finite size effects. Lee and Teitel and Harris et al. have previously showed a $`z=2`$ collapse using FFH scaling for their numerical data ($`L=24`$) and JJA data, respectively. However, we do not believe that that work should not be taken as evidence for $`z=2`$ because the scaling collapse was not convincing, they did not vary $`z`$ to optimize the collapse, and it was later shown that the data collapse better for much larger values of $`z`$.

Finally, we note that we calculated the values of the $`J`$-$`E`$ exponent $`\alpha `$ in order to study the predictions of Minnhagen et al. about its temperature dependence. In accord with the work of others, our work indicates that the prediction of Ref.  for the temperature dependence of $`\alpha `$ (usually denoted $`\alpha _{PM}`$) is in better agreement with the data than the original prediction for the temperature dependence of $`\alpha `$ (sometimes denoted $`\alpha _{AHNS}`$.) This agreement further verifies our $`L=100`$ $`J`$-$`E`$ data.

### B Vortex Correlation Length for $`T<T_c`$

Our method for studying the vortex correlation length is to investigate the value of the $`J`$-$`E`$ exponent $`\alpha `$ as a function of $`E`$ by looking at the slope of our $`J`$-$`E`$ curves. More specifically, we study the derivative $`d\mathrm{log}J/d\mathrm{log}E`$, which is equal to $`\alpha `$ when $`JE^\alpha `$. A similar examination by Medvedyeva et al. (see their Figure 2) for system sizes of $`L=8`$, using the analogous $`V`$-$`I`$ curves in the 2D RSJ model, revealed that the slope started at a low value for small $`I`$, peaked, and then descended to a value near 1 at large values of $`I`$. For each of the isotherms with temperatures below the critical temperature, they found the peak to occur at roughly the same value of $`I`$. Because the scaling variable is either $`I\xi `$ when $`\xi L`$ or $`IL`$ when $`L\xi `$, the fact that the peak in $`d\mathrm{log}V/d\mathrm{log}I`$ does not move is an indication that finite size events dominate this temperature regime (and the scaling parameter is $`EL`$). Above $`T_c`$, the authors found the peak to move quickly to larger values of $`I`$, indicating that $`I\xi `$ is the principal scaling variable since $`\xi `$ is decreasing quickly in this regime.

For our $`L=24`$ data, $`d\mathrm{log}J/d\mathrm{log}E`$ behaves differently below $`T_c`$ than does the analogous quantity of Ref. . As one can see in Figure 7, the locations of the peaks as a function of $`E`$ move to smaller values as one increases the temperature from 0.09 to about 0.14. In the temperature range 0.15-0.22, the peaks develop into plateaus.

The peaks and plateaus indicate that there are two competing length scales. On the large length scale side (small $`E`$) is the temperature-dependent correlation length $`\xi _{}(T)`$ while at small length scales (large $`E`$) is the lattice size, which enters because of a saturation of vortex density. A plateau signifies that $`\xi `$ is much larger than the saturation length scale while a peak indicates that $`\xi `$ is approaching the lattice size. The value for the left hand size of the plateau is determined by $`\xi `$ or $`L`$, whichever is smaller. In Figure 7, the point at which $`\xi `$ becomes larger than $`L`$ occurs at about $`T=0.19`$, above which the plateau begins at $`E0.13`$ on the low $`E`$ side. That the value of $`E`$ for which $`d\mathrm{log}J/d\mathrm{log}E`$ peaks or begins to plateau depends on temperature for a broad temperature range below $`T_c`$ is an indication that the scaling variable for $`T<0.19<T_c`$ is $`E\xi `$ and that, therefore, $`\xi _{}(T)`$ is finite below $`T_c`$.

The lattice spacing affects the $`J`$-$`E`$ curves for $`E0.7`$, as can be seen in Figure 7, and determines a maximum value of $`E`$ at which $`d\mathrm{log}J/d\mathrm{log}E`$ peaks. Therefore, the discrete lattice size can also result in the position of the peak in $`d\mathrm{log}J/d\mathrm{log}E`$ not moving. As a consequence, that the peaks would be located all at the same approximate value of $`E`$, as found in Ref. , should not be taken as evidence strictly for finite size effects.

The values of $`E`$ that correspond to the influence of the lattice spacing and the lattice size can be predicted by simple arguments. (See for example Ref. .) The approximate energy for a vortex pair is $`(\mathrm{log}RER)/ϵ`$, which increases for small $`R`$ and then decreases for large $`R`$, thereby peaking at $`R=1/E`$. (This equation explains why one probes different length scales by changing the value of $`E`$ as discussed in Section IIIA.) Since the lattice size corresponds to a value $`1`$, one expects its effects to manifest themselves for $`E1`$, which is what we observe in Figure 7. Similarly, one expects finite size effects to manifest themselves for $`E1/L0.04`$. Indeed, in Figure 7, all of the isotherms are ohmic for $`E0.5`$, in agreement with our prediction.

## IV Discussion and Summary

Using system sizes up to $`L=200`$ in Monte Carlo simulations of the 2DCG, we conclude $`z=2.0\pm 0.2`$ using Fisher-Fisher-Huse scaling and finite size scaling, reinforcing the results of others that $`z=2`$. (See, e.g., Refs. .) We have also presented evidence that $`\xi _{}`$, the correlation length below $`T_c`$, is finite.

Previous analysis of experimental data on 2D superconductors, JJA’s, and superfluids, and of numerical data on the 2DCG by two of these authors using FFH scaling indicated $`z5.6`$. As stated here and in Ref. , it appears that this large value is due to finite size effects. Because a majority of the $`V`$-$`I`$ data from the literature was examined in Ref. , it appears that nearly all experimental $`V`$-$`I`$ measurements on 2D superconductors and JJA’s are dominated by finite size effects. Why the value $`z6`$ should result from finite size effects remains a mystery, and should be studied further.

Because these results indicate that finite size effects are ubiquitous in 2D superconductors and 2D JJA’s, the $`V`$-$`I`$ and $`R`$ data from the literature that were used to verify KTB behavior in the conventional manner needs to be re-examined to account for the finite size effects. In the presence of strong finite size effects like those in these samples, two aspects of the conventional method fail. First, the small $`I`$ behavior of the $`V`$-$`I`$ isotherms would not reveal any quick changes in $`\alpha `$ near $`T_c`$ because finite size effects dominate at small $`I`$. In other words, as is well known, finite size effects wash out the expected jump in $`\alpha `$. Secondly and more importantly, the resistance formulas are no longer valid in the presence of finite size effects. That previous authors could fit their resistance data to the more commonly used Minnhagen formula is most likely due to the number of fitting parameters. Finally, we stress that the most effective means for studying the temperature dependence of $`\alpha `$ should begin by extracting $`\alpha (T)`$ from the $`V`$-$`I`$ curves in 2D superconductors and JJA’s for the same value of $`I`$.

Clearly, in order for this data to be re-interpreted, a better understanding of finite size effects is needed since previous checks for finite size effects in the data of Ref.  were negative. For example, those authors found that the value of $`\lambda _{}(0)`$ derived from a fit to $`\lambda _{}(T)`$ was four times larger than the expected value. Furthermore, for that same data, where finite size effects are now known to dominate, the resistance was found to fit better to the Kosterlitz resistance formula than a resistance formula derived for finite-size dominated Ohmic resistance. Finally, in Ref. , it was found that the length scale at which the crossover from finite-size-induced Ohmic behavior to non-Ohmic behavior was expected to occur for the data of Ref.  did not correspond with the observed value.

###### Acknowledgements.

Acknowledgement is made by SWP to the donors of The Petroleum Research Fund, administered by the ACS, for support of this research.


text:
# Does Dark Matter Exist?

## 1. Introduction

Discrepancies between the apparent baryon content of galaxies and galaxy clusters and their dynamically inferred masses are generally attributed to dark matter (DM). A still larger mass discrepancy exists out to the horizon scale, since Big Bang nucleosynthesis requires a low baryonic mass fraction in the universe (Tytler et al. 2000), which nevertheless appears to have flat geometry (de Bernardis et al. 2000; Hanany et al. 2000). The missing energy density required by conventional cosmology is now supposed to reside not only in non-baryonic particles, but also partly as vacuum energy.

The currently popular $`\mathrm{\Lambda }`$CDM model for the evolution of structure in the universe is quite successful on large scales (e.g. Bahcall et al. 1999), but its predictions on galaxy scales disagree with observed properties of galaxies. A number of lines of evidence suggest that the problem is more fundamental than simply our inadequate understanding of the process of galaxy formation.

## 2. Mass Distributions in the Inner Parts of Galaxies

Before listing the problems of DM, we briefly review two properties of the mass distribution within galaxies that exacerbate some of the difficulties. Our conclusions are controversial, but we find strong evidence for both points.

First, putative galaxy halos seem to have cores of finite, as opposed to cusped, central densities. Density estimates from rotation curves of low-luminosity and low surface brightness (SB) galaxies are most revealing, since these galaxies have the largest mass discrepancies and therefore halos for which compression by baryonic infall is least important. Van den Bosch et al. (2000) stress that the inner slope in HI data can be underestimated because of “beam smearing,” and Swaters et al. (2000) indeed find that optical data often indicate a steeper rise; but their cases either require, or are consistent with, a finite density core – whatever M/L is ascribed to the luminous component. Other well-resolved rotation curves of low-luminosity galaxies (Rubin et al. 1985; Courteau 1997; Sofue et al. 1999; Blais-Ouellete et al. 2000a,b; Palunas & Williams 2000) generally also reveal a gentle rise, indicating that their halos do indeed have low-density cores.

Second, the mass in the inner parts of bright galaxies is dominated by the luminous disk and bulge. “Maximum disk” models are strongly suggested by the impressive match between the shapes of inner rotation curves and predictions from the luminous matter alone (Kalnajs 1983; Kent 1986; Broeils & Courteau 1997; Palunas & Williams 2000), and other evidence is reviewed by Bosma (1999) and Freeman (this meeting). Furthermore, recent work on barred galaxies (Debattista & Sellwood 1998, 2000; Weiner et al. 2000) has shown that any DM halo makes a negligible contribution to the inner rotation curve even after the formation of the disk and bulge, and that this conclusion for barred galaxies also holds for their unbarred counterparts.

## 3. Problems Created by the Dark Matter Hypothesis

The problems we list here mostly relate to properties of galaxies. The first ten are well established and are ordered such that those which are generic to any form of DM come before those specific to particular DM models. The last two are less well established observationally but are potentially just as thorny. The physics of galaxy formation is undeniably messy, creating scope for counter-arguments that may avert individual difficulties. We do not have space here to evaluate the many such arguments that have been advanced, but no combination of such ideas comes anywhere near to resolving all these problems.

1. The “disk-halo conspiracy” (Bahcall & Casertano 1986) describes the absence of a feature in galaxy rotation curves at which the dominant source of central attraction changes from luminous matter to dark. Many galaxies are now known in which the rotation curve does drop somewhat at the edge of the visible disk (e.g. Casertano & van Gorkom 1991; Verheijen 1997; Bosma 1999), but it is extremely rare for the drop to exceed about 10%. Blumenthal et al. (1986) showed that a featureless rotation curve is expected if DM dominates galaxies right to their centers, but it is much harder to understand why the circular orbital speed from the luminous matter, which dominates the inner region (see §2), should be so similar to that from the DM at larger radii. For any galaxy dominated by stars in its center, initial conditions for the dark and luminous matter must be finely tuned to produce a flat rotation curve.
2. Extreme low-SB galaxies lie on the same Tully-Fisher relation (TFR) derived from high-SB galaxies, with somewhat greater scatter (Zwaan et al. 1995; Sprayberry et al. 1995; McGaugh et al. 2000). Thus we observe similar circular speeds in all galaxies of a given luminosity, no matter how widely the luminous material is spread. This amazing result requires that the overall M/L of the galaxy rises with decreasing SB in just the right way so as to preserve a tight relation between total luminosity and circular speed. Either the true M/L of the stellar population changes with surface brightness, which seems unlikely (de Blok & McGaugh 1997), or the DM fraction rises as the luminous surface density declines. The needed variations would be minor if DM dominated in all galaxies, but since stars dominate the mass in the inner parts of high-SB galaxies (§2), eliminating the SB dependence again requires careful tuning.
3. Sanders (1990), Milgrom (1998) and McGaugh (1999a) show that mass discrepancies begin to be detectable only when the acceleration drops below $`10^8\text{ cm s}^2`$. Any DM model must reproduce this characteristic acceleration scale over a wide range of galaxy sizes, but none has yet done so convincingly.
4. Aside from the disk M/L, DM halo fits to rotation curves generally employ two extra parameters: e.g. the core radius and asymptotic velocity, or the scale radius and concentration index. Actual galaxy rotation curves do not require all this freedom, however, since they can be fitted with only the disk M/L as a parameter (e.g. Sanders & Verheijen 1998) when a modified gravitational force law of the MOND form (Milgrom 1983) is employed. Any DM model must therefore contain a physical mechanism that relates the halo parameters to the luminous mass distribution.
5. A merging hierarchy causes the cooled baryonic fraction to lose angular momentum to the halo, making disks that are too small (Navarro & White 1994; Navarro & Steinmetz 1997). The predicted angular momentum of the disk is at least an order of magnitude less than that observed. The problem is only partially ameliorated (MacLow & Ferrara 1999; Navarro & Steinmetz 2000b) if some process (usually described as “feedback from star formation”) prevents most of the gas from cooling until after the galaxy is assembled. While this difficulty is best known within the CDM context, merging protogalaxies in any hierarchical structure formation model with DM will involve chronic angular momentum loss to the halos (e.g. Barnes 1992).
6. Every collapsed halo should manifest the same peak phase space density, $`f_{\mathrm{max}}`$, if DM is collisionless, was initially homogeneously distributed, and had an initially finite $`f_{\mathrm{max}}`$ (Liouville’s theorem). Further, the finite central density of galaxy halos (§2) both suggests that halos collapsed from material having in initially finite $`f_{\mathrm{max}}`$ (infinite initial phase space density forms cusped halos) and also allows $`f_{\mathrm{max}}`$ to be estimated easily. The spectacular variation of $`f_{\mathrm{max}}`$ between galaxies found by Sellwood (2000) and Dalcanton & Hogan (2000) indicates that DM cannot be a simple collisionless particle.
7. High-resolution simulations that follow the formation and evolution of individual galaxy halos in CDM find strongly cusped density profiles (Moore et al. 1998; Klypin et al. 2000) even before the baryonic component cools and settles to the center. No observational evidence requires halos to have the predicted cusps. Further, the “concentration index” has a wide range (Bullock et al. 1999), but most fits to rotation curves yield values well below the predicted range in all types of galaxy (§2) – even the Milky Way (Navarro & Steinmetz 2000a).
8. Navarro & Steinmetz (2000a) describe their failure to predict the zero-point of the TFR as a “fatal problem for the $`\mathrm{\Lambda }`$CDM paradigm.” They show that no matter what M/L is assumed for the disk, the predicted circular speed at a given luminosity is too high because the halo density is too high.
9. Simulations produce numerous sub-clumps within large DM halos (Klypin et al. 1999; Moore et al. 1999). The clumps are more numerous than the numbers of observed satellite galaxies, and may threaten the survival of a thin disk in the host galaxy.
10. The TFR discrepancy is even worse, since CDM predicts $`LV^3`$ (Dalcanton, Spergel & Summers 1997; Mo, Mao & White 1998), whereas Verheijen (1997) stresses that when $`V`$ is interpreted as the circular velocity of the flat part of the rotation curve, the true relation is very nearly $`LV^4`$. Any mechanism which systematically boosts luminosity as a function of mass must also reproduce the very small scatter in the TFR.
11. The DM halos that form in simulations are generally tri-axial (Dubinski & Carlberg 1991; Warren et al. 1992), but become nearly oblate in their inner parts when a disk is added (Dubinski 1994). Current constraints on halo shapes (Sackett 1999) are generally thought to be consistent with these predictions. However, the halo of NGC 2403 seems to become more nearly axisymmetric at larger radii (Schoenmakers 1998), opposite to the CDM expectation, and Franx et al. (1994) find IC 2006 to be impressively round at $`6R_e`$. Much more data are needed to determine whether this behavior is typical or anomalous.
12. The first precision measurements of the microwave background power spectrum at sub-degree scales show a second acoustic peak greatly suppressed compared to the first (de Bernardis et al. 2000; Hanany et al. 2000), producing an uncomfortable fit to standard cosmological models (e.g. Lange et al. 2000; Tegmark et al. 2000). McGaugh (1999b) has pointed out that unforced acoustic oscillations, as might be expected in the absence of dark matter potential wells, give such a peak height ratio.

Only three of these problems hinge on the properties mentioned in §2: If disks in bright galaxies are significantly sub-maximal, problem 1 would largely go away and 7 would be weakened (while 2 would be altered, not solved). If halos have mild density cusps, problem 6 would go away and 7 would again be weakened.

## 4. Discussion

The first five problems on our list are the most generic and intractable in DM models. Problems 1, 2 & 4 do not merely arise from mismatches between theory and observation, but require fine tuning in any DM scenario. Any physical process to account for them must somehow couple the dark and luminous components. If the only interaction between the baryons and DM is gravitational, it is hard to imagine how such a mechanism could arise.

Problems 7–10 are well-rehearsed difficulties of CDM. The possibility that 7 & 8 could be solved by a more realistic treatment of baryonic processes alone is remote since they relate to the total mass distribution which, in the simulations, is dominated by DM. Even if the baryonic mass fraction is significant, it is unclear what process could loosen tight DM concentrations; for example, Debattista & Sellwood (2000) show that an intolerable degree of dynamical friction has a very mild effect on the halo density profile. Problems 9 & 10 could conceivably be solved by some kind of complex baryonic physics (i.e. star formation), but this solution to 10 would require more fine tuning.

Another approach to problems 7–10 is to modify the assumed properties of the DM particles. Warm DM (e.g. Colombi, Dodelson & Widrow 1996; Sommer-Larsen & Dolgov 1999; Hogan 1999) is designed to aid problems 7 & 9 in two ways: WDM particles stream freely in the early universe suppressing small-scale power. They also have larger velocity spreads in halos of greater volume density, because of Liouville’s theorem, thereby precluding strong density gradients. However, fermionic WDM falls foul of problem 6, while bosonic warm DM, which has a small fraction of very high phase space density material (e.g. Madsen 2000), would probably again form cusped halos as in CDM. Self-interacting dark matter (Spergel & Steinhardt 2000) can in principle solve problems 7 & 9 (Burkert 2000; Davé et al. 2000) but only if the cross section is tuned; it may also have difficulties with its prediction of spherical halos. We emphasize that hypothetical modifications to the properties of DM particles have been motivated mostly by problems 7 & 9, and any collisional effect removes 6, but none address the much more difficult problems 1–5.

## 5. The Alternative

If the problems listed in §3 cannot be solved in the context of an acceptably simple DM model, the only logical alternative is to modify gravity in the weak-field regime. No satisfactory theory that achieves this has yet been proposed.

Remarkably, we already know the effective force-law for galaxies (McGaugh 1999a): MOND (Milgrom 1983) is a highly successful phenomenological description of the dynamical data on spiral galaxies which effortlessly resolves problems 1 through 4, 6, 10 & 11 above. Thus it disposes of those most intractable in DM models (although problem 3 is “resolved” only by introducing a new fundamental physical scale). Problem 2 was a notable prediction of MOND. It seems not unreasonable to imagine that the other problems might also be resolved once a theory of galaxy formation is formulated in some modified gravity theory.

But modifying the law of gravitation is not a panacea; in particular it faces serious challenges in reproducing the successes of standard cosmology. The following is a list of some more obvious difficulties.

1. A modified theory must be cast as a generally covariant theory to reproduce Einstein gravity in its well-tested regimes. So far, attempts to fit MOND into a natural metric framework have been largely unsuccessful (Sanders 1997).
2. It must establish a cosmology in which the scale factor evolution yields the first CMB acoustic peak at its observed angular scale.
3. Without DM, Silk damping (Silk 1968) at the time of recombination would eradicate primordial fluctuations on scales of galaxies and smaller. A new source of fluctuations to seed the growth of galaxies would therefore be required (Sanders 1998).
4. In the DM picture, the tiny fluctuations observed in the microwave background have grown to the observed current large scale structure only because the DM density variations started growing before recombination. A cosmology without DM requires other means to increase the growth factor since the time of recombination (Sanders 1998).
5. Observed gravitational lensing must arise solely from the baryonic matter distribution.

One other issue might be added to both this list of problems and that in §3. Current observational data strongly indicate a cosmological constant (e.g. Jaffe et al. 2000), while conventional physical theories offer no explanation for its tiny value (Weinberg 2000) and it requires us to live at a special time in the evolution of the universe. It is conceivable that a dramatic alteration of gravitation and cosmology might provide a different reason for accelerated expansion.

Future signatures of the absence of DM might include a still more strongly suppressed third acoustic peak in the CMB. In the absence of any driving force from dark matter, its amplitude should be further suppressed relative to the second peak due to Silk damping, while most conventional models predict an enhanced third peak from driven baryon compression. Another possible test is the imprint of these peaks on the large-scale structure power spectrum (Eisenstein & Hu 1998): without a dominant contribution of DM to the gravitational potential, substantial oscillations in the power spectrum of large-scale structure might survive.

## 6. Conclusions

Dark matter models are facing increasingly serious difficulties on galaxy scales as the quality of both the simulations and the data constraining the properties of galaxy halos accumulates. Those who suggest that a revised law of gravity in the ultra-weak field limit might be more natural (Milgrom 1983; Sanders 1997; McGaugh & de Blok 1998; Mannheim 2000; and others) are regarded as radicals. However, theories incorporating both a fine-tuned cosmological constant and dark matter particles with ad hoc properties have lost the compelling aesthetic simplicity of the original CDM model.

The problem of formulating a generally covariant, alternative theory of gravity upon which a competing model of structure formation can be built is itself a severe challenge. Without such a theory, we cannot say whether the other problems on our preliminary list in §5 will turn out to be harder to resolve than those of §3, and it is possible that the list would grow to become still more daunting if such a theory could be developed. But the difficulties facing DM models and the success of the MOND phenomenology suggest that a serious examination of alternative gravity solutions is called for.

### Acknowledgments.

We thank Paul Steinhardt and Michael Turner for critical readings of a draft version. Comments from Stacy McGaugh, Bob Sanders, David Spergel, Matthias Steinmetz and Scott Tremaine were also helpful. This work was supported by NASA LTSA grant NAG 5-6037 and NASA ATP grant NAG 5-7015. AK is a Cotrell Scholar of the Research Corporation.

## References

Bahcall, J. N. & Casertano, S. 1986, ApJ, 308, 347

Bahcall, N., Ostriker, J. P., Perlmutter, S. & Steinhardt, P. J. 1999, Science, 284, 1481

Barnes, J. E. 1992, ApJ, 393, 484

Blais-Ouellette, S, Carignan, C., Amram, P. & Côté, S. 2000a, AJ, 118, 2123

Blais-Ouellette, S, Amram, P. & Carignan, C. 2000b, astro-ph/0006449

Blumenthal, G. R., Faber, S. M., Flores, R. & Primack, J. R. 1986, ApJ, 301, 27

Bosma, A. 1999, in Galaxy Dynamics, eds. D. Merritt, J. A. Sellwood & M. Valluri, ASP Conference Series 182, (San Francisco: Astronomical Society of the Pacific) p. 339

Broeils, A. H. & Courteau, S. 1997, in Dark and Visible Matter in Galaxies, eds. M. Persic & P. Salucci, ASP Conference Series 117, (San Francisco: Astronomical Society of the Pacific) p. 74

Bullock, J. S. et al. 1999, astro-ph/9908159

Burkert, A. 2000, ApJ, to appear (astro-ph/0002409)

Casertano, S. & van Gorkom, J. H. 1991, AJ, 101, 1231

Colombi, S., Dodelson, S. & Widrow, L. M. 1996, ApJ, 458, 1

Courteau, S. 1997, AJ, 114, 2402

Dalcanton, J. & Hogan, C. J. 2000, astro-ph/0004381

Dalcanton, J., Spergel, D. N. & Summers, J. J. 1997, ApJ, 482, 659

Davé, R., Spergel, D. N., Steinhardt, P. J. & Wandelt, B. D. 2000, ApJ, submitted (astro-ph/0006218)

de Bernardis, P., et al. 2000, Nature, 404, 955

Debattista, V. P. & Sellwood, J. A. 1998, ApJ, 493, L5

Debattista, V. P. & Sellwood, J. A. 2000, ApJ, 544, to appear (astro-ph/0006275)

de Blok, W. J. G. & McGaugh, S. S. 1997, MNRAS, 290, 533

Dubinski, J. 1994, ApJ, 431, 617

Dubinski, J. & Carlberg, R. G. 1991, ApJ, 378, 496

Eisenstein, D. & Hu, W. 1998, ApJ, 496, 605

Franx, M., van Gorkom, J. H. & de Zeeuw, T. 1994, ApJ, 436, 642

Hanany, S., et al. 2000, astro-ph/0005123

Hogan, C. J. 1999, astro-ph/9912549

Jaffe, A. H., et al. 2000, astro-ph/0007333

Kalnajs, A. J. 1983, in Internal Kinematics and Dynamics of Galaxies IAU Symposium 100, ed. E. Athanassoula (Dordrecht: Reidel) p. 87

Kent, S. M. 1986, AJ, 91, 1301

Klypin, A., Kravtsov, A. V., Valenzuela, O. & Prada, F. 1999, ApJ, 522, 82

Klypin, A., Kravtsov, A. V., Bullock, J. S. & Primack, J. R. 2000, astro-ph/0006343

Lange, A. E., et al. 2000, astro-ph/0005004

MacLow, M-M. & Ferrara, A. 1999, ApJ, 513, 142

Madsen, J. 2000, astro-ph/0006074

Mannheim, P. D. 2000, gr-qc/0001011

McGaugh, S. 1999a, in Galaxy Dynamics, eds. D. Merritt, J. A. Sellwood & M. Valluri, ASP Conference Series 182, (San Francisco: Astronomical Society of the Pacific) p. 528

McGaugh, S. 1999b, ApJ, 523, L99

McGaugh, S. S. & de Blok, W. J. G. 1998, ApJ, 499, 41 & 66

McGaugh, S. S., Schombert, J. M., Bothun, G. D. & de Blok, W. J. G. 2000, ApJ, 533, 99

Milgrom, M. 1983, ApJ, 270, 365

Milgrom, M. 1998, astro-ph/9810302

Mo, H. J., Mao, S. & White, S. D. M. 1998, MNRAS, 295, 319

Moore, B., Governato, F., Quinn, T., Stadel, J. & Lake, G. 1998, ApJ, 499, L5

Moore, B., Ghigna, S., Governato, F. & Lake, G. 1999, ApJ, 524, L19

Navarro, J. F. & Steinmetz, M. 1997, ApJ, 478, 13

Navarro, J. F. & Steinmetz, M. 2000a, ApJ, 528, 607

Navarro, J. F. & Steinmetz, M. 2000b, ApJ, 538, 477

Navarro, J. F. & White, S. D. M. 1994, MNRAS, 267, 401

Palunas, P. & Williams, T. B. 2000, AJ, (to appear)

Rubin, V. C., Burstein, D., Ford, W. K., Jr. & Thonnard, N. 1985, ApJ, 289, 81

Sackett, P. D. 1999, in Galaxy Dynamics, eds. D. Merritt, J. A. Sellwood & M. Valluri, ASP Conference Series 182, (San Francisco: Astronomical Society of the Pacific) p. 393

Sanders, R. H. 1990, A&A Rev., 2, 1

Sanders, R. H. 1997, ApJ, 480, 492

Sanders, R. H. 1998, MNRAS, 296, 1009

Sanders, R. H. & Verheijen, M. A. W. 1998, ApJ, 503, 97

Schoenmakers, R. 1998, in Galactic Halos, ed. D. Zaritsky, ASP Conference Series 136, (San Francisco: Astronomical Society of the Pacific) p. 240

Sellwood, J. A. 2000, ApJ, 540, L1

Silk, J. 1968, ApJ, 151, 459

Sofue, Y., Tutui, Y., Honma, A., Takamiya, T., Koda, J., & Takeda, Y. 1999, ApJ, 523, 136

Sommer-Larsen, J. & Dolgov, A. 1999, astro-ph/9912166

Spergel, D. N. & Steinhardt, P. J. 1999, astro-ph/9909386

Sprayberry, D., Bernstein, G. M., Impey, C. D. & Bothun, G. D. 1995, ApJ, 438, 72

Swaters, R. A., Madore, B. F. & Trewhella, M. 2000, ApJ, 531, L107

Tegmark, M., Zaldarriaga, M. & Hamilton, A. J. S. 2000, astro-ph/0008167

Tytler, D., O’Meara, J. M., Suzuki, N. & Lubin, D. 2000, Physica Scripta, T85, 12 (astro-ph/0001318)

van den Bosch, F. C, Robertson, B. E., Dalcanton, J. J. & de Blok, W. J. G, 2000, AJ, 119, 1579

Verheijen, M. A. W. 1997, Ph.D. thesis, University of Groningen

Warren, M. S., Quinn, P. J., Salmon, J. K. & Zurek, W. H. 1992, ApJ, 399, 405

Weinberg, S. 2000, astro-ph/0005265

Weiner, B. J., Sellwood, J. A. & Williams, T. B. 2000, ApJ, to appear, astro-ph/0008205

Zwaan, M. A., van der Hulst, J. M., de Blok, W. J. G. & McGaugh, S. S. 1995, MNRAS, 273, L35


text:
# Resolving the Mystery of X-ray Faint Elliptical Galaxies: Chandra X-ray Observations of NGC 4697

## 1 Introduction

Observations with the Einstein X-ray Observatory showed that elliptical and S0 galaxies are luminous sources of X-ray emission (e.g., Forman, Jones, & Tucker 1985). At least for the X-ray luminous early-type galaxies (defined as those having a relatively high ratio of X-ray to optical luminosity $`L_X/L_B`$), the bulk of the X-ray luminosity is from hot ($``$$`10^7`$ K) interstellar gas. However, there is a large dispersion in the X-ray luminosities of early-type galaxies of a given optical luminosity. We will refer to galaxies which have a very low $`L_X/L_B`$ ratio as “X-ray faint.” In these X-ray faint ellipticals, much of the hot interstellar gas may have been lost in galactic winds or by ram pressure stripping by ambient intracluster or intragroup gas. The source of the bulk of the X-ray emission in X-ray faint galaxies is uncertain.

In general, X-ray faint galaxies exhibit significantly different spectral properties than their X-ray bright counterparts. The X-ray bright galaxies are dominated by thermal emission at $`kT0.8`$ keV due to their hot interstellar medium (ISM). On the other hand, the X-ray faint galaxies exhibit two distinct spectral components: first, they have a hard $``$5–10 keV component, most easily seen in ASCA spectra (Matsumoto et al. 1997), which is roughly proportional to the optical luminosity of the galaxy. This suggests that the hard component is due to low-mass X-ray binaries (LMXBs) like those seen in the bulge of our Galaxy and M31. X-ray faint galaxies also have a very soft ($``$0.2 keV) component, whose origin is uncertain (Fabbiano, Kim, & Trinchieri 1994; Pellegrini 1994; Kim et al. 1996). Suggested stellar sources for the very soft emission in X-ray faint ellipticals include active M stars, RS CVn binaries, or supersoft sources, but none of these appears to work quantitatively (Pellegrini & Fabbiano 1994; Irwin & Sarazin 1998b). It is possible that the soft X-rays are due to warm (0.2 keV) ISM. Recently, we proposed that the very soft emission in X-ray faint ellipticals is due to the same LMXBs responsible for the hard emission (Irwin & Sarazin 1998a,b). However, the origin of the very soft component in X-ray faint ellipticals remains something of a mystery.

With the superb spatial resolution of the Chandra X-ray Observatory, it should be possible to resolve the emission from nearby X-ray faint early-type galaxies into LMXBs, if this is indeed the source of their emission. Recently, we used a deep ROSAT HRI observation to detect a number of discrete X-ray sources in the X-ray faint elliptical NGC 4697 (Irwin, Sarazin, & Bregman 2000, hereafter ISB). However, the bulk of the X-ray emission was not resolved. We also simulated a 40 ksec Chandra observation of NGC 4697, and showed that it should be possible to detect $``$100 LMXBs if they provide the bulk of the emission. Here, we present the results of exactly this observation. At a distance of 15.9 Mpc (Faber et. al. 1989; assuming a Hubble constant of 50 km s-1 Mpc-1), NGC 4697 is the closest normal, optically luminous, X-ray faint elliptical galaxy. Given its proximity, NGC 4697 is an ideal target for detecting the LMXB population. It is sufficiently X-ray faint that diffuse ISM emission should not bury the fainter LMXBs. The purposes of the Chandra observation are to resolve and study the LMXB population of NGC 4697, to determine the source of both the hard and soft spectral components, and to detect or place strong limits on any residual diffuse (possibly gaseous) emission.

## 2 X-ray Observation

NGC 4697 was observed on 2000 January 15-16 with the Chandra ACIS-S3 detector for a total exposure of 39,434 s. The S3 chip was chosen because of its good low energy response. No background flares occurred during this observation. This observation was processed at a time when the standard pipeline processing introduced a boresight error of about 8″ in the absolute positions of X-ray sources. The absolute positions were corrected by using the optical identifications and positions from USNO-A2.0 (Monet et al. 1998) of X-ray sources. The absolute positions are accurate to about 0.5″ near the center of the S3 image, with larger errors further out. Background was taken from the outer regions of the same chip and blank sky backgrounds. More details concerning the X-ray observation, data analysis, variability of the sources, and spectral fitting are given in Sarazin, Irwin, & Bregman (2000, hereafter Paper II).

## 3 X-ray Image

The Chandra S3 chip X-ray image is shown in Figure 2. The central portion of the Chandra image (adaptively smoothed) is shown in Figure 3. These images show the basic result of the Chandra observation: much of the emission from the galaxy is resolved into individual point sources of X-rays (§ 4). However, Figure 3 indicates that there is also some unresolved, diffuse emission. The diffuse emission is centrally concentrated, and its X-ray surface brightness increases by about a factor of 30 into the center of the galaxy. Both the point sources and the diffuse emission are elongated in the same direction as this E6 galaxy (a position angle of 67).

We determined the portion of the X-ray emission due to resolved sources and diffuse emission, both for the entire Chandra  X-ray band (0.3–10 keV) and for three narrower bands: hard H (2-10 keV), medium M (1-2 keV), and soft S (0.3-1 keV). Here, we report on the counts coming from within the elliptical optical isophote containing one half of the optical light (“one effective radius,” semimajor axis $`95^{\prime \prime }`$, semiminor axis $`55^{\prime \prime }`$). First, the population of resolved sources was determined as discussed below (§ 4). Then, the total emission was determined, and the source flux subtracted to give the amount of unresolved emission.

Within one effective radius, 61% of the X-ray emission is resolved into individual X-ray sources for the total band. In the hard, medium, and soft bands, the resolved fractions are 81%, 79%, and 48%, respectively. Thus, it is clear that the bulk of the total X-ray emission from NGC 4697 comes from point sources, which are likely to be LMXBs.

A portion of the unresolved emission must also come from LMXBs which fall below our threshold for reliable detection of resolved sources. Indeed, if one sets the detection threshold for sources lower, one finds many more 2-$`\sigma `$ fluctuations than expected just from Poisson statistics. If the observed luminosity function of the resolved sources in NGC 4697 (Fig. 4 below) is extended down to $`L_X10^{36}`$ ergs s-1, the approximate lower limit for LMXBs in globular clusters (Hertz & Grindlay 1983), the contribution of LMXBs to the total band X-ray emission from NGC 4697 increases to 74% within one effective radius. The X-ray spectrum of the diffuse emission is somewhat softer than that of the resolved sources. If one decomposes the unresolved emission into two components, one with the same spectrum as that of the sum of the resolved sources plus a softer component, one finds that the LMXBs would contribute about 77% of the total band emission.

The remaining $``$23% of the emission would come from a more diffuse component with a soft ($``$0.3 keV) spectrum. This soft component has a very extended spatial distribution (Paper II), and is most likely due to interstellar gas. Additional evidence for ISM emission comes from the adaptively smoothed X-ray image (Fig. 3). If the diffuse X-ray emission is from a large number of faint stellar sources, then it should appear smooth and symmetrical (like the optical image). The diffuse emission in Figure 3 is smooth and fairly symmetric in the inner regions, but there is an extended asymmetrical feature at the outer regions to the northeast. Such asymmetries are more easily produced in the distribution of gas (e.g., by ram pressure) than in the distribution of stellar sources.

## 4 Resolved Sources

The discrete X-ray source population on the ACIS S3 image (Fig. 2) was determined, using a wavelet detection algorithm with the source detection threshold set at $`10^6`$, which implies that $``$1 false source (due to a statistical fluctuation in the background) would be detected in the entire S3 image. The minimum detectable flux was about $`2.6\times 10^4`$ cts s-1 ($`L_X=5.0\times 10^{37}`$ ergs s-1 at the NGC 4697 distance) in the 0.3–10 keV band in the part of the chip covered by the galaxy. Fluxes were corrected for exposure and the instrument PSF. A total of 90 sources were found; a complete list with positions and other properties is given in Paper II.

Some of the detected sources are likely to be unrelated foreground or (more likely) background objects. Based on the source counts in Brandt et al. (2000) and Mushotzky et al. (2000), we would expect about 10–15 unrelated sources in our observation. The unrelated sources should mainly be found at larger distances from the optical center of NGC 4697, while the sources associated with NGC 4697 should be concentrated to the center of the galaxy. Within 2′ of the center of NGC 4697 (roughly the region covered by Fig. 3), $``$3 of the $``$60 detected sources would be expected to be unrelated to NGC 4697.

Eight of the X-ray sources coincide (within 1″) with candidate globular clusters in NGC 4697. CXOU J124834.4-055014 (Paper II, Src. 64) is identified with globular cluster 33 in Hanes (1977), while CXOU J124846.8-054852 (Paper II, Src. 72, ISB Src. 11) is associated with Hanes globular cluster 24. Six additional X-ray sources (CXOU J124827.0-054925, CXOU J124826.1-054729, CXOU J124830.8-054836, CXOU J124834.2-054926, CXOU J124835.9-054551, CXOU J124841.5-054736) agree with the positions of candidate globulars kindly provided by Kavelaars (2000). Given the positional agreement of the X-ray sources with the globular candidates, it is likely that all 8 of these identifications are real. However, at the distance to NGC 4697, globular clusters are not resolved in ground-based optical images, and the candidate globulars were identified by luminosities and (possibly) colors. As a result, as many as half of them might be unrelated faint optical objects, rather than globular clusters.

The count rates for the sources were converted into unabsorbed luminosities (0.3-10 keV) assuming that all of the sources were at the distance of NGC 4697. Using the best-fit Chandra X-ray spectrum of the resolved sources (Paper II), the resulting X-ray luminosities range from about $`5\times 10^{37}`$ to $`2.5\times 10^{39}`$ ergs s-1 (calibration uncertainties might affect the luminosities by $``$25%). The cumulative luminosity function of all of the sources is shown in Figure 4. Note that the luminosity function appears to have a knee at $`L_X3\times 10^{38}`$ ergs s-1. In fitting the luminosity function, we corrected for unrelated sources using the deep source counts (Brandt et al. 2000; Mushotzky et al. (2000), assuming that the shape of the number versus flux relation is that derived from ROSAT deep fields (Hasinger et al. 1998). A single power-law (plus the background contribution) fit to the luminosity function could be rejected at the $`>`$95% confidence level. However, a broken power-law fit was successful:

$$\frac{dN}{dL_{38}}=N_o\left(\frac{L_{38}}{L_b}\right)^\alpha ,$$
(1)
with $`\alpha =\alpha _l`$ for $`L_{38}L_b`$ and $`\alpha =\alpha _h`$ for $`L_{38}>L_b`$. Here, $`L_{38}`$ is the X-ray luminosity (0.3–10 keV) in units of $`10^{38}`$ ergs s-1. The best fit, determined by the maximum-likelihood method, gave $`N_o=8.0_{5.2}^{+7.1}`$, $`\alpha _l=1.29_{0.49}^{+0.36}`$, $`\alpha _h=2.76_{0.39}^{+1.81}`$, and a break luminosity of $`L_b=3.2_{0.8}^{+2.0}\times 10^{38}`$ ergs s-1 (Figure 4). The errors (at the 90% confidence level) were determined by Monte Carlo simulations.

The break luminosity is similar to the Eddington luminosity for spherical accretion onto a 1.4 $`M_{}`$ neutron star. This suggests that the sources with luminosities above this break are accreting black holes, while those below the break are predominantly neutrons stars. This would imply that NGC 4697 contains $``$15 luminous X-ray binary systems containing black holes. If the more luminous of these systems ($`L_X10^{39}`$ ergs s-1) are limited by the Eddington luminosity, they must contain fairly massive ($`M8M_{}`$) black holes.

We studied the crude spectral properties of the resolved sources by using hardness ratios (full results and spectral analysis in Paper II). Hardness ratios or X-ray colors have the advantage that they can be applied to weaker sources. Figure 4 plots H31 vs. H21. For comparison, the hardness ratios (H21,H31) are $`(0.38,0.57)`$ for all of the emission, $`(0.69,0.82)`$ for the unresolved emission, and $`(0.14,0.37)`$ for the sum of the sources, all within one effective radius.

There are three moderate luminosity sources with hardness ratios of $`(1,1)`$, which means that they have no detectable emission beyond 1 keV. These three sources are almost certainly supersoft sources (e.g., Kahabka & van den Heuvel 1997). There are three sources, one of which is the brightest source in the field, with values of $``$$`(1,1)`$, all of which are located $`>`$3$`\stackrel{}{\mathrm{.}}`$5 from the center of NGC 4697. These are probably unrelated, strongly absorbed AGNs. There are eight sources near $`(0,1)`$ which have essentially no hard emission. Six of these sources are at large radii ($`>`$2$`\stackrel{}{\mathrm{.}}`$2), which suggests that this population is also unrelated to NGC 4697. Most of the sources lie in a diagonal swath centered at about $`(0.15,0.40)`$. These values are similar to but slightly harder than the integrated colors for the entire galaxy, but considerably harder than the values for the unresolved emission. The X-ray sources identified with globular cluster candidates have hardness ratios which are similar to other LMXBs, with the exception of CXOU J124834.4-055014, whose hardness ratios suggest that it is actually a background AGN.

## 5 Conclusions

Our high spatial resolution Chandra observation of the X-ray faint elliptical galaxy NGC 4697 resolves most of the X-ray emission (61% within one effective radius) into point sources. A total of 90 individual sources are detected, of which $``$80 are low mass X-ray binaries (LMXBs) associated with NGC 4697. Much of the emission is resolved even in the softest band. The dominance of LMXBs implies that this and other X-ray faint ellipticals have lost most of their interstellar gas.

Of the remaining unresolved emission, it is likely that about half is from fainter LMXBs, while about 23% of the total emission is probably from interstellar gas with $`kT0.3`$ keV. Faint, diffuse X-ray emission in the Chandra image is asymmetrical and very extended, which also suggests that a portion of the emission is due to hot gas. Ram pressure may have produced the asymmetric feature.

Three of the resolved sources appear to be supersoft sources. Eight of the resolved sources in the outer parts of NGC 4697 (about 25%) are coincident with candidate globular clusters. On the other hand, the candidate globular clusters contain about 0.1% of the optical light of the galaxy in this region. This indicates that globulars have an unusually large fraction of X-ray binaries (by a factor of $``$250), as is seen in our Galaxy (e.g., Hertz & Grindlay 1983).

The X-ray luminosities (0.3–10 keV) of the resolved LMXBs range from $``$$`5\times 10^{37}`$ to $``$$`2.5\times 10^{39}`$ ergs s-1. The luminosity function has a “knee” at $`3.2\times 10^{38}`$ ergs s-1, which is roughly the Eddington luminosity of a 1.4 $`M_{}`$ neutron star (NS). This knee might provide a standard candle which could be used to determine distances to galaxies. This knee may separate accreting NS and black hole (BH) binaries. If the brightest sources in NGC 4697 are Eddington limited, they must contain fairly massive BHs.

Our detection of a large population of binaries with NSs and massive BHs provides perhaps the most direct evidence that this elliptical galaxy (or its progenitors) once contained a large number of massive main sequence stars. The population of LMXBs provides a tool to study the high mass end of the initial mass function of early-type galaxies, long after the massive main sequence stars have died. The detection of NS and massive black holes in NGC 4697 also provides the first evidence for the current existence of massive stars (albeit degenerate) in a normal, optically luminous elliptical galaxy.

We are extremely grateful to JJ Kavelaars for providing his unpublished list of globular clusters in NGC 4697. Bill Harris, JJ Kavelaars, and Arunav Kundu also gave helpful comments on the globular cluster population of NGC 4697. Support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Number GO0-1019X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. J. A. I. was supported by Chandra Fellowship grant PF9-10009, awarded through the Chandra Science Center.