text:
# Temperature Variations and Abundance Determinations in Planetary Nebulae

## 1. Overview

Photoionization models for chemically homogeneous gaseous nebulae of constant density predict an almost constant temperature and consequently very often observers assume a constant temperature to determine chemical abundances. There is growing observational evidence that indicates the presence of large temperature variations in gaseous nebulae, in contradiction with the models mentioned above. Consequently the relevance of temperature variations in the abundance determinations has to be considered in detail.

In section 2 we review the observational evidence in favor of large temperature variations. In section 3 we review different possible explanations for the observed variations. In section 4 we discuss the implications of the temperature variations in the abundance determinations, these implications depend on the source of the temperature variations.

Recent reviews on the temperature structure of gaseous nebulae are those of Peimbert, M. (1995, 2002), Esteban (2002), Stasińska (2002), and Liu (2002a, 2002b).

## 2. Temperature Variations

### 2.1. Definitions

The average temperature, $`T_0`$, and the mean square temperature fluctuation, $`t^2`$, are given by

$$T_0(N_e,N_i)=\frac{T_e(𝐫)N_e(𝐫)N_i(𝐫)𝑑V}{N_e(𝐫)N_i(𝐫)𝑑V},$$
(1)
and

$$t^2=\frac{(T_eT_0)^2N_eN_i𝑑V}{T_0^2N_eN_i𝑑V},$$
(2)
respectively, where $`N_e`$ and $`N_i`$ are the electron and the ion densities of the observed emission line and $`V`$ is the observed volume (Peimbert 1967).

To determine $`T_0`$ and $`t^2`$ we need two different methods to derive $`T_e`$: one that weighs preferentially the high temperature regions and one that weighs preferentially the low temperature regions. For example the temperature derived from the ratio of the \[O III\] $`\lambda \lambda `$ 4363, 5007 lines, $`T_{(4363/5007)}`$, and the temperature derived from the ratio of the Balmer continuum to $`I(H\beta )`$, $`T_{(\mathrm{Bac}/H\beta )}`$, that are given by

$$T_{(4363/5007)}=T_0\left[1+\frac{1}{2}\left(\frac{90800}{T_0}3\right)t^2\right],$$
(3)
and

$$T_{(\mathrm{Bac}/\mathrm{H}\beta )}=T(\mathrm{Bac})=T_0(11.70t^2),$$
(4)
respectively. It is also possible to use the intensity ratio of a collisionally excited line of an element $`p+1`$ times ionized to a recombination line of the same element $`p`$ times ionized, this ratio is independent of the element abundance and in the low density limit depends only on the electron temperature. Two examples of this method are:

$$T_{(\mathrm{O}\mathrm{II}rec/\mathrm{O}\mathrm{III}coll)}=T_{(4649/5007)}=f_1(T_0,t^2)$$
(5)
and

$$T_{(\mathrm{C}\mathrm{II}rec/\mathrm{C}\mathrm{III}coll)}=T_{(4267/1909)}=f_2(T_0,t^2),$$
(6)
by combining these ratios with a temperature determined from the ratio of two collisionally excited lines, like $`T_{(4363/5007)}`$, it is also possible to derive $`T_0`$ and $`t^2`$.

From the ratio of two He I recombination lines it is possible to obtain relationships of the type:

$$T_{(\mathrm{He}\mathrm{I},i/j)}=f_{ij}(T_0,t^2,N_e,\tau _{3889}),$$
(7)
and from a set of 5 or more helium lines (which should include $`\lambda \lambda `$ 3889, 6678, and 7065) it is possible to obtain $`T_0`$, $`t^2`$, $`N_e`$, and $`\tau _{3889}`$ (see Peimbert, Peimbert, & Ruiz 2000). Sometimes the errors in these determinations are quite large, in such case it is possible to combine the information from the helium lines with the one from $`T_{(4363/5007)}`$ to obtain better determinations of $`T_0`$ and $`t^2`$ (see Peimbert, Peimbert, & Luridiana 2002).

### 2.2. Errors in the temperature determinations

Different types of errors are present in the determinations of the electron temperatures, and in some cases these errors have led to erroneous temperature determinations. Before accepting that temperature variations are present in a given object, it is necessary to rule out any possible sources of error.

Rola & Stasińska (1994) and Rola & Pelat (1994) have argued that there is a systematic overestimation of the line intensities of lines with low signal to noise ratios that leads to the overestimation of the $`T_{(4267/1909)}`$ temperatures. This effect is present in many objects and has to be taken into consideration, but for many well observed objects it is not the main cause of the derived large $`T_{(4363/5007)}`$$`T_{(4267/1909)}`$ differences (e.g., Peimbert, Torres-Peimbert, & Luridiana 1995).

The discrepancy between the C abundances derived from the C II $`\lambda `$ 4267 recombination line and those derived from the collisionally excited 1906+1909 C III lines has been known for at least 20 years and the idea that the large intensity of the C II $`\lambda `$ 4267 line was due to an incomplete knowledge of the relevant atomic physics has been mentioned on several occasions. Nevertheless similar excesses in the intensities of recombination lines are observed in different multiplets of the same ion and in other elements (Liu et al. 1995; Mathis, Torres-Peimbert, & Peimbert 1998).

Viegas & Clegg (1994) suggested that the differences between $`T_{(4649/5007)}`$ and $`T`$(Bac) could be due to density fluctuations diminishing the $`I(4959)`$ and $`I(5007)`$ values by collisional de-excitation yielding a spuriously high $`T_{(4649/5007)}`$ value. This is indeed the case for the planetary nebula M2-29, where from ground based observations of the whole object a $`T_{(4363/5007)}`$ of 24,000 K and a $`N_e=2000`$ cm-3 are derived, while from HST observations it is found that it has at least two density components (one with $`N_e10^4`$ cm-3 and another with $`N_e10^6`$ cm-3, the $`T_e`$ in the low density component is of 8800 K); the high temperature derived from the ground observations is spurious and yields O/H$`=5.1\times 10^5`$, a factor of ten smaller than the value derived adopting a two density component model (Peña, Torres-Peimbert, & Ruiz 1991; Torres-Peimbert et al. 2002). A similar result has been found for Mz 3 by Zhang & Liu (2002).

### 2.3. Observations of Temperature Variations

From uniform density chemically homogeneous photoionization models it has been found that $`t^2`$ is in the 0.000 to 0.025 range, with typical values around 0.005 (e.g., Gruenwald & Viegas 1995).

From observations and by combining equation 3 with any of the equations 4 to 7 it is possible to determine $`t^2`$ and $`T_0`$. These combinations indicate values of $`t^2`$ in the 0.00 to 0.15 range with typical values around 0.04. The errors in the best $`t^2`$ determinations given in the literature are typically in the 0.01 to 0.02 range.

Some of the best determinations of $`t^2`$ and $`T_0`$ are: from equations 3 and 4 those by Peimbert (1971), Liu & Danziger (1993), Liu et al. (1995, 2000, 2001); from equations 3 and 5 those by Peimbert, Storey, & Torres-Peimbert (1993), Liu et al. (1995, 2000, 2001); from equations 3 and 6 those by Rola & Stasińska (1994), Peimbert, et al. (1995); and from equations 3 and 7 those by Peimbert et al. (1995, 2000, 2002).

Liu et al.(2001, see also Peimbert, A. 2002), presented a good correlation between $`Log(\mathrm{O}^{++}/\mathrm{H}^+)_{RL}/(\mathrm{O}^{++}/\mathrm{H}^+)_{CEL}`$ and $`T_{(4363/5007)}`$$`T`$(Bac) and mention that this correlation strongly supports the idea that the temperature variations are real. In Table 1 we present the $`T_{(4363/5007)}`$ and $`T(\mathrm{Bac})`$ temperatures derived by Liu et al. (2000, 2001) for three of the objects that show very large temperature variations. From these values we have computed $`T_{(4649/5007)}`$ and from equations 3 and 4 we have derived $`t^2`$(3,4), and from equations 3 and 5 we have derived $`t^2`$(3,5). From this table it can be seen that $`t^2`$(3,4) and $`t^2`$(3,5) are not only proportional to each other but that they are in good agreement with each other. This implies that most of the temperature differences among $`T(\mathrm{Bac})`$, $`T_{(4363/5007)}`$, and $`T_{(4649/5007)}`$ are due to real temperature variations in the nebulae, thus an explanation for the origin of these variations has to be sought.

## 3. What Causes the Temperature Variations?

In what follows we present seven different mechanisms that can produce large temperature variations. It is clear that all of them should be considered. Their relative importance probably varies for each object, and in many cases, the ones that show $`t^2<0.01`$ values, all of them are negligible.

### 3.1. Shadowed Regions

To explain the intensities of lines of low degree of ionization Mathis (1976) suggested the existence of regions shadowed from the central stars, these regions would be ionized by diffuse radiation and would have a considerably smaller temperature than those regions ionized directly from the central stars of H II regions.

### 3.2. Chemical Inhomogeneities

Some PNe, like A30 and A78, present important chemical inhomogeneities; therefore the idea that temperature variations can be due to various degrees of chemical inhomogeneities has been proposed. Chemically inhomogeneous photoionization models have been presented by different authors (e.g., Torres-Peimbert, Peimbert, & Peña 1990; Kingdon & Ferland 1998; Pequignot et al. 2002a, 2002b; Tylenda 2002).

### 3.3. Density Variations

Extreme density variations are present in most PNe, as can be seen from optical images; these clumps, shells, and filaments, can produce temperature variations by themselves under the assumption of photoionization heating only. Mihalszky & Ferland (1983) have presented chemically homogeneous photoionization models, including density variations in the $`10^2`$ to $`10^{4.3}\mathrm{cm}^3`$ range, and find values of $`t^2`$ of about 0.005, which could explain only a small fraction of the largest discrepancies observed.

### 3.4. Deposition of Mechanical Energy

Shocks have been proposed as a mechanism to account for temperature fluctuations (Peimbert, Sarmiento, & Fierro 1991; Peimbert et al. 1995 and references therein). Peimbert et al. (1995) have found that those objects that show higher velocity dispersions are the ones that show higher temperature variations, supporting the presence of shocks. In a PN a shock front is expected to be formed at the interface between a fast stellar wind from the central star and the slowly expanding circumstellar envelope. The fast stellar wind would carve out a hot bubble in the interior of the nebula that would emit primarily in X-rays. Clear evidence for such hot bubbles has been found recently for a small number of PNe from X-ray imaging (Kastner et al. 2000; Kastner, Vrtilek, & Soker 2001; Chu et al. 2001).

Mellema (1997) finds strong temperature variations from hydrodynamical simulations of bipolar PNe. The dense gas near the equator has temperatures around 9000 K while the gas near the poles has temperatures between 20,000 and 50,000 K. The reasons are two: at the poles the nebula has a higher shock velocity, leading to more heating; and a lower density, leading to less efficient cooling. Bohigas (2002) finds that in the outer regions of 9 PNe of Type I the $`T(\mathrm{O}^{++})/T(\mathrm{N}^+)`$ is large and the H$`\alpha `$/\[S II\] line ratio is small, a typical combination in shock excited plasmas.

### 3.5. Deposition of Magnetic Energy

Ferland, in the meeting on Ionized Gaseous Nebulae (Mexico City 2000), suggested that magnetic reconnection can provide localized temperature variations in ionized plasmas. According to García-Segura, López, & Franco (2001) the multiple, regularly spaced concentric shells around some PNe could be due to the effects a of a solar-like magnetic cycle, with periodic polarity inversion, in the slow wind of an asymptotic giant branch (AGB) star. Presumably these shells of alternating polarity could give rise to magnetic reconnection processes once that they are compressed in the formed PN (i.e. in the swept-up shell). Magnetic fields of the order of milligauss have been measured in the torus surrounding K3-35 (Miranda et al. 2002); also the presence of a magnetic field has been inferred in the planetary nebula OH0.9+1.3 (Zijlstra et al. 1989) and in the pre-planetary nebula IRAS17150$``$3224 (Hu et al. 1993).

### 3.6. Dust Heating

Stasińska & Szczerba (2001) have analyzed the effects of photoelectric heating by dust grains in photoionization models of PNe. They have shown that this process is particularly important if filamentary PNe contain a population of small grains. The temperature structure of such dusty and filamentary nebulae would produce large temperature variations.

### 3.7. Decrease of the Ionizing Flux

The ionizing flux from the central stars of PNe is not constant. As a first approximation the ionization increases, then reaches a maximum and finally decreases. During the decrease and in those directions where the nebula is density bounded there would be ionized regions that would become isolated from the stellar ionizing photons, these regions would cool first and would recombine later (e.g., Spitzer 1978). A typical recombination time is $`10^5/N_e`$ years while the cooling time is an order of magnitude smaller, therefore cool regions with a relatively high degree of ionization are expected to be present.

This effect might be important for PNe produced by massive stars, because the central stars evolve relatively fast. One possible example of this mechanism might be provided by N66 in the LMC. In this object the ionizing flux has diminished by a factor of two in the last ten years (Peña et al. 1997; Peña, Hamann, & Ruiz 2002) and the RL/CEL abundance ratio is of about an order of magnitude (Tsamis et al. 2002). Other PNe central stars have presented important flux changes in timescales of a few years (e.g., Kostiakova 2002 and references therein). It is possible that this mechanism might account for the very low values of $`T`$(Bac), $``$3000 K, found by Luo & Liu (2002) in the outer parts of NGC 7009.

It would be possible also to isolate regions from the stellar ionizing photons if the amount of material in the line of sight increases with time, even if the stellar ionizing radiation remains constant. It has been suggested that dust clouds moving tangentially to the line of sight are responsible for the 1981–1985 and 1996–1997 fadings of V651 Mon, the central star of NGC 2346 (Costero et al. 1986; Kato, Nogami, & Baba 2001). In addition to dust clouds a similar effect would be caused by tangentially moving gaseous clouds of higher density than the ambient density.

## 4. Discussion

Temperature variations higher than those predicted by chemically homogeneous photoionization models of constant density are definitively present in a good fraction of PNe.

Chemically homogeneous photoionization models of constant density provide a good first approximation to the temperature structure and to the chemical abundances only for about half of the well observed nebulae. However, for the other half, it is paramount to have observations of high accuracy to determine their temperature structure which is needed to obtain accurate abundances of heavy elements relative to hydrogen.

Spatially resolved values of $`T_0`$ and $`t^2`$ are needed to decide among the various possibilities mentioned above. For example, if $`T_0`$ decreases in the outer regions, the temperature variations might be due to a decrease of the ionizing flux; in the opposite case the temperature variations could be due to shock heating.

A very important observational clue to study this problem is provided by Garnett & Dinerstein (2001) who find a large range in the difference between the O++ abundances derived from O II and those derived from \[O III\], spanning from no difference up to a factor of six. The size of this discrepancy is anti-correlated with nebular surface brightness; compact, high-surface-brightness nebulae have the smallest discrepancies.

It is important to determine the cause of the temperature variations to be able to constrain the models of stellar evolution and of galactic chemical evolution. If the temperature fluctuations due to chemical inhomogeneities dominate, the proper abundances to use are those derived from forbidden lines under the assumption that $`t^20.01`$; alternatively if the temperature fluctuations due to the combination of the other six mechanisms mentioned above dominate, the proper abundances to use are those provided by the recombination lines.

The C/H values derived from recombination lines in PNe are in better agreement with chemical evolution models of the Galaxy than the values derived from collisionally excited lines under the assumption that $`t^2=0.000`$ (Carigi 2002). Alternatively the abundances derived from the infrared lines seem to be in favor of chemical inhomogeneities, but these abundances could be higher if large density fluctuations are present.

The presence of abundance gradients in the Galaxy and the values of $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$ and $`\mathrm{\Delta }Y/\mathrm{\Delta }O`$ derived from galactic and extragalactic H II regions indicate that there are significant temperature variations in H II regions and that they are not due to chemical inhomogeneities (see Peimbert & Peimbert 2002 for a review).

#### Acknowledgments.

We are grateful to G. García-Segura, A. Peimbert, and M. Peña, for several excellent suggestions, and to our Australian colleagues for an outstanding meeting.

## References

Bohigas, J. 2002, these proceedings

Carigi, L. 2002, in Ionized Gaseous Nebulae, RevMexAACS, 12, 234

Chu, Y.-H., Guerrero, M. A., Gruendl, R. A., Williams, R. M., & Kaler, J. B. 2001, ApJ, 553, L69

Costero, R., Tapia, M., Echevarría, J., Roth, M., Quintero, A., & Barral, J. F. 1986, RMexAA, 13, 149

Esteban, C. 2002, in Ionized Gaseous Nebulae, RevMexAACS, 12, 56

García-Segura, G., López, J. A., & Franco, J. 2001, ApJ, 560, 928

Garnett, D. R., & Dinerstein, H. L. 2001, RevMexAASC, 10, 13

Gruenwald, R., & Viegas, S.M. 1995, A&A, 303, 535

Hu, J. Y., Slijkhuis, S., Rieu, N.-Q., & de Jong, T. 1993, A&A, 273, 185

Kastner, J. H., Soker, N., Vrtilek, S. D., & Dgani, R. 2000, ApJ, 545, L57

Kastner, J. H., Vrtilek, S. D., & Soker, N. 2001, ApJ, 550, L189

Kato, T., Nogami, D., & Baba, H. 2001, PASJ, 53, 901

Kingdon, J. B., & Ferland, G. J. 1998, ApJ, 506, 323

Kostyakova, E. B. 2002, in Ionized Gaseous Nebulae, RevMexAACS, 12, 167

Liu, X.-W. 2002a, in Ionized Gaseous Nebulae, RevMexAACS, 12, 70

Liu, X.-W. 2002b, these proceedings

Liu, X.-W., & Danziger, I. J. 1993, MNRAS, 263, 256

Liu, X.-W., Luo, S.-G., Barlow, M. J., Danziger, I. J., & Storey, P. J. 2001, MNRAS, 327, 141

Liu, X.-W., Storey, P. J., Barlow, M. J., & Clegg, R. E. S. 1995, MNRAS, 272, 369

Liu, X.-W., Storey, P. J., Barlow, M. J., Danziger, I. J., Cohen, M., & Bryce, M. 2000, MNRAS, 312, 585

Luo, S.-G., & Liu, X.-W 2002, these proceedings

Mathis, J. S. 1976, ApJ, 207, 442

Mathis, J. S., Torres-Peimbert, S., & Peimbert, M. 1998, ApJ, 495, 328

Mellema, G. 1997, in IAU Symposium 180, eds. H. J. Habing & H. J. G. L. M. Lamers, p. 262

Mihalszky, J. S., & Ferland, J. 1983, PASP, 95, 284

Miranda, L. F., Gómez, Y., Anglada, G., & Torrelles, J. M. 2001, Nature, 414, 284

Peimbert, A. 2002, in preparation

Peimbert, A., Peimbert, M., & Luridiana, V. 2002, ApJ, 565, 668

Peimbert, M. 1967, ApJ, 150, 825

Peimbert, M. 1971, Bol. Obs. Tonantzintla y Tacubaya, 6, 29

Peimbert, M. 1995, in The Analysis of Emission Lines, eds. R. E. Williams, & M. Livio (Cambridge), 165

Peimbert, M. 2002, in Ionized Gaseous Nebulae, RevMexAACS, 12, 275

Peimbert, M., & Peimbert, A. 2002, RevMexAASC, in press

Peimbert, M., Peimbert, A., & Ruiz, M. T. 2000, ApJ, 541, 688

Peimbert, M., Sarmiento , A., & Fierro, J. 1991, PASP, 103, 815

Peimbert, M., Storey , P. J., & Torres-Peimbert, S. 1993, ApJ, 414, 626

Peimbert, M., Torres-Peimbert, S., & Luridiana, V. 1995, RevMexAA, 31, 131

Peña, M., Hamann, W.-R., Koesterke, L., Maza, J., Méndez, R. H., Peimbert, M., Ruiz, M. T., & Torres-Peimbert, S. 1997, ApJ, 491, 233

Peña, M., Hamann, W.-R., & Ruiz, M. T. 2002, these proceedings

Peña, M., Torres-Peimbert, S., & Ruiz, M. T. 1991, PASP, 103, 865

Péquignot, D., Amara, M., Liu, X.-W., Barlow, M. J., Storey, P. J., Morisset, C., Torres-Peimbert, S., & Peimbert, M. 2002a, in Ionized Gaseous Nebulae, RMexAASC, 12, 142

Péquignot, D., Liu, X.-W., Barlow, M. J., Storey P. J., & Morisset, C. 2002b, these proceedings

Rola, C., & Pelat, D. 1994, A&A, 287, 676

Rola, C., & Stasińska, G. 1994, A&A, 282, 199

Spitzer L. Jr. 1978, in Physical Processes in the Interstellar Medium (New York: Wiley-Interscience), p. 139

Stasińska, G. 2002, in Ionized Gaseous Nebulae, RevMexAACS, 12, 62

Stasińska, G., & Szczerba, R. 2001, A&A, 379, 1024

Torres-Peimbert, S., Dufour, R. J., Peimbert, M., & Peña, M. 2002, in preparation

Torres-Peimbert, S., Peimbert, M., & Peña, M.: 1990, A&A, 233, 540

Tsamis, Y. G., Barlow, M. J., Liu, X.-W., & Danziger, I. J. 2002, these proceedings

Tylenda, R. 2002, these proceedings

Viegas, S. M., & Clegg, R. E. S. 1994, MNRAS, 271, 993

Zhang, Y., & Liu, X.-W. 2002, these proceedings

Zijlstra, A. A., Te Lintel Hekkert, P., Pottasch, S. R., Caswell, J. L., Ratag, M., & Habing, H. J. 1989 A&A, 217, 157


format: text

text:
# 𝑎𝑏-plane tunneling and Andreev spectroscopy of superconducting gap and pseudogap in (Bi,Pb)2Sr2Ca2Cu3O10+δ and Bi2Sr2CaCu2O8+δ

## 1. Introduction

Along with the usual coherence gap $`\mathrm{\Delta }_s`$, in the spectrum of quasiparticle excitations in high-Tc superconductors there appears a gap $`\mathrm{\Delta }_p`$ (pseudogap), which persists above the superconducting transition temperature $`T_c`$ . Pseudogap has the same $`d`$-symmetry as $`\mathrm{\Delta }_s`$, but disappears (more accurately: becomes indistinct) at some temperature $`T^{}>T_c`$ . The relationship between the pseudogap and superconductivity is far from clear . One of the reasons appears to be that the most popular methods of investigating the excitation spectrum in cuprates, like tunneling and angle-resolved photoemission (ARPES), cannot distinguish between $`\mathrm{\Delta }_p`$ and $`\mathrm{\Delta }_s`$ without recourse to various theoretical models. However, it is known that in the process of Andreev reflection of an electron from the normal metal-superconductor (N-S) interface, a Cooper pair is created in the superconductor . This occurs only in the presence of nonzero energy gap $`\mathrm{\Delta }_s`$ in the superconductor. In other words, the process of Andreev transformation of an electron-hole pair into a Cooper pair is possible only for a reflection from the superconducting order parameter $`\mathrm{\Delta }_s`$. In marked contrast, the tunneling effect is sensitive to any singularity in the quasiparticle excitation spectrum . Therefore, the tunneling characteristics at $`T<T_c`$ in general depend on joint contributions of the energy gap and pseudogap.

$`d`$-wave symmetry of the energy gap introduces some additional complications. Dominant contribution to the junction conductivity in classical (Giaever) tunneling comes from electrons with wave vectors forming a narrow, only a few degrees wide cone . Accordingly, tunnel junctions yield information on the gap anisotropy $`\mathrm{\Delta }(𝐤)`$, and the gap revealed in tunneling experiments can be expressed as $`\mathrm{\Delta }(𝐤)=[\mathrm{\Delta }_s^2(𝐤)+\mathrm{\Delta }_p^2(𝐤)]^{1/2}`$ . In the case of Andreev reflection from a clean N-S interface, the situation is different. The incident electron is not scattered, but reflected back along the same trajectory. This is true for any angle of incidence. It can be said that all incident electrons participate in Andreev reflection on equal rights. Therefore, measurement of a single Andreev N-S junction in principle is sufficient to determine the maximal value of the superconducting gap $`\mathrm{\Delta }_s(𝐤)`$.

In this paper we employ the above discussed characteristic features of tunneling and Andreev spectroscopy to investigate the temperature dependence of the energy gaps $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_s`$ in $`\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_2\mathrm{O}_{8+\delta }`$ (Bi$`2212`$) and (Bi,Pb)2Sr2Ca2Cu3O10+δ \[(BiPb)$`2223`$\] cuprates. The well studied Bi$`2212`$ has two CuO2 layers per unit cell and strong incommensurate modulation in the BiO layer , which complicates the interpretation of tunneling and ARPES data. The substitution of Bi by Pb in the (BiPb)$`2223`$ compound completely erases the superstructure in the BiO layers.

Tunneling measurements were carried out on “break junctions” with the barrier surface practically normal to the crystallographic axes in the base $`ab`$-plane of the material. In the $`c`$ direction, the influence of BiO layer on the tunneling spectra is much more pronounced. Andreev experiments were performed on S-N-S junctions. In both cases we retained only the samples showing pure tunneling or Andreev characteristics. Temperature dependencies of the energy gaps obtained in the two types of experiments were completely different, and testify to fundamental differences between the ”superconducting” gap $`\mathrm{\Delta }_s`$ and the $`a,b`$-axis quasiparticle gap $`\mathrm{\Delta }`$.

## 2. Sample preparation

The tunnel junctions were elaborated from Bi$`2223`$ and Bi$`2212`$ single crystals. Textured $`(\mathrm{Bi}_{1.6}\mathrm{Pb}_{0.4})\mathrm{Sr}_2\mathrm{Ca}_2\mathrm{Cu}_2\mathrm{O}_{10+\delta }`$ and $`\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_2\mathrm{O}_8`$ samples in the form of $`10\times 1\times 0.1\mathrm{mm}^3`$ rectangular bars were prepared  by compacting powdered (BiPb)$`2223`$ and Bi$`2212`$ compounds, respectively at $`3040\mathrm{kbar}`$ between two steel anvils. The powder was contained between two thin copper wires, whose deformation provided uniform pressure distribution in sample volume. In this manner the powder was compacted into dense plane-parallel bars about $`0.1`$mm thick. The bars were then pre-annealed at $`T=845^{}`$C for $`16`$ h, compressed again, and finally annealed at $`T=830^{}`$C for $`14`$ h, obtaining a well pronounced texture. Usually the Bi2212 samples were slightly overdoped and exhibited critical temperature $`T_c=8084`$K. The doping level of the oxygen content was controlled by annealing optimally-doped samples in flowing gas adjusted for different partial pressures of oxygen. The samples emerging from this procedure were highly textured, composed of tightly packed microcrystals aligned in one direction. Sample quality was controlled by transport measurements. We used for further processing only the samples, which were showing critical current density $`\mathrm{J}_\mathrm{c}(\mathrm{T}=4.2\mathrm{K})>410^4\mathrm{A}/\mathrm{cm}^2`$. The superconducting transition temperature $`T_c`$ was determined from the midpoint of the resistive $`R(T)`$ transition (see Fig. 1).

The S-I-S and S-N-S junctions were made by breaking specially prepared Bi$`2212`$ and Bi$`2223`$ samples. Each sample was hermetically sealed by insulating resin and glued to an elastic steel plate, which was then bent until a crack occurred, running across the sample width and detected by monitoring the sample resistance. The hermetic seal remained unbroken in this process. After relieving the external load, the sample returned to its initial position with the crack closed and the microcrystals once again tightly pressed into each other on the line of the fracture. The best alignment is expected in the sample region in which the shear deformation was minimal. This is apparently one of the reasons, why such a procedure results in the realization of one effective junction of the microcrystal-microcrystal type. The selection, from among the competing junctions, of a single junction with minimal tunneling resistance is further assisted by the nature of the tunneling effect, which decreases exponentially with the barrier thickness. Small sample thickness ($`<100\mu `$m) and relatively large size of the microcrystals ($`>10\mu `$m) are also important factors, enhancing junction quality. Such break junctions on microcrystals were found to be particularly effective in the investigation of high-Tc superconductors . The typical normal-state resistance of junctions used in the present study was between a few ohms and a few tens of ohms, and were remarkably stable.

The surface of our Bi$`2212`$ and (BiPb)$`2223`$ break-junctions was perpendicular to the CuO2 plane and the direction of tunneling formed only a very small angle $`\alpha `$ with one of the crystallographic axes ($`a`$ or $`b`$) in this plane, as testified by the presence of the Andreev bound states, seen in the tunneling S-I-S characteristics as a characteristic peak of conductivity at zero bias (cf. Fig. 2). Numerical calculations, based on a simplified theoretical model  and taking into account $`d`$-wave mechanism of pairing, show that the appearance of such narrow zero-bias peak in tunneling conductance occurs at $`\alpha 6^{}`$. In high quality break junctions the zero bias conductance peak (ZBCP) was reported to coexist with the Josephson effect , but we have to rule out this possibility, because of wrong signature: ZBCP was insensitive to magnetic field and did not reflect on the I-V characteristics. The spectra $`\sigma (V)=dI/dV`$ show the quasiparticle peaks at $`2\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is defined as a quarter of the peak-to-peak separation (Fig. 2). We use this $`\mathrm{\Delta }`$ value as a measure of the gap, since there is no exact method of extracting the energy gap from the tunneling spectra, given that the exact functional form of the density of states for high-$`T_c`$ superconductors is not known.

In general, the type of the junction was determined ex post facto from their conductance $`\sigma (V)`$ spectra. We retained for further investigation only the junctions conforming to either S-I-S or S-N-S types. For example, the $`\sigma (V)`$ curve in Fig. 2 reveals all the characteristic features of a superconducting tunnel S-I-S junction: an almost flat region around zero bias followed by a sharp increase in the tunneling current, peaking around $`\pm 60`$meV $`(2\mathrm{\Delta })`$; at still higher bias voltages $`V`$ the conductance depends parabolically on $`V`$. The junction shown in Fig. 3, on the other hand, behaves as a typical Andreev S-N-S junction. First, there is a low resistance region at low bias voltages, seen as a broad pedestal spanning the coordinate origin. The next indication is the excess current, which was observed in all S-N-S junctions included in this study. Finally, the differential conductivity of the junction at e$`V>2\mathrm{\Delta }_s`$ coincides with the normal state conductivity at $`T>T_c`$ \[see inset in Fig. 3(b)\], i.e. for $`T>T_c`$ practically all bias voltage is applied directly to the junction.

One may enquire about the mechanism, which might produce in an apparently random manner either S-I-S or S-N-S junctions. The insulating layer in S-I-S junctions is most probably caused by oxygen depletion. As to the normal barrier, we speculate that the CuO2 planes are more hard to fracture than the buffer layers. After breaking the sample, they penetrate slightly into the buffer layers (see left inset in Fig. 4). In this manner, the coupling between the CuO2 planes belonging to the separated sample parts would be stronger than the normal coupling across the buffer layers, and it could assist to create a constriction, which would act as a normal 3D metal. This hypothesis is in agreement with the scanning microscope study of Bi-2212 single crystal break-junctions of the fracture surfaces, which revealed rough, but stratified, fracture surfaces .

## 3. Experimental Results

Temperature dependence of the energy gap $`\mathrm{\Delta }_s(T)`$ obtained from Andreev S-N-S measurements for Bi$`2212`$ exhibited a BCS-like form (see Fig. 3). We have used two methods to determine $`\mathrm{\Delta }_s`$ for Andreev junctions. The first one is shown in Fig. 4 and relies on measuring the distance between the points of maximal slope changes of $`\sigma (V)`$ plot is taken as the measure of $`4\mathrm{\Delta }_s`$. The details of the second one are shown in top inset in Fig. 4. The rationale for both methods is in recent calculations , based on the Klapwijk, Blonder and Tinkham  treatment of multiple Andreev reflections between two superconductors, which indicate that $`2\mathrm{\Delta }_s`$ is determined by the separation of extrema in $`d\sigma /dV`$. The results obtained by both methods are plotted together in Fig. 3(b). As seen, these results differ slightly, but both outline essentially the same $`\mathrm{\Delta }_s(T)`$ dependence.

The $`\mathrm{\Delta }(T)`$ gap dependence, determined from tunneling measurements performed on the same compounds, diverged considerably from the BCS relation (Fig. 5). In fact, $`\mathrm{\Delta }(T)`$ depends on temperature very weakly for $`TT_c`$. According to ARPES investigation , such behavior of $`\mathrm{\Delta }(T)`$ in Bi2212 near optimal doping is expected for the $`a`$ (or $`b`$) direction in the CuO2 plane. This result agrees with our assumption about the direction of tunneling in our Bi$`2212`$-I-Bi$`2212`$ junctions. As mentioned above, a further confirmation is provided by the presence of Andreev bound state, seen in the spectra of S-N-S and S-I-S junctions as a characteristic peak of conductivity at zero bias (cf. Fig. 3 and Fig. 5). According to the ARPES data , near optimal doping the $`\mathrm{\Delta }(T)`$ gap becomes temperature dependent only when $`\alpha `$ is of the order of $`15^{}`$. For technological reasons, the formation of break junctions with crystal broken at such an angle is not probable. As a result, the tunneling characteristics at $`T>T_c`$ relate to the gap in $`(100)`$ or $`(010)`$ direction.

In full agreement with the ARPES results , with increasing temperature the gap $`\mathrm{\Delta }`$ of Bi$`2212`$ becomes filled with quasiparticle excitations, and the conductance peaks at $`2\mathrm{\Delta }`$ become less distinct. The distance between the still discernible conductance peaks does not decrease, and the $`\mathrm{\Delta }(T)`$ gap is seen to continue into the region $`T>T_c`$. Similar behavior is observed also for the (BiPb)$`2223`$ compound.

The temperature dependence of the proper coherent gap $`\mathrm{\Delta }_s(T)`$ behaves in a completely different manner (Fig. 3). With increasing temperature, the gap narrows, and at $`T=T_c`$ it closes completely. High curvature of the Andreev conductance dip at $`eV2\mathrm{\Delta }_s`$ is evidence both of the good quality of the investigated junctions and of the long lifetime of quasiparticles in the gap region. This was confirmed by the analysis of spectra of the normal metal - constriction - superconductor (N-c-S) junctions . For Bi$`2212`$ Andreev N-c-S junction, the Blonder-Tinkham-Klapwijk  parameter $`Z`$ used to obtain theoretical fit was small, $`Z0.5`$, a value characteristic for very clean N-S contacts.

For energies beyond the gap $`\mathrm{\Delta }`$ value, tunneling in the $`ab`$-plane of (BiPb)$`2223`$ S-I-S junction revealed the so-called dip and hump structures, as shown in Fig. 2 and Fig. 6. In Fig. 6, the voltage axis is normalized to the voltage $`eV_p=\mathrm{\Delta }`$, and the conductance axis is normalized to the background; the spectra are shifted vertically for clarity. The dip and hump features roughly scale with the gap $`\mathrm{\Delta }`$ for various oxygen doping levels (see Fig. 7). There is, however, a slight deviation of the data from the straight $`\mathrm{\Delta }`$ line.

## 4. Theoretical implications

The considerable interest in the pseudogap investigation is stimulated to a great extent by the theoretical models of high-Tc superconductivity, in which pseudogap appears as a precursor of the superconducting gap , e.g. bipolaron model . In another group of models, the appearance of pseudogap is related to some sort of magnetic pairing . However, the domains of applicability of these models are not very strictly defined and it is quite possible that the pseudogap (like high temperature superconductivity) is caused by several simultaneously acting mechanisms.

For example, in the Emerson-Kilverson-Zachar (EKZ) theory  the crucial role in the formation of high temperature superconductivity is ascribed to the separation of spin and charge, arising as a result of partitioning of CuO2 planes into narrow conducting and dielectric stripes. “Pairing” at $`T^{}>T_c`$ in the EKZ model means the formation of a spin gap. A wide spin gap (or pseudogap $`\mathrm{\Delta }_p`$) is indeed formed in the space limited, hole-free region, such as the region between the conducting stripes. Phase-coherent (i.e. really superconducting) state is created only at $`T<T_c`$. The model well explains the smooth transition of pseudogap into tunneling gap $`\mathrm{\Delta }`$ when the temperature decreases below $`T_c`$. However, the observed dependence on temperature of the order parameter gap $`\mathrm{\Delta }(T)`$ at $`T<T_c`$ is fundamentally different from that of the gap $`\mathrm{\Delta }_s(T)`$, as shown in Fig. 3 and Fig. 5. It is not clear how in the phase fluctuation picture appears the BCS-like $`\mathrm{\Delta }_s(T)`$ dependence. Such a situation would be possible e.g. in the generation of charge (and spin) density waves with the superconducting gap and pseudogap competing for the same region of the Brillouin zone . Then the transition to the superconducting state could occur in the presence of a pseudogap in normal excitations, opening e.g. in the electron-hole channel (i.e. a pseudogap, which would not transform directly into the superconducting gap, as in the Emery-Kivelson model).

There are numerous experiments, which confirm the essentially different nature of the superconducting gap $`\mathrm{\Delta }_s`$ and gap (pseudogap) $`\mathrm{\Delta }`$ . The most convincing are intrinsic $`c`$-axis tunneling experiments (in stacked layers) . However, they yield different results from the point contact, scanning tunneling spectroscopy (STM), and break junction experiments: the hump was observed at the energy of $`2\mathrm{\Delta }`$, instead of $`4\mathrm{\Delta }`$. The authors note a similarity between the observed $`c`$-axis pseudogap and Coulomb pseudogap for tunneling into a two-dimensional electron system. In our case, the tunneling and Andreev reflection were realized in $`ab`$-plane, and together with the $`\mathrm{\Delta }(T)`$ dependence (Fig. 4), we clearly observed the peak-dip-hump structure (Figs. 2 and 3). Position of dip and hump for SIS junctions was at $`3\mathrm{\Delta }`$ and $`4\mathrm{\Delta }`$ (Fig. 2). This suggests that the observed dip-hump structure may originate from short-range magnetic correlations in the $`ab`$-plane . Then the gap $`\mathrm{\Delta }`$ would be the fermionic excitation gap and $`\mathrm{\Delta }_s`$ – the mean-field order parameter. It should be emphasized, finally, that the observed $`\mathrm{\Delta }_s(T)`$ dependence exhibits non-BCS behavior at $`T0`$ (Fig. 3).

In summary, our $`ab`$-plane tunneling and Andreev spectroscopy studies of normal and slightly overdoped (BiPb)$`2223`$ and Bi$`2212`$ compounds show presence both of a superconducting energy gap $`\mathrm{\Delta }_s`$, corresponding to the $`d`$-wave Cooper pairing, and a dip-hump structure at $`3\mathrm{\Delta }`$ and $`4\mathrm{\Delta }`$ (for the SIS junction). This suggests that high-energy pseudogap, which is associated with the dip and hump, could be magnetic in origin. The gap $`\mathrm{\Delta }`$ is nearly temperature independent and becomes blurred above $`T_c`$, being continuously transformed with increasing temperature into the pseudogap. In contrast, the order parameter gap $`\mathrm{\Delta }_s(T)`$ has a strong temperature dependence and for $`T0`$ reveals a non-BCS mean field behavior. Our findings are in general agreement with those of Deutscher , although it must be emphasized again that we have considered the slightly overdoped case.

## Acknowledgments

This work was supported by Polish Government (KBN) Grant No PBZ-KBN-013/T08/19.


text:
# Anomalous Heat Conduction in a Di-atomic One-Dimensional Ideal Gas

## Abstract

We provide firm convincing evidence that the energy transport in a one- dimensional gas of elastically colliding free particles of unequal masses is anomalous, i.e. the Fourier Law does not hold. Our conclusions are based on the analysis of the dependence of the heat current on the number of particles, of the internal temperature profile and on the Green-Kubo formalism.

PACS number: 05.45.-a

Understanding the dynamical origin for the validity of the Fourier law of heat conduction in deterministic one-dimensional particle chains is one of the oldest and most frustrating problems in non-equilibrium statistical physics. Due to some very basic unresolved issues the problem has been a source of many recent publications.

In the absence of analytical results, these papers are mainly oriented towards a numerical analysis of the problem. However, due to the delicate nature of the questions under discussion, numerical results sometimes lead to different conclusions. This is the case, for example, of the 1d hard-point particles with alternating masses for which opposite conclusions have been reached . This disagreement is not extremely surprising since this system lies in the foggy region which separates clear, regular integrable systems, from the totally chaotic, deterministic, motion. Indeed this system has zero Lyapounov exponent and therefore it lacks of the exponential local instability which characterizes chaotic systems. On the other hand it has been shown  that such systems can exhibit Gaussian diffusive behaviour and, more recently , an example has been shown of a system with zero Lyapounov exponent which however obeys the Fourier law. From the point of view of a general theoretical understanding, the fact that the alternating mass problem lies in this critical region renders particularly important to establish, beyond any reasonable doubt, its conducting properties. This is what we set up to do in the present paper.

We consider a one-dimensional gas of interacting particles with coordinates $`q_n`$ and momenta $`p_n`$ for which the hamiltonian can be written in the form

$$H=\underset{n=0}{\overset{N1}{}}h_n,h_n=\frac{p_n^2}{2m_n}+V(q_{n+1}q_n).$$
(1)
The energy current from site $`n`$ to site $`n+1`$ is defined as $`j_n=\{h_{n+1},h_n\}`$ and satisfies the continuity equation $`(d/dt)h_n=\{H,h_n\}=j_nj_{n1}`$. In particular we focus our attention on the ideal gas model of elastically colliding particles, $`V(q>0)=0,V(q<0)=\mathrm{}`$, with alternating masses, $`m_{2n1}=m_1=\sqrt{r}`$, $`m_{2n}=m_2=1/\sqrt{r}`$, where the ratio $`r=m_1/m_2`$ serves as a model parameter. We have mainly considered the value $`r=(\sqrt{5}1)/2`$; however all the reported numerical results have been checked also for several other values of $`r`$ ($`0<r<1`$) where we found no qualitative distinction.

We place our system of $`N`$ particles between two stochastic maxwellian heat reservoirs at temperatures $`T_\mathrm{L}`$ and $`T_\mathrm{R}`$ (see  for a description of the reservoir model). We chose the temperatures of the reservoirs, $`T_\mathrm{L}=1`$ and $`T_\mathrm{R}=2`$, and measure the long-time averaged heat current $`J=lim_t\mathrm{}(1/t)_0^t𝑑t^{}J(t^{})`$ versus the system size $`N`$, where $`J=(1/N)_{n=1}^{N1}j_n`$. Here we want to note that strict equivalence between our definition of the heat current (which simply accounts for the energy transferred during collisions), and the ‘free particle’ current $`j_n=m_nv_n^3/2`$ used by some authors, e.g  (which does not ’feel’ the collisions), is nontrivial in general .

In order to ensure that our results are not affected by finite size effects we have put particular care in using an efficient numerical scheme which allows to reach high N values. Our algorithm, developed for the first time in Ref., searches in a partially ordered tree (heap) of pre-computed candidates pairs for the next collision and, due to this, its requires only $`\mathrm{log}_2N`$ computer operations per collision. As a consequence, we were able to simulate very long chains and we have obtained reliably converged results for lattices with sizes $`N`$ up to 30000. Convergence has been controlled by checking the constancy of the finite-time-averaged heat current $`(1/t)_0^t𝑑t^{}J(t^{})`$, and to this end simulations for the largest system sizes had to be carried on up to $`10^{12}`$ pair collisions. It is also clear that convergence problems suggest to keep far away form the $`r`$ values too close to one or to zero. The analysis made in  indicated that the range $`0.15<r<0.6`$ was the most effective in attenuating solitary pulses and the value $`r=0.2`$ was chosen. In the present paper we take the somehow ’standard’ choice $`r=(\sqrt{5}1)/2`$.

Validity of Fourier law implies the scaling $`JTN^1`$. Our numerical results shown in fig. 1 clearly demonstrate instead a different power-law behaviour namely $`JN^\alpha `$ with $`\alpha 0.745\pm 0.005`$ over a very large range in $`N`$. We have also found that the scaling exponent $`\alpha `$ does not change appreciably with the mass ratio $`r`$. For example for ten times smaller value of $`r`$ the asymptotic scaling only sets in later (i.e. for larger values of $`N`$, see fig. 1). The possibility of a slow convergence to the asymptotic value might be at the origin of the slightly different numerical values for $`\alpha `$ found in previous numerical experiments ($`\alpha 0.65`$ in Hatano , $`\alpha 0.83`$ in Dhar ). Since the model under consideration is energy scaling we do not expect any dependence of the exponent $`\alpha `$ on the reservoirs temperatures.

The above results therefore solve the existing controversy and clearly show that the alternating mass, 1d hard point gas does not obey Fourier heat law. We turn now to the analysis of other quantities which, besides providing additional confirmation of the above conclusions, illuminate interesting aspects of the heat conduction problem. A quantity of main interest is the internal local temperature profile $`T_n=p_n^2/(2m_n)`$ in the non-equilibrium steady state for the system placed in between the heat reservoirs. First we notice that the temperature profile in discrete index variable $`n`$ is different than the temperature profile in position variable $`q_n`$  since the inverse density $`dq/dn=q_{n+1}q_n`$ is non-uniform in non-equilibrium, in fact it is simply proportional to the temperature due to constancy of pressure . Now, in case of a Fourier law, the thermal conductivity $`\kappa `$ scales with temperature like $`\kappa \sqrt{T}`$. Therefore, from $`\sqrt{T}(dT/dn)dn/dq=\mathrm{const}`$ one obtains the temperature profile $`T_n^{\mathrm{kin}}=(T_\mathrm{L}^{1/2}+(T_\mathrm{R}^{1/2}T_\mathrm{L}^{1/2})n/N)^2`$. Extensive numerical simulations showed that the temperature profile in our model converges, for sufficiently large $`N`$, to a well-defined scaling function $`T_n^{\mathrm{scal}}=f(n/N)`$ which is slightly, but significantly different, from the kinetic temperature profile $`T_n^{\mathrm{kin}}`$. This is another evidence for the anomalous heat transport and for the non-validity of the Fourier law in our system. It should be remarked that the convergence to the temperature profile predicted by kinetic theory observed in , which has indeed been considered as surprising by the author himself, actually does not take place.

A standard theoretical analysis of transport laws is based on Kubo formulae . However, applicability of Kubo formula in momentum conserving cases, i.e. for translationally invariant systems like model (1), is not very clear. This is particularly critical in view of a recent claim  that Kubo formula diverges for a momentum conserving lattice with non-vanishing pressure. For this latter type of models we have an additional difficulty in applying the Kubo formalism, namely, as we show below, the result depends not only on the temperature gradient, but also on other thermodynamic properties of the initial non-equilibrium state – i.e. the isobaric case (constant pressure profile) or the isodense case (constant density profile). There is no a-priori argument which favours either of these two options: the choice depends on the specific physical situation of interest. For instance the steady-state heat current simulation considered above (figs. 1,2) is clearly described by the isobaric state. Since we want to consider both situations we need to revise the derivation of Kubo formula by following the time evolution of a general non-equilibrium initial state in an isolated system with periodic boundary conditions $`q_Nq_0+N,p_Np_0`$. To this end, we prepare the initial state in a local-equilibrium state described by the inverse temperature profile $`\beta _n`$ and by an additional thermodynamic potential $`\gamma _n`$

$$\rho _{\mathrm{neq}}=Z_{\mathrm{neq}}^1\mathrm{exp}\left(\underset{n}{}\beta _nh_n\underset{n}{}\gamma _n(q_{n+1}q_n)\right)$$
(2)
This (small) additional term is necessary in order to equilibrate the pressure in the isobaric case. Notice that $`\gamma _n`$ is undetermined up to an arbitrary additive constant due to a gauge invariance of the second term in (2). In order to determine the gradient of the $`\gamma `$-potential which is necessary to keep the physical pressure $`\varphi `$ constant ($`n`$-independent), we compute the generalized pressure $`\varphi _l`$

$`\beta _l\varphi _l`$ $`=`$ $`{\displaystyle \frac{}{a}}\mathrm{ln}Z_l(a)|_{a=0},`$ (3)
$`Z_l(a)`$ $`=`$ $`{\displaystyle e^{_n(\beta _nV(q_{n+1}q_n+a\delta _{ln})\gamma _n(q_{n+1}q_n+a\delta _{ln}))}𝑑\stackrel{}{q}}`$ (4)
By a simple trick, a shift of one variable $`q_lq_l+a`$ in the integral $`Z_l(a)`$, we find $`Z_l(a)Z_{l1}(a)`$ and therefore

$$\beta _l\varphi _l=\beta _{l1}\varphi _{l1}=\mathrm{const}.$$
(5)
Writing the physical pressure (force) as $`\varphi =V^{}(q_{n+1}q_n)_{\mathrm{neq}}=\varphi _n+\gamma _n/\beta _n`$, multiplying by $`\beta _n`$, and taking the first difference $`f_n:=f_nf_{n1}`$ we obtain the required ’gradient’

$$\gamma _n=\varphi \beta _n.$$
(6)
In the following we consider two specific cases: (i) The initial isodense state with $`q_{n+1}q_n_{\mathrm{neq}}=\mathrm{const}`$. This is obtained by putting $`\gamma _n0`$. (ii) The initial isobaric state with uniform pressure profile. This is obtained by specifying the $`\gamma `$-potential according to eq. (6). We note again that the isobaric state (ii) is the one which is relevant for the steady non-equilibrium state of a system in contact with heat reservoirs. Carefully repeating the first few steps in the derivation of the Green-Kubo formula (following Ref.) we arrive at the very general linear response formula

$$A(t)A(t_0)_{\mathrm{neq}}=_{t_0}^t𝑑t^{}A(t^{})\underset{n}{}(\beta _nj_n+\gamma _nv_n)_{\mathrm{eq}}$$

where $`v_n=\dot{q}_n`$ are the particles’ velocities. In the last step we have assumed that we are close to equilibrium ($`\beta _n`$ and $`\gamma _n`$ small) so that the RHS can be evaluated in the corresponding equilibrium state $`_{\mathrm{eq}}`$. Let us now consider the periodic temperature profile $`\beta _n=\beta +ϵ\mathrm{sin}(2\pi kn/N)`$, and compute the total heat that has been transported in time $`t`$ from warm to cold regions of the lattice, namely $`Q(t)=_0^t𝑑t^{}J_k(t)`$ where $`J_k:=N^1_{n=0}^{N1}\mathrm{cos}(2\pi kn/N)j_n`$. Inserting $`Q`$ for $`A`$ and using (6) (case (i) is obtained by formally setting $`\varphi =0`$) we obtain

$$Q(t)_{\mathrm{neq}}=\frac{2\pi ϵ}{N}_0^t𝑑t^{}(tt^{})J_k(t)(J_k+\varphi V_k)_{\mathrm{eq}}.$$
(7)
where $`V_k:=N^1_{n=0}^{N1}\mathrm{cos}(2\pi kn/N)v_n`$. We see that the growth of the transported heat $`Q(t)_{\mathrm{neq}}`$ is given by the generalized correlation function $`C_k(t)=C_{JJ}(t)+\varphi C_{JV}(t)`$, where $`C_{JJ}(t)=J_k(t)J_k_{\mathrm{eq}}`$ and $`C_{JV}(t)=J_k(t)V_k_{\mathrm{eq}}`$. In the isodense case (i) expression (7) reduces to the usual current autocorrelation function only.

In the case of Fourier law, we expect initial linear growth $`Q(t)_{\mathrm{neq}}ϵ\kappa t`$, while for the ballistic heat transport we expect quadratic growth $`Q(t)_{\mathrm{neq}}t^2`$ (this has been confirmed numerically for the integrable gas of equal masses $`r=1`$). However, in a generic system with momentum conservation and non-vanishing pressure $`\varphi 0`$, like our diatomic hard-point gas, we find qualitatively different behavior in cases (i) and (ii). For example, due to the result , $`C_{JJ}(t)`$ has a plateau $`\varphi ^2`$ and the transport is ballistic in the isodense case, while in the isobaric case the second term, $`C_{JV}(t)`$, compensates for the plateau and yields a much slower increase of the transported heat. In this latter case, independent numerical computations of $`Q(t)_{\mathrm{neq}}`$ and of $`C_k(t)`$ for $`N`$ up to $`32768`$ shown in fig. 3 give $`Q(t)_{\mathrm{neq}}t^\nu `$ with $`\nu 1.255`$, which is still clearly super-diffusive, and directly validate the formula (7).

In fig. 4 we show the generalized correlation function $`C_k(t)`$ for $`k=0`$ and $`k=1`$ separately. Note that the results for $`k=1`$ exhibits some oscillations for longer times due to finite size effects, namely due to periodicity of the lattice, while $`k=0`$ gives the spatially homogeneous correlation function which has the same long-time behavior with weaker finite size effects \[however, the case $`k=0`$ is not strictly related to Kubo formula (7)\]. We see that in the time range where $`C_0(t)`$ and $`C_1(t)`$ match, the asymptotic decay is compatible with $`C_k(t)t^\mu `$ with the exponent $`\mu =2\nu =0.745`$ consistent with eq. (7), and satisfying $`\mu =\alpha `$.

This results can be interpreted in the following way. In the isodense initial state the initial temperature gradients drive the heat ballistically in terms of sound waves . On the other hand, in the isobaric initial state, we have density gradients which drive the heat in the opposite direction and almost exactly compensate for the effect of temperature gradients so that the net effect is a sub-ballistic, but still super-diffusive, energy transport. In order to illustrate the mechanism of ballistic energy transport we show in fig. 5 the spatio-temporal current-current correlation function $`c_{jj}(n,t)=j_0j_n(t)_{\mathrm{eq}}`$ which exhibits clear ballistic tongues along the lines $`n=\pm c_st`$ where $`c_\mathrm{s}=1.78`$.

In this paper we have discussed the thermal conducting properties of a one dimensional hard point gas with alternating masses. This problem has a long history and recently several numerical results have appeared which lead to contrasting conclusions on whether Fourier law is obeyed or not. In the latter case, different values have been found for the rate of divergence of coefficient of thermal conductivity as a function of the particles number. Indeed, the slow convergence properties which sometimes characterize this problem suggest particular care in the interpretation of numerical findings.

Here we have presented accurate numerical analysis, made possible by a powerful integration scheme, which allows us to establish definite convincing evidence that the system does not obey the Fourier law. Moreover, by considering a typical mass ratio $`r=(\sqrt{5}1)/2`$, we have found that the asymptotic scalings: (i) steady-state heat current between heat baths $`JN^\alpha `$, (ii) heat transported within a non-equillibrium isobaric initial state in isolated system $`Q(t)t^{2\alpha }`$, and (iii) generalized current-velocity correlation in equilibrium state $`C_k(t)t^\alpha `$, are described by just one exponent $`\alpha =0.745`$.

This work has been partially supported by EU Contract N0. HPRN-CT-1999- 00163(LOCNET network) and by MURST (Prin 2000, Caos e localizzazione in Meccanica Classica e Quantistica. TP acknowledges financial support by the Ministry of Science, Education and Sport of Slovenia.


text:
# Single-particle density matrix and superfluidity in the two-dimensional Bose Coulomb fluid

## I INTRODUCTION

The fluid of point charges interacting via a logarithmic Coulomb potential in a strictly two-dimensional (2D) space is a basic model in statistical mechanics with main relevance to the theory of vortex fluctuations in superfluid or superconducting films. Magro and Ceperley have discussed the bosonic ground state of this model system. They first used a sum-rule argument from the work of Pitaevskii and Stringari to show that the zero-point fluctuations associated with long-wavelength plasmons rule out the presence of a Bose-Einstein condensate even at zero temperature. They proceeded to evaluate the single-particle density matrix $`\rho (r)`$ by quantum Monte Carlo (QMC) methods and showed that the ground state exhibits algebraic off-diagonal long-range order: the decay of $`\rho (r)`$ with increasing distance $`r`$ is through the power law $`r^{r_s/4}`$, where $`r_s`$ is a coupling strength parameter determined by the areal density $`n`$.

As Magro and Ceperley emphasize, although the ground state of the 2D Bose Coulomb fluid (2D-BCF) is not condensed, superfluidity is nonetheless possible. If the model is indeed superfluid, it may then provide an ideal system in which to study the differences between superfluidity and Bose-Einstein condensation.

A number of structural and collective dynamical properties have been calculated for the 2D-BCF within approximate theoretical schemes, but no other studies seem to have been made of its momentum distribution and one-body density matrix beyond that of Magro and Ceperley. In the present paper we evaluate these properties for the charged fluid in the weak-coupling regime corresponding to high density, by adopting a hydrodynamic reduction of its Hamiltonian as previously proposed by Popov for 2D neutral Bose fluids (see also Meng). A crucial point of our approach is to allow for a difference between the long-wavelength dispersion relations of single-particle and collective excitations at finite temperature $`T`$, through the use of distinct $`f`$-sum rules on the superfluid and the total particle number densities. We recover by this approach the results of Magro and Ceperley for the power-law decay of $`\rho (r)`$ in the ground state and obtain their extension to finite temperature. We show that within this approach the 2D-BCF is indeed superfluid at $`T=0`$ and that in the weak-coupling regime a slowly declining superfluid fraction persists as the temperature is raised. In fact, the superfluid fraction enters to determine the value of the exponent for the power-law decay of $`\rho (r)`$. The presence of the superfluid density in the power-law decay of the density matrix in a neutral 2D Bose gas was already argued by Popov.

Our present progress in understanding the properties of the 2D-BCF borrows from the theory of a quasicondensate state in the neutral 2D Bose gas, that was developed by a number of authors and has also been used within a Bogoliubov approach to describe the 2D fluid of charged bosons with $`e^2/r`$ interactions. In essence, the local properties of a quasicondensate, over distances that are shorter than the phase correlation length, are the same as those of a genuine condensate. However, in our treatment of the 2D-BCF we dispense with the notion of a quasicondensate fraction and base our arguments on the $`f`$-sum rule for the superfluid fraction.

The paper is organized as follows. In Section II, after introducing the Hamiltonian of the charged fluid with its neutralizing background and recalling the $`f`$-sum rules for a superfluid from the work of Hohenberg and Martin, we present the hydrodynamic Hamiltonian as expressed in terms of density and phase operators. Following Meng we diagonalize it in the momentum representation and obtain a dispersion relation for single-particle excitations which has the Bogoliubov form but contains the superfluid density in place of the condensate density. In Section III we use the hydrodynamic Hamiltonian to determine the power-law decay of the one-body density matrix from the behavior of the correlations between phase fluctuations entering the single-particle Green’s function. In Section IV we evaluate $`\rho (r)`$ by using the resummation method adopted by other authors to deal with a fluid with a quasicondensate, and show that the power-law decay obtained from the single-particle Green’s function is recovered. We then report in Section V numerical results for the superfluid density from Landau’s formula and for the one-body density matrix as functions of temperature at various values of the coupling strength $`r_s`$ in the weak-coupling regime. Section VI concludes the paper with a brief summary. In an Appendix we show that the standard Bogoliubov approach would indeed yield that the condensate fraction is zero in the ground state of the 2D-BCF.

## II HYDRODYNAMIC HAMILTONIAN AND SINGLE-PARTICLE EXCITATIONS

The 2D-BCF is described by the Hamiltonian

$$H=𝑑𝐫\psi ^{}(𝐫)\left[\frac{^2}{2m}\mu \right]\psi (𝐫)+\frac{1}{2}𝑑𝐫𝑑𝐫^{}\psi ^{}(𝐫)\psi ^{}(𝐫^{})V(|𝐫𝐫^{}|)\psi (𝐫^{})\psi (𝐫)$$
(1)
(having set $`\mathrm{}=1`$), where $`\psi (r)`$ is the field operator, $`\mu `$ is the chemical potential and $`V(r)`$ is the interaction potential. This is the solution of the 2D Poisson equation $`^2V(r)=2\pi e^2\delta (𝐫)`$, yielding $`V(r)=e^2\mathrm{ln}(r/l_0)`$ where $`l_0`$ is a reference length that we shall take as $`l_0=(me^2)^{1/2}`$. The Fourier transform of the potential is $`V_𝐤=2\pi e^2/k^2`$, and in the following we shall set $`V_{𝐤=0}=0`$ on account of a uniform background ensuring global charge neutrality. The coupling strength is measured by the dimensionless parameter $`r_s`$, which is defined by $`r_s=(2me^2/\pi n)^{1/2}`$.

A special role will be played in the following by a sum rule involving the superfluid density $`n_s`$. We first recall that the usual $`f`$-sum rule involving the particle density $`n`$ can be recast through the continuity equation into a sum rule on the longitudinal current-current response function $`\chi _{JJ}(𝐤,\omega )`$,

$$_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{\pi }\frac{\mathrm{Im}\chi _{JJ}(𝐤,\omega )}{\omega }=\frac{n}{m}.$$
(2)
For a charged fluid this is equivalent to the well-known plasmon sum rule on the longitudinal electrical current density. In addition to Eq. (2), a sum rule involving the response function of the superfluid velocity $`v_s`$ holds in a superfluid, which reads

$$\underset{𝐤0}{lim}_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{\pi }\frac{\mathrm{Im}\chi _{v_sv_s}(𝐤,\omega )}{\omega }=\frac{1}{mn_s}.$$
(3)
The role of this relation as a long-wavelength $`f`$-sum rule on the superfluid response becomes evident from its similarity to Eq. (2) when it is rewritten as a sum rule for the superfluid particle current density $`𝐉_s=n_s𝐯_s`$.

We are now ready to proceed to a hydrodynamic reduction of the Hamiltonian (1), following the work of Popov and Meng on neutral superfluids. This is obtained by making the transformation

$$\psi (𝐫)=\sqrt{\widehat{\rho }(𝐫)}\mathrm{exp}[i\mathrm{\Phi }(𝐫)],\psi ^{}(𝐫)=\mathrm{exp}[i\mathrm{\Phi }(𝐫)]\sqrt{\widehat{\rho }(𝐫)}$$
(4)
for the components of the field operators which correspond to wave number $`k`$ below a cut-off $`k_o`$. Here $`\widehat{\rho }(𝐫)`$ and $`\mathrm{\Phi }(𝐫)`$ are the particle density and phase operators. The transformed Hamiltonian can be brought to a quadratic form by considering only small fluctuations in $`\widehat{\rho }(𝐫)`$ around a constant value $`\rho _0`$,

$$\widehat{\rho }(𝐫)=\rho _0+\eta (𝐫).$$
(5)
A term linear in $`\eta (𝐫)`$ can be dropped by setting $`\mu =V_{𝐤=0}=0`$, and we obtain the hydrodynamic Hamiltonian

$$H_h=𝑑𝐫\left[(8mn_s)^1\left(\eta (𝐫)\right)^2+(n_s/2m)\left(\mathrm{\Phi }(𝐫)\right)^2\right]+\frac{1}{2}𝑑𝐫𝑑𝐫^{}V(|𝐫𝐫^{}|)\eta (𝐫)\eta (𝐫^{}).$$
(6)
Here, the constant $`\rho _0`$ has been taken equal to $`n_s`$ in order to satisfy the sum rule (3), as we shall explicitly demonstrate at the end of this Section. As usual, the superfluid velocity is $`𝐯_s=\mathrm{\Phi }(𝐫)/m`$. We stress that the Hamiltonian in Eq. (6) describes only long-wavelength fluctuations of small amplitude in the superfluid density and velocity field.

The Hamiltonian (6) can be diagonalized by first expanding the density fluctuation and phase operators in the form

$$\eta (𝐫)=\underset{𝐤(k<k_0)}{}q_𝐤\mathrm{exp}(i𝐤𝐫),\mathrm{\Phi }(𝐫)=\underset{𝐤(k<k_0)}{}p_𝐤\mathrm{exp}(i𝐤𝐫),$$
(7)
which yields

$$H_h=\underset{𝐤(k<k_0)}{}\left\{\left[(k^2/8mn_s)+\frac{1}{2}V_𝐤\right]q_𝐤q_𝐤+(n_sk^2/2m)p_𝐤p_𝐤\right\}.$$
(8)
We next set $`q_𝐤=(2\alpha _𝐤)^{1/2}(a_𝐤+a_𝐤^{})`$ and $`p_𝐤=i\alpha _𝐤(2\alpha _𝐤)^{1/2}(a_𝐤a_𝐤^{})`$ with $`[a_𝐤,a_𝐤^{}^{}]=\delta _{𝐤,𝐤^{}}`$, where $`\alpha _𝐤=2mE_𝐤/(n_sk^2)`$ and

$$E_𝐤=\sqrt{\frac{2\pi n_se^2}{m}+\left(\frac{k^2}{2m}\right)^2}.$$
(9)
This finally yields

$$H_h=\underset{𝐤(k<k_0)}{}E_𝐤a_𝐤^{}a_𝐤.$$
(10)
Equation (9) has the same form as the Bogoliubov spectrum for a weakly coupled charged Bose gas, with the superfluid density $`n_s`$ taking the place of the condensate density. It will be used in the next Section as the dispersion relation of single-particle excitations in the evaluation of the phase-phase correlations from the single-particle Green’s function. Equation (9) should also be contrasted with the dispersion relation of collective excitations, which within a weak-coupling theory is

$$\omega _𝐤=\sqrt{\frac{2\pi ne^2}{m}+\left(\frac{k^2}{2m}\right)^2}.$$
(11)
This relation tends to the 2D plasma frequency $`\omega _p=(2\pi ne^2/m)^{1/2}`$ at long wavelengths, in accord with the sum rule reported in Eq. (2). Notice that in the limit $`T0`$, where we expect that $`n_sn`$, Eqs. (9) and (11) give the same result. According to the Gavoret-Nozières theorem, single-particle and collective excitations coincide at long wavelengths in a Bose-condensed fluid at zero temperature.

We conclude this Section by showing that the sum rule (3) is satisfied in our approach. The response function $`\chi _{v_sv_s}(𝐤,\omega )`$ for the superfluid velocity is related to the phase-phase response by $`\chi _{v_sv_s}(𝐤,\omega )=k^2\chi _{\varphi \varphi }(𝐤,\omega )/m^2`$. At long wavelengths we get from Eq. (8)

$$\underset{k0}{lim}\chi _{v_sv_s}(𝐤,\omega )=\frac{k^2V_𝐤}{2m^2\mathrm{\Omega }_p}\left(\frac{1}{\omega \mathrm{\Omega }_p+i\epsilon }\frac{1}{\omega +\mathrm{\Omega }_p+i\epsilon }\right)$$
(12)
where $`\epsilon =0^+`$ and $`\mathrm{\Omega }_p=(2\pi n_se^2/m)^{1/2}`$. Insertion of Eq. (11) into the integral on the LHS of Eq. (3) leads immediately to the desired result.

## III ASYMPTOTIC BEHAVIOR OF THE SINGLE-PARTICLE GREEN’S FUNCTION

We show in this Section that the one-body density matrix of the 2D-BCF at temperature $`T`$ below a critical temperature $`T_c`$ has a power-law decay, as a consequence of the correlations between phase fluctuations in the superfluid. We follow the method proposed for the neutral 2D gas in the work of Popov (see also Fisher and Hohenberg), which derives the power law by evaluating the single-particle Green’s function at low momenta from the mean-square fluctuations of the phase $`\mathrm{\Phi }(𝐫,\tau )`$.

More precisely, the phase fluctuations determine the single-particle Green’s function $`G(𝐫,\tau ;𝐫_1,\tau _1)`$ in the low-momentum regime (below the cut-off momentum $`k_0`$) according to

$$G(𝐫,\tau ;𝐫_1,\tau _1)\mathrm{exp}\{\frac{1}{2}[\mathrm{\Phi }(𝐫,\tau )\mathrm{\Phi }(𝐫_1,\tau _1)]^2\}.$$
(13)
From Eq. (19.16) in Chapter 6 of Popov’s book and using the dispersion relation for single-particle excitations given in Eq. (9) we find

$$\mathrm{\Phi }(𝐤,\omega )\mathrm{\Phi }(𝐤,\omega )\frac{V_𝐤}{\omega ^2+n_sV_𝐤k^2/m}$$
(14)
so that

$`{\displaystyle \frac{1}{2}}[\mathrm{\Phi }(𝐫,\tau )\mathrm{\Phi }(𝐫_1,\tau _1)]^2`$ $``$ $`{\displaystyle \frac{k_BT}{2}}{\displaystyle \underset{𝐤(k<k_0)}{}}{\displaystyle \underset{\omega }{}}{\displaystyle \frac{V_𝐤}{\omega ^2+n_sV_𝐤k^2/m}}`$ (15)
$`\times `$ $`|\mathrm{exp}[i(𝐤𝐫\omega \tau )]\mathrm{exp}[i(𝐤𝐫_1\omega \tau _1)]|^2.`$ (16)
The summation over Matsubara frequencies in Eq. (15) can be carried out explicitly for $`\tau _1=\tau ^+`$, and for $`r|𝐫𝐫_1|\mathrm{}`$ we find

$$\frac{1}{2}[\mathrm{\Phi }(𝐫,\tau )\mathrm{\Phi }(𝐫_1,\tau _1)]^2=\frac{r_s}{4}\left(\frac{n}{n_s}\right)^{1/2}\mathrm{coth}\left(\frac{(n_s/n)^{1/2}}{Tr_s}\right)_{1/r}^{1/L}\frac{dk}{k}[1J_0(kr)].$$
(17)
Here, the temperature $`T`$ is in units of $`e^2/k_B`$, $`L`$ is a length scale of order $`1/k_0`$, and $`J_0(x)`$ is the zero-order Bessel function. Finally, the expression for the mean-square phase fluctuation has the form

$$\frac{1}{2}[\mathrm{\Phi }(𝐫,\tau )\mathrm{\Phi }(𝐫_1,\tau _1)]^2\alpha \mathrm{ln}(r/L),$$
(18)
where the quantity $`\alpha `$ is given by

$$\alpha =\frac{r_s}{4}\left(\frac{n}{n_s}\right)^{1/2}\mathrm{coth}\left(\frac{(n_s/n)^{1/2}}{Tr_s}\right).$$
(19)
In conclusion, the one-body density matrix of the 2D-BCF decays to zero with the power law $`r^\alpha `$ as a consequence of the logarithmic correlations between phase fluctuations. In the limit of zero temperature the decay follows the law $`r^{r_s/4}`$, which is the result found in the QMC study of Magro and Ceperley.

Having reached their result from the finite-temperature formula given in Eq. (18), we can draw two main conclusions: (i) the off-diagonal order reported in the QMC work of Magro and Ceperley proves that the 2D-BCF is a superfluid at $`T=0`$; and (ii) a quantum simulation study of the one-body density matrix at finite temperature would provide for the 2D-BCF an alternative approach to the superfluid density, to be compared with the methods that are currently available for its calculation (see e.g. the method based on the relation between the superfluid density and the mean square winding number). Notice also that, as discussed in the book of Forster, the equivalence between the definition of superfluid density that is provided by the single-particle excitation spectrum and the definition that may be extracted from transverse current correlations and rotation experiments remains a question of continued interest.

As a final point, we should remark that our calculation has not included the possibility of vortex formation in the superfluid as the temperature is raised towards the critical temperature for superfluidity. Equation (18) should therefore be expected to become invalid near $`T_c`$.

## IV ALTERNATIVE APPROACH TO THE ONE-BODY DENSITY MATRIX

We present in this Section an alternative calculation of the asymptotic behavior of the one-body density matrix, following an approach which is closely akin to the resummation method adopted in the book of Popov and in the work by Kagan and coworkers for neutral 2D Bose fluids in a quasicondensate state.

We start from the definition $`\rho (𝐫,𝐫^{})=\psi ^{}(𝐫)\psi (𝐫^{})`$ and use the expression of the field operators in terms of the hydrodynamic density and phase fluctuation operators to find

$$\rho (𝐫,𝐫^{})=n_s\frac{1}{2}\delta (𝐫𝐫^{})+\frac{1}{4n_s}\eta (𝐫)\eta (𝐫^{})+n_s\mathrm{\Phi }(𝐫)\mathrm{\Phi }(𝐫^{}),$$
(20)
to lowest order in the density and phase fluctuations. The transformations carried out in Section II yield

$$\mathrm{\Phi }(𝐫)\mathrm{\Phi }(𝐫^{})=\underset{𝐤(k<k_0)}{}\alpha _𝐤\left[f(E_𝐤)+\frac{1}{2}\right]\mathrm{exp}[i𝐤(𝐫𝐫^{})]$$
(21)
and a similar expression for $`\eta (𝐫)\eta (𝐫^{})`$, with $`f(E_𝐤)=a_𝐤^{}a_𝐤=[\mathrm{exp}(E_𝐤/k_BT)1]^1`$. The final result is

$$\rho (r)=n_s+\underset{𝐤(k<k_0)}{}\left\{\left(\frac{1}{4n_s\alpha _𝐤}+n_s\alpha _𝐤\right)\left[f(E_𝐤)+\frac{1}{2}\right]\frac{1}{2}\right\}\mathrm{exp}[i𝐤𝐫]$$
(22)
where $`r=|𝐫𝐫^{}|`$. This expression has the same form as in the standard Bogoliubov approach, except that the superfluid density $`n_s`$ replaces the condensate density $`n_0`$. Equation (13) can be re-obtained by an approximate resummation of the phase fluctuations to infinite order.

We write the density matrix in the form

$$\rho (r)=\stackrel{~}{\rho }(r)\mathrm{exp}[\mathrm{\Lambda }(r)]$$
(23)
where $`\mathrm{\Lambda }(r)`$ collects all the terms that are responsible for the slow asymptotic decay of $`\rho (r)`$. We find

$`\mathrm{\Lambda }(r)`$ $`=`$ $`{\displaystyle \frac{d^2k}{(2\pi )^2}[1\mathrm{cos}(𝐤𝐫)]\frac{V_𝐤}{2E_𝐤}\left[2f(E_𝐤)+1\right]}`$ (24)
$`=`$ $`{\displaystyle \frac{r_s}{4}}\left({\displaystyle \frac{n}{n_s}}\right)^{1/2}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle \frac{1J_0(xR)}{xg(x)}}\left\{1+{\displaystyle \frac{2}{\mathrm{exp}[2(n_s/n)^{1/2}g(x)/(r_sT)]1}}\right\}`$ (25)
from Eq. (21), where $`R=2[n_s/(nr_s^2)]^{1/4}(r/l_0)`$ and $`g(x)=(1+x^4)^{1/2}`$. We also find

$`{\displaystyle \frac{\stackrel{~}{\rho }(r)}{n}}`$ $`=`$ $`1{\displaystyle \frac{1}{n}}{\displaystyle \frac{d^2k}{(2\pi )^2}[1\mathrm{cos}(𝐤𝐫)]\left\{\frac{k^2}{2mE_𝐤}\left[f(E_𝐤)+\frac{1}{2}\right]\frac{1}{2}\right\}}`$ (26)
$`=`$ $`1r_s\left({\displaystyle \frac{n_s}{n}}\right)^{1/2}{\displaystyle _0^{\mathrm{}}}x𝑑x[1J_0(xR)]\left\{{\displaystyle \frac{x^2}{2g(x)}}+{\displaystyle \frac{x^2}{g(x)(\mathrm{exp}[2(n_s/n)^{1/2}g(x)/(r_sT)]1)}}{\displaystyle \frac{1}{2}}\right\}.`$ (27)
In these equations the temperature is scaled in units of $`e^2/k_B`$.

A power-law decay of $`\rho (r)`$ can now be demonstrated analytically from Eqs. (22)-(24). The function $`J_0(xR)`$ in the integrand in Eq. (23) provides a lower limit of integration going as $`r^1`$, while the upper limit is set by the cut-off wave number $`k_01/L`$ for the applicability of the hydrodynamic Hamiltonian. The asymptotic calculation of the integrals yields

$$\rho (r)\stackrel{~}{\rho }_0(r/L)^\alpha $$
(28)
where the value of the exponent $`\alpha `$ is given by Eq. (18) and we have defined

$$\stackrel{~}{\rho }_0=\underset{r\mathrm{}}{lim}\stackrel{~}{\rho }(r).$$
(29)
In the numerical calculations that are reported in the next Section we find that the ratio $`\stackrel{~}{\rho }_0/n`$ becomes larger than unity at very low temperature, so that this quantity cannot be interpreted as a quasicondensate fraction.

## V NUMERICAL RESULTS

In this Section we report some illustrative calculations of the one-body density matrix of the 2D-BCF, based on Eqs. (22) - (24). We need for this purpose to first evaluate the superfluid density as a function of temperature, for which we adopt the Landau theory based on damping of superfluid flow by emission of collective excitations. That is,

$$\frac{n_s}{n}=1(2nmk_BT)^1\underset{𝐤0}{}k^2\frac{\mathrm{exp}(\beta \omega _𝐤)}{[\mathrm{exp}(\beta \omega _𝐤)1]^2}$$
(30)
with $`\beta =1/(k_BT)`$ and $`\omega _𝐤`$ given by Eq. (11). This yields

$$\frac{n_s}{n}=1\frac{2}{T}_0^{\mathrm{}}𝑑x\frac{x^3\mathrm{exp}[2g(x)/(r_sT)]}{\{\mathrm{exp}[2g(x)/(r_sT)]1\}^2}$$
(31)
with $`T`$ in units of $`e^2/k_B`$. The results for the superfluid fraction $`n_s(T)/n`$ are shown in Figure 1 for some values of $`r_s`$ in the weak-coupling regime.

It is worth noting that, as a consequence of the behavior shown by the superfluid fraction in Figure 1, the single-particle excitation energy $`E_𝐤`$ given in Eq. (9) starts for $`T=0`$ at the collective excitation energy $`\omega _𝐤`$ given in Eq. (11) and decreases with increasing temperature till its leading long-walength term vanishes at the critical temperature $`T_c`$. The plasmon gap in the collective excitation spectrum of Eq. (11) remains instead constant with increasing temperature. This behavior of the plasmon gap in the single-particle spectrum of a charged boson fluid has previously been found in the 3D case within the Hartree-Fock-Bogoliubov approximation, with small deviations from the predictions of the Gavoret-Nozières theorem being present at $`T=0`$.

Our results for the one-body density matrix are reported in Figures 2 and 3, first at $`T=0`$ for several values of $`r_s`$ (Figure 2) and then at $`r_s=0.1`$ for several values of $`T`$ (Figure 3). In each of these Figures the right-hand panel shows a logarithmic plot of the numerical results for $`\rho (r)/n`$ and a linear fit based on Eq. (25), using the value of the exponent $`\alpha `$ given by Eq. (18). The length $`L`$ in Eq. (25), which is not determined by the asymptotic calculation, appears in a logarithmic scale as an additive constant and has been obtained from the linear fitting procedure at large $`r`$. Evidently the power-law prediction is in excellent agreement with the numerical results over an extended range of values for the reduced distance .

## VI SUMMARY

In summary, we have studied superfluidity in a weakly interacting 2D fluid of charged bosons with $`\mathrm{ln}(r)`$ interactions as a function of temperature, using a hydrodynamic reduction of the Hamiltonian that has allowed us to treat the system as a superfluid in the absence of a Bose-Einstein condensate at zero temperature. In the Appendix we show that the absence of a condensate at $`T=0`$ also follows for this system from the standard Bogoliubov approach.

We have shown that the assumption of complete superfluidity for this system at $`T=0`$ accounts quantitatively for the QMC data of Magro and Ceperley on the power-law decay of the one-body density matrix. In turn, computer studies at finite temperature would allow one to determine the behavior of the superfluid fraction and provide a full quantitative test of our theory.

###### Acknowledgements.

This work was partially supported by MIUR under the PRIN-2000 Initiative.

## APPENDIX. ABSENCE OF BOSE-EINSTEIN CONDENSATION IN THE 2D-BCF

We show here that within the Bogoliubov approximation there is no Bose-Einstein condensate in the 2D-BCF both at zero and at finite temperature. In contrast, in the 2D fluid of charged bosons with $`e^2/r`$ interactions a condensate is present at $`T=0`$ and the notion of quasicondensation can be used at $`T0`$.

If we assume that the condensate fraction $`n_0/n`$ is non-zero, then it is easy to show that the expression (9) for the single-particle excitation energy becomes

$$E_𝐤=\sqrt{\frac{2\pi n_0e^2}{m}+\left(\frac{k^2}{2m}\right)^2}$$
(A.1)
in the Bogoliubov approach. The corresponding equation for the condensate fraction reads

$$\left(1\frac{n_0}{n}\right)\left(\frac{n_0}{n}\right)^{1/2}=\frac{r_s}{2}_0^{\mathrm{}}x𝑑x\left\{\frac{h(x)}{2x^2g(x)}1+\frac{[x^2g(x)]^1}{\mathrm{exp}[2(n_0/n)^{1/2}g(x)/(r_sT)]1}\right\}$$
(A.2)
where $`h(x)=1+2x^4`$.

The first term in the square brackets on the RHS of Eq. (A.2) gives a contribution of order $`x^1`$ to the integrand as $`x0`$, and hence the integral diverges. We conclude that $`n_0/n`$ must vanish at any temperature, in agreement with the general sum-rule argument given by Magro and Ceperley. The property $`n_0/n=0`$ is consistent, on the other hand, with the Bogoliubov - des Gennes equations for the Bogoliubov amplitudes $`u_𝐤`$ and $`v_𝐤`$. In this case these equations yield $`v_𝐤=0`$ and $`u_𝐤=1`$, so that the particle density is related to the (now finite) chemical potential by

$$n=\underset{𝐤0}{}\left\{\mathrm{exp}\left[\beta \left(\frac{k^2}{2m}\mu \right)\right]1\right\}^1.$$
(A.3)