text:
# Stellar Masses of Luminous Compact Blue Galaxies at Redshifts 𝑧=0.4-1.2

## 1 Introduction

Luminous Compact Blue Galaxies (LCBGs) refers to a heterogenous class of starburst systems with typical luminosities around $`L^{}`$ that become very common at intermediate redshifts (Koo et al. 1994, 1995; Guzmán et al. 1996, 1997; Phillips et al. 1997; Mallén-Ornelas et al. 1999; Hammer et al. 2001). LCBGs are particularly interesting for studies of galaxy formation and evolution since they have evolved more than any other galaxy class in the last $``$8 Gyrs (Mallén-Ornelas et al. 1999), and have been identified as a major contributor to the observed enhancement of the global star formation rate density of the universe at $`z1`$ (Guzmán et al. 1997; Gruel 2002). Despite their cosmological importance, the nature of LCBGs and their relation to today’s galaxy population still remain largely unknown. The most comprehensive study of LCBGs at intermediate redshift to date –only 45 objects– is that of Phillips et al. (1997) and Guzmán et al. (1997), who concluded that the LCBG class is populated by a mixture of starbursts. About $``$60% of galaxies in their sample are classified as “HII-like” since they are spectroscopically similar to today’s population of vigorously star-forming HII galaxies. The remaining $``$40% are classified as “SB disk-like” since they are more evolved star-forming systems similar to local spirals with a central starburst, and giant irregular galaxies.

Given the diverse nature of the LCBG population, they are unlikely to evolve into one homogeneous galaxy class. There are two popular scenarios. Koo et al. (1994) and Guzmán et al. (1996) have suggested that some subset of the most compact HII-type LCBGs at intermediate redshifts may be the progenitors of local low-mass spheroidal galaxies (or dwarf ellipticals), such as NGC 205. Alternatively, Phillips et al. (1997) and Hammer et al. (2001) have suggested that SB disk-like LCBGs may actually be more massive disks forming from the center outward to become present-day $`L^{}`$ spirals. At the heart of this debate is the question: are the masses of LCBGs comparable to that of today’s massive or low-mass galaxies? Prior measures of Keck velocity widths and HST sizes show that LCBGs have mass-to-light ratios that are only $`0.11.0M_{}/L_B`$, i.e., about 10 times smaller than a typical $`L^{}`$ galaxy today (Koo et al. 1994; Guzmán et al. 1996; Phillips et al. 1997). Since ionized gas emission line widths may not reflect the true gravitational potential, an independent estimate of galaxy mass is desirable. Stellar mass is an obvious choice, since HI and CO emission line widths are too difficult to measure in galaxies beyond the local universe.

Several authors have pointed out that the near-IR luminosity of a galaxy is a robust estimator of its stellar mass (e.g., Aragón-Salamanca et al. 1993; Alonso-Herrero et al. 1996; Charlot 1998). The main uncertainty in the inferred stellar masses arises from the age of the stellar population (Rix & Rieke 1993). Instead of using H$`\alpha `$ to constrain the strength and age of a young burst (Gil de Paz et al. 2000; Östlin et al. 2001; Pérez-González et al. 2003a, 2003b), an alternative approach when studying galaxies at higher redshifts is to use rest-frame far-UV fluxes (Brinchmann & Ellis 2000, hereafter BE00; Papovich, Dickinson & Ferguson 2001). Here we use a modified version of the method described in BE00 to estimate the stellar masses of a sample of LCBGs at intermediate redshift. Our goal is to test previous claims that LCBGs have, on average, an order of magnitude smaller masses than typical galaxies today with similar luminosities. In section 2 we introduce the photometric data set used in this analysis. Section 3 describes the technique developed to derive stellar masses. The results are summarized in Section 4. For consistency with previous work, throughout this paper we adopt Ho = 50 km s-1 Mpc-1 and qo = 0.05.

## 2 The Data

The initial galaxy sample consists of 40 LCBGs at redshifts $`0.4<z<1.2`$ selected from the Flanking Fields around the Hubble Deep Field North, as described in Phillips et al. (1997). A wide variety of information for this sample was made publicly available by the Deep Evolutionary Extragalactic Probe (DEEP) collaboration, including F814W surface photometry, (VF606W \- F814W) colors, Oxygen and Balmer emission line equivalent widths, velocity widths, star formation rates, and virial masses (Phillips et al. 1997; Guzmán et al. 1997). In addition, deep optical Un, G, R photometry for the entire HDF+FF area is available through the Caltech Faint Galaxy Redshift Survey (Hogg et al. 2000), while deep near-IR photometry in the notched HK’ filter is listed in the Hawaii Flanking Field Catalog (Barger et al. 1999). After cross-correlating the various catalogs, we ended up with a final sample of 36 objects with deep Un, G, R, F814W, and HK’ photometry. Two thirds of the objects are classified as HII-like, while the remaining one third are classified as SB disk-like. Limiting magnitudes for the photometry in each of these bands are: Un = 25, G = 26, R = 25.5, F814W = 25.0, HK’ = 20.7. As described in Phillips et al. (1997), all LCBGs in our sample are brighter than F814W=23.74 mag and their half-light radii are smaller than 0.5 arcsec. For F814W and HK’, an estimate of the total magnitude was derived from the total flux within a 3 arcsec diameter aperture. For Un, G, and R, total magnitudes were derived after correcting the magnitudes within 1.7 arcsec diameter apertures for flux outside the aperture considering the point source case, as described in Hogg et al. (2000). The largest errors in the combined data set are due to the comparatively shallower near-IR data (typically $`\sigma _{HK^{}}0.2`$ mag).

## 3 Stellar Mass Estimates

The resulting photometric data set covers Un through HK’ in the observed frame. The red end of this range corresponds to $``$0.9$`\mu `$m in the rest-frame of the highest redshift considered here, and to $``$1.3$`\mu `$m for a more typical $`z0.7`$ galaxy in our sample. The blue end corresponds to rest-frame $``$0.18$`\mu `$m and $``$0.23$`\mu `$m, respectively. Gil de Paz & Madore (2002) have shown that this filter combination is well suited to derive physical properties of galaxies at $`z0.7`$ from broad-band data alone, including star formation timescale, age, metallicity, and stellar mass.

Similarly to the stellar mass estimates method described in BE00, we adopt a grid of evolutionary synthesis models with simple star formation histories (Bruzual & Charlot 2002) and select the most appropriate evolutionary history by error-weighted fitting of the optical and IR photometry. The models have three main parameters: initial mass function (IMF), star formation rate (SFR) history, and metallicity. There are, however, some differences between our method and that described by BE00. Firstly, instead of k-correcting the observed magnitudes for our LCBG sample, we chose to shift all model spectra to the redshift of each of our objects. Redshifted model spectra were convolved with the empirical filter functions corresponding to the actual Un, G, R, F814W, and HK’ broadband filters used in the observations to derive the model magnitudes. These magnitudes were then corrected for extinction assuming not only a range of values for the color excess E(B-V) a la BE00, but also a variety of extinction laws including those derived for the Milky Way (MW; Mathis 1990), Large Magellanic Cloud (LMC; Bouchet et al. 1985), LMC with wavelength-dependent extinction correction (LMC-modified; Mas-Hesse & Kunth, private communication), Small Magellanic Cloud (SMC; Bouchet et al. 1985), and the so-called Calzetti’s law (C94; Calzetti et al. 1994).

Secondly, there is growing evidence that LCBGs, both in the local and distant universe, have a composite stellar population that can be well represented by a young burst superimposed on an older underlying population (Guzmán et al. 1998; Bergvall & Östlin 2001; Östlin et al. 2001; Gruel 2002). We use a composite “underlying+burst” model to fit the observed broadband magnitudes. The burst and underlying components are defined in terms of the assigned star formation rate histories, e.g. instantaneous burst versus constant or exponentially decaying SFR. For simplicity, both populations are assumed to have the same IMF and metallicity. The best fit model is found using a two-step, iterative procedure. Firstly, a single burst model is fitted to the observed Un, G, and R optical bands alone, i.e., the extinction-corrected model magnitudes are compared to the observed magnitudes and the best fit is selected using a simple, weighted, least-squares minimization technique.777No correction for nebular continuum were made to the Un, G, and R magnitudes since its contribution is only significant for burst ages of a few $`10^6`$ yr which are ruled out by the H$`\beta `$ equivalent widths measured for our sample galaxies (Guzmán et al. 1997). The correction for nebular line emission on these broadband magnitudes is estimated to be $`<10`$% (Bergvall & Östlin 2002) and were not applied. Weights were scaled with the flux in each band using the limiting magnitudes to provide the reference minimum weight. Since the observed Un, G, and R bands map approximately the rest-frame UV region of the spectrum at the redshifts of our galaxies, light is entirely dominated by the youngest stars in the burst. The fitted model will provide an initial value for the age of the burst. Secondly, we fit a “underlying+burst” model to all observed magnitudes, including F814W and HK’. The burst component is now restricted to vary by less than 10 timesteps (or a factor $``$2) in age from the initial value. Stellar masses are estimated by normalizing the combined “underlying+burst” spectral energy distribution to the observed magnitudes modulo the IMF. Table 1 summarizes the model ingredients.

Over 800 different model combinations were investigated. As shown by Gil de Paz & Madore (2002), the comparison of broad-band photometry with the predictions of evolutionary synthesis models to derive galaxy physical properties yields, in general, highly degenerate results. A first set of “best fit” models were selected on the basis of having the lowest residuals, i.e., the lowest standard deviations for the five differences between model and observed magnitudes.888During the model fitting procedure, it was found that offsets needed to be applied to the observed magnitudes. We have not been able to identify the source of these offsets. They were derived iteratively as part of the fitting technique. Extensive Montecarlo simulations were performed to ensure that our code could identify and recover artificially introduced zero-point offsets in fake galaxy catalogs that best mimick the properties of our LCBG sample. We are confident that the uncertainty in the derived stellar masses is dominated by model degeneracy and random photometric errors, as explained in the text, and not by offset errors, which contribute only 0.1-0.15 dex to the mass estimates. The median residuals of the “best fit” models are less than 0.15 mag. The largest median residuals are $``$0.35 mag, which correspond to models that use the MW extinction law with E(B-V)=1. On this basis alone, we were able to reject the following cases: (i) models with constant SFR for the underlying component; (ii) models with MW and LMC-modified extinction laws; and (iii) models with LMC and SMC laws and E(B-V)$`>`$0.5. The residuals were found to be only slightly dependent on the choice of IMF, metallicity, and star formation timescales. Since there is growing evidence that the IMF in star-forming and starburst regions is universal and consistent with Salpeter’s (Leitherer 1998), we adopt a Salpeter IMF with lower and upper mass cutoffs of 0.08 and 125 M, respectively, for both the burst and underlying components. We also adopt a metallicity $`Z=0.4Z_{}`$, consistent with the range of metallicities observed in distant LCBGs (Guzmán et al. 1996; Kobulnicki & Zaritski 1999). Finally, we set the star formation timescales to be instantaneous for the burst, and exponential with $`\tau =1`$ Gyr for the underlying component.

The more restrictive choice of parameters reduces significantly the number of “best fit” models. The solutions are, however, still degenerate since various combinations of E(B-V), ranging from 0 to 0.5 mag, and extinction laws (LMC, SMC, and C94), all produce a low residuals fit to the data. Since the timescales for the burst and underlying components in the models are now fixed, we find that the dominant source of degeneracy is the age-extinction degeneracy (Gil de Paz & Madore 2002). This is in the sense that an older burst with less dust-extinction may yield as good fit to the data as a younger but more extincted burst. If the burst is assumed to be less than $`10^8`$ yrs old, as suggested by the relatively strong \[OII\]$`\lambda `$3727 and H$`\beta `$ equivalent widths (Guzmán et al. 1997), then we find only one “best-fit” model that fulfills all the conditions described in our analysis. According to this model, LCBGs are consistent with being experiencing, on average, a 9% burst of star formation that is only 13 Myrs old. This young component is superimposed on an older, redder population that is, on average, 2 Gyrs old. These results are qualitatively in good agreement with the burst strength and age estimates based on mass-to-light ratios and color analysis for a different sample of LCBGs (Guzmán et al. 1998). The derived amount of extinction, E(B-V)=0.5 mag, is also in good agreement with the average value measured in starburst galaxies similar to LCBGs (Gil de Paz et al. 2000; Hammer et al. 2001; Rosa-González et al. 2002). We emphasise, however, that the results on the stellar mass measurements are remarkably independent on the final model choice, as it is discussed below.

Figure 1 shows the distribution of stellar masses for our LCBG sample. The median value of the overall distribution is $`5\times 10^9`$ M. This value is very similar to that derived for nearby LCBGs (Östlin et al. 2001). Note, however, that this value is about $`10`$ times smaller than the average stellar mass in Hammer et al.’s (2001) sample of compact galaxies at intermediate redshifts. This difference is most likely due to the intrinsic higher luminosity of Hammer et al.’s sample, and the fact that these authors adopted a different K-band mass-to-light ratio. The median values for the HII-like and SB disk-like LCBGs are $`5\times 10^9`$ M and $`2\times 10^{10}`$ M, respectively. Although non-significant in this small sample, these values suggest that SB disk-like LCBGs tend to be more massive than HII-like LCBGs. The mean value for the subsample of 17 LCBGs with blue luminosities around L (i.e., $`M_B^{}\pm 1`$) is $`1.8\pm 1\times 10^{10}`$ M. This value is one order of magnitude lower than the stellar mass expected for a local $`L^{}`$ galaxy. Assuming $`M_K^{}=25.1`$ (for Ho= 50 km s-1 Mpc-1; Mobasher, Sharples & Ellis 1993), and a K-band mass-to-light ratio of 1 M/LK,⊙ (Héraudeau & Simien 1997), the estimated stellar mass is $`2\times 10^{11}`$ M. Thus, we conclude that the L CBGs at intermediate redshifts are about 10 times less massive than a typical L galaxy today. This result is in good agreement with previous claims based on virial mass estimates (Koo et al. 1994; Guzmán et al. 1996; Phillips et al. 1997). Our conclusion is not affected by the final choice of the model we adopted here. The median value for the stellar mass estimates derived from the other “best-fit” models differ in a factor $``$2 from the median value given above. The highest values correspond to models using the C94 extinction law, which yield a median stellar mass of $`1.3\times 10^{10}`$ M. We also note that a Salpeter IMF yields systematically higher masses, although only by $`0.2\pm 0.05`$ dex. The use of a different IMF would strengthen our conclusion. Most interestingly, we find that LCBGs at intermediate redshifts span the same range of stellar masses and stellar mass-to-light ratios as that characteristic of Lyman-break galaxies at redshifts $`z3`$ (Papovich, Dickinson & Ferguson 2001).

Figure 2 shows the comparison of the stellar masses with virial masses derived by Phillips et al (1997). The virial masses used for this figure have been corrected for a systematic underestimate of the galaxy gravitational potential when emission-line velocity widths are used. This effect amounts to a $``$30% increase in observed velocity widths (Guzmán et al 1996; Rix et al 1997; Phillips et al 1997; Pisano et al 2001). No corrections for inclination were made. Rotationally-supported, face-on LCBGs are likely to have virial masses significantly lower than their stellar masses. The 1-$`\sigma `$ error bars have been calculated by propagating the magnitude errors, the standard deviations of the fit residuals, and the quoted errors in half-light radii and velocity widths. The largest source of uncertainty is associated to the errors in the shallower HK’ photometry. The average error in the stellar mass measurements amounts to 0.36 dex (in logM). Excluding the two most uncertain measurements, the average of the estimated errors drop to 0.27 dex, i.e., a factor of $``$2. This value is comparable to the median error estimated by BE00 in their error analysis. The average random error for the virial mass measurements is estimated to be 0.18 dex. The median virial-to-stellar mass ratio is 0.39 dex. Thus the stellar mass estimates of LCBGs are, on average, $``$2 times smaller than the virial mass estimates. This result is also in good agreement with the comparison between stellar masses and dynamical masses made by BE00 and Papovich, Dickinson & Ferguson (2001). The scatter in the virial-to-stellar mass ratio is slightly larger than expected from the error analysis, although the exclusion of the two most deviant points results in 0.54 dex, or a scatter of a factor $``$3. No systematic differences in the virial-to-stellar mass ratios are detected between the SB disk-like and HII-like LCBGs. We note, however, that Figure 2 suggests a trend in the sense that the lowest stellar mass LCBGs tend to have the highest virial-to-stellar mass ratios.

## 4 Conclusions

We have presented a new method for determining robust stellar masses for LCBGs up to redshifts $`z1.2`$. Our method expands on previous work in the field by including a two-component “burst+underlying” model population fitted to a combination of rest-frame UV, optical and near-IR broad-band photometry in an iterative fashion. Although specifically tailored for starburst galaxies with a composite stellar population, this method can also be used to estimate the stellar masses of the field galaxies in general of known redshift. According to our best-fit model, LCBGs at intermediate redshifts can be described, on average, as starburst systems which are undergoing a 13 Myrs old burst involving 9% of the galaxy mass superimposed on an underlying, older population (2 Gyrs). This model requires a modest amount of extinction (E(B-V)=0.5) with a LMC extinction law. We have demonstrated via comparisons and simulations that the stellar masses can be estimated with a precision of log $`\mathrm{\Delta }`$Mstar $``$0.3 dex, and are not significantly affected by the main assumptions needed to constrain the degeneracy intrinsic to this method. The median stellar mass for LCBGs in our sample is $`5\times 10^9`$ M, or $``$2.5 times smaller than previous mass estimates based on the virial theorem using scale lenghts and emission-line velocity widths. The new, independent mass estimates are consistent with previous claims that LCBGs at intermediate redshifts are about 10 times less massive than a typical L galaxy today.

RG is grateful to the Institut d’Astrophysique de Paris and the Laboratoire d’Astrophysique at the Observatoire Midi-Pyrenee for their hospitality and financial support for this project. RG also acknowledges funding from NASA grants HF01092.01-97A and LTSA NAG5-11635. G.O. acknowledges support from the Swedish Natural Sciences Research Council and the STINT foundation. We thank Roser Pelló for providing k-corrections used for an independent test of the quality of the photometric data, and Jean Michel Desert for helping in the crosscorrelation of the various catalogs.

## 5 References

* Alonso-Herrero A., Aragón-Salamanca A., Zamorano J. & Rego M. 1996, MNRAS, 278, 417
* Aragón-Salamanca A., Ellis R.S., Couch W.J. & Carter D. 1993, MNRAS, 262, 794
* Barger A.J., Cowie L.L., Trentham N., Fulton E., Hu E.M., Songaila A. & Hall D. 1999, AJ, 117, 102
* Bergvall N. & Östlin G. 2002, A&A 390, 891
* Bouchet et al. 1985, A&A, 149, 330
* Brinchmann J. & Ellis R.S. 2000, ApJL, 536, 77
* Bruzual G. & Charlot S. 2002, in prep.
* Calzetti et al., 1994, ApJ 429, 582
* Charlot S. 1998, in Benvenuti P. et al, eds., AIP Conf. Proceedings Vol. 408, The Ultraviolet Universe at Low and High Redshift: Probing the Progress of Galaxy Evolution. American Institute of Physics, New York, p. 403
* Gil de Paz A., Aragón-Salamanca A., Gallego J., Alonso-Herrero A., Zamorano J. & Kauffmann G. 2000, MNRAS, 316, 357
* Gil de Paz A. & Madore B.F. 2002, AJ, 123, 1864
* Gruel N. 2002, Ph.D. thesis, Observatoire de Paris
* Guzmán, R., Koo, D. C., Faber, S. M., Illingworth, G. D., Takamiya, M., Kron, R., & Bershady, M. A. 1996, ApJ, 460, L5
* Guzmán R., Gallego J., Koo D. C., Phillips A. C., Lowenthal J. D., Faber S. M., Illingworth G. D. & Vogt N. P. 1997, ApJ, 489, 559
* Guzmán R., Koo D. C., Jangren, A., Bershady M., Faber S. M. & Illingworth G. D., 1998, ApJ, 495, L13
* Hammer, F., Gruel, N., Thuan, T. X., Flores, H., & Infante, L., 2001, ApJ, 550, 570
* Héraudeau P. & Simien F. 1997, A&A, 326, 897
* Hogg, D.W., Pahre M.A., Adelberger K.L., Blandford R., Cohen J.G., Gautier T.N., Jarrett T., Neugebauer G. & Steidel C.S. 2000, ApJS, 127, 1
* Kobulnicky H. D., & Zaritsky D. 1999, ApJ, 511, 118
* Koo, D. C., Bershady, M. A., Wirth, G. D., Stanford, S. A., & Majewski, S. R. 1994, ApJL, 427, L9
* Koo, D. C. Guzmán, R., Faber, S. M., Illingworth, G. D., Bershady, M. A., Kron, R., & Takamiya, M. 1995, ApJ, 440, L49
* Leitherer 1998, in Gilmore G., & Howell D. eds., 38th Herstmonceux Conference, The Stellar Initial Mass Function. ASP Conference Series, Vol. 142, p.61
* Leitherer et al. 1999, ApJS, 123, 3
* Mallén-Ornelas, G., Lilly, S. J., Crampton D., & Schade D. 1999, ApJ, 518, 83
* Mathis 1990, ARA&A, 38
* Mobasher B., Sharples R.M. & Ellis R.S. 1993, MNRAS, 263, 560
* Östlin G., Amram P., Bergvall N., Masegosa J., Boulesteix J. & Márquez I. 2001, A&A 374, 800
* Papovich C., Dickinson M. & Ferguson, H.P. 2001, ApJ, 559, 620
* Pérez-González P.G., Gil de Paz A., Zamorano J., Gallego J., Alonso-Herrero A. & Aragón-Salamanca A. 2003a, MNRAS, 338, 508
* Pérez-González P.G., Gil de Paz A., Zamorano J., Gallego J., Alonso-Herrero A. & Aragón-Salamanca A. 2003b, MNRAS, 338, 525
* Phillips A. C., Guzmán R., Gallego J., Koo D. C., Lowenthal J. D., Vogt N. P., Faber S. M. & Illingworth G. D. 1997, ApJ, 489, 543
* Pisano D.J., Kobulnicky H. D., Guzmán R., Gallego J. & Bershady, M. A. 2001, AJ, 122, 1194
* Rix H.W. & Rieke M.J. 1993, ApJ, 418, 123
* Rix H.W., Guhathakurta P., Colles M. & Ing K. 1997, MNRAS, 285, 779
* Rosa-González D., Terlevich E. & Terlevich R. 2002, MNRAS, 332, 283


format: text

text:
# Negative heat capacity of sodium clusters

## I Introduction

Negative microcanonical heat capacity in atomic and molecular clusters was theoretically predicted by considering simple models of the distribution of local minima that characterize the potential energy landscape (PEL) of clusters Bixon . In that study, it was found that for high values of the parameter involving the ratios of the vibrational frequencies corresponding to the global and local isomers, the caloric curve displays an S-shaped loop, with a negative heat capacity in the vicinity of the melting point Bixon . In another study on the solid-liquid transition of clusters Labas , it was shown that in microcanonical simulations of Lennard-Jones clusters, an increase in total energy causes a temperature reduction. This effect was related to the broadening of the cluster kinetic-energy distribution toward lower energy values Labas .

Although the existence of negative heat capacity in physical systems like stars or star clusters Thirring ; Bell , and in fragmenting nuclei Gross ; Agostino  is well-known, this peculiar effect gained a lot of interest in the field of atomic and molecular clusters due to recent experimental results where a negative heat capacity was measured for a 147-atom sodium cluster Schmidt . In this study, the photofragmentation mass spectra was used to measure the internal energy of free, mass selected clusters with known temperature. These measurements were used to determine the microcanonical caloric curve of the Na+147 that shows the characteristic S-shaped (backbending) feature, indicating a negative heat capacity Schmidt . The negative value of the microcanonical heat capacity was interpreted by considering that a finite system upon melting tries to avoid partly molten states and prefers to convert some of its kinetic energy into potential energy Schmidt ; Schmidt1 . This peculiarity has been attributed to the non-additivity of the total energy of a cluster with finite size Schmidt ; Schmidt1 .

Microcanonical heat capacities of metal clusters have been theoretically investigated using constant-energy molecular dynamics (MD) with many-body potentials Gallego ; Calvo ; Jell  and an orbital-free version of the first-principles MD method Agua1 ; Agua2 . In these studies, heat capacities of fcc transition and noble metal clusters with up to 23 atoms were calculated to characterize their melting-like transition Gallego . In another study, on the melting of sodium clusters Calvo , the microcanonical caloric curve of the Na55 cluster was obtained. However, in not one of these studies was a signature of a negative heat capacity found. Similar results, indicating the non-existence of a negative heat capacity in constant-energy orbital-free first-principles MD simulations of larger sodium clusters (Na55, Na92, and Na142), were obtained Agua2 . Nevertheless, in such calculations the simulation time employed was too short to obtain converged results Agua2 . On the other hand, in microcanonical MD simulations of AlN, $`N`$ = 7, 13, 55, and 147, clusters, a negative heat capacity was obtained for the larger Al147 cluster Jell .

In the present work, motivated by the availability of experimental techniques that allow the measurement of the microcanonical heat capacity and other thermal properties of mass selected metal clusters Haber1 ; Haber2 ; Haber3 , we theoretically investigate the behavior of the heat capacity of sodium clusters in the size range of 13-147 atoms. In our approach, constant-energy MD simulations are performed using a phenomenological many-body potential that mimics the metallic bonding of sodium clusters. This approximation allows us to use simulation times of the order of $``$ 50 ns, in order to obtain converged averages of the microcanonical heat capacity and other cluster thermal properties. Our main objective is to gain additional insights into the conditions that determine if a cluster has a negative heat capacity. The main finding of this work shows that the width of the distribution function of the kinetic energy and the spread of the distribution of the potential energy minima (isomers), characterizing the PEL, are useful features to determine the signature of the cluster heat capacity. In section II, we provide the theoretical background on which this study is based. The results and their discussion are given in section III, and section IV contains a summary and the conclusions of this work.

## II Theoretical background

The heat capacity and temperature of sodium clusters as a function of the cluster total energy are calculated through constant-energy MD simulations using the microcanonical expressions derived in Refs. Pearson, ; Sugano, :

$$\frac{C}{Nk_B}=\left[NN(1\frac{2}{3N6})KK^1\right]^1$$
(1)
$$T=\frac{2K}{(3N6)k_B},$$
(2)
where $`K`$ is the kinetic energy of the cluster, $`k_B`$ is the Boltzmann constant, and $``$$``$ denotes a time average. In these formulas, 3$`N`$ was changed to 3$`N`$-6, the number of degrees of freedom of the system, since the calculations are performed for a non-translating and non-rotating cluster in a three-dimensional space (the position of the center of mass was fixed and the total momentum was held to zero during the simulations).

In our implementation of the constant-energy MD method, the Newton’s equations of motion are solved with the Verlet algorithm Verlet  using a time step of 2.4 fs, which provides total energy conservation within 0.001 $`\%`$. A typical calculation consists in heating up a cluster from its lowest-energy solid-like configuration until it transforms into a liquid-like cluster. To simulate this procedure the cluster total energy is increased in a step-like manner by scaling up the atomic velocities. For each initial condition the cluster was equilibrated during 104 time steps and the time averages of the physical quantities are calculated using 107 time steps. This averaging time is increased by a factor of 2 when the cluster is in the region of the solid-to-liquid transition in order to ensure the calculation of fully converged averages.

To model the metallic bonding in sodium clusters we used the many-body Gupta potential Gupta , which is based on the second moment approximation of a tight-binding hamiltonian Cleri . Its analytical expression is given by:

$$V=\underset{i=1}{\overset{N}{}}V_i$$
(3)
$$V_i=A\underset{ji}{}e^{p\left(\frac{r_{ij}}{r_0}1\right)}\xi \left(\underset{ji}{}e^{2q\left(\frac{r_{ij}}{r_0}1\right)}\right)^{\frac{1}{2}}$$
(4)
where $`r_0`$, $`A`$, $`\xi `$, $`p`$, and $`q`$ are adjustable parameters Cleri . For sodium clusters these parameters have been fitted to band structure calculations Li . Their values are: $`A`$=0.01595 eV, $`\xi `$=0.29113 eV, $`r_0`$=6.99 bohr, $`p`$=10.13, and $`q`$=1.30 Li . This phenomelogical many-body potential has been used to study the melting-like transition in sodium clusters of different sizes using Monte Carlo (MC) Calvo  and constant-energy MD simulations Reyes . A good qualitative agreement has been obtained between structural and thermal properties calculated using the Gupta potential Calvo ; Reyes  and those generated from first-principles methods Agua1 ; Agua2 . An additional advantage in using this potential is that it allows simulation times of the order of 50 ns, necessary to obtain fully converged time averages in the melting-like transition region.

## III Results and Discussion

The microcanonical heat capacities of the NaN, $`N`$ = 13, 20, 55, 135, 142, and 147, clusters, calculated using Eq. (1), are displayed in Fig. 1. For $`N`$ = 13, 20, and 55, they are continuous functions of the cluster total energy showing a maximum value that is characteristic of a melting-like transition Calvo ; Agua1 ; Agua2 ; Reyes . On the other hand, the heat capacity of the larger clusters ($`N`$ = 135, 142, and 147) shows two discontinuity points and a continuous negative-valued interval between them. This peculiar behavior in the heat capacity is related with a backbending loop in the caloric curve (temperature as a function of the total energy) Bixon ; Labas ; Wales1 ; Wales . In fact, our calculated microcanonical caloric curves of Na135, Na142, and Na147 show the backbending loop at the same energies where the heat capacity takes negative values (see the caloric curves shown in Fig. 2 of Ref. Reyes, ). In previous studies, the negative slope (backbending loop) of the microcanonical caloric curve has been attributed to a peculiar behavior of the cluster entropy as a function of energy that shows a dent with inverted curvature in the region of the solid-liquid transition Bixon ; Labas ; Schmidt ; Schmidt1 ; Wales .

In the present work, we analyze the behavior of the microcanonical heat capacity of sodium clusters from a different perspective. First, we consider Eq. (1) and obtain the condition to have a negative value in the heat capacity:

$$Z_EKK^1>\frac{3N6}{3N8}.$$
(5)
Figure 2 display the values of $`Z_E`$ (black dots) as a function of the cluster total energy $`E`$ that were calculated from a time average, using the MD trajectories. In the same scale the threshold value $`Z_c(3N6)/(3N8)`$ for each cluster size is given. In Fig. 2 , it can be graphically seen how the relative difference between $`Z_E`$ and $`Z_c`$ changes with the cluster size. For the three smaller clusters (see panels (a), (b), and (c) in Fig. 2 ) the $`Z_E`$ values do not overcome the threshold value, whereas for the three larger clusters there is a range of total energy where $`Z_E`$ satisfy the condition to have negative heat capacities (see panels (d), (e), and (f) in Fig. 2 ).

In order to investigate what determines a negative value of the heat capacity, we consider the quantity $`Z_E`$. This is the product of the averages of the kinetic energy and of the inverse of this quantity, and therefore, its value will depend on the distribution function of the kinetic energy, $`g_E(K)`$. The average of any function of the kinetic energy $`f(K)`$ can be obtained through the following expression:

$$f_E_g=f(K)g_E(K)𝑑K.$$
(6)
Since the distribution function of the kinetic energy, $`g_E(K)`$, determines the behavior of $`Z_E`$, it is useful to analyze $`g_E(K)`$ at different values of the cluster total energy $`E`$. The calculation of this quantity is straightforward from the constant-energy MD simulations. Figure 3 shows $`g_E`$ as a function of the normalized mean deviation $`\delta K`$ = ($`K`$ \- $`K`$)/ $`K`$ for three different energies, corresponding to the cases where the cluster is in the solid- and liquid-like phases, and at the middle of the melting-like transition. The analysis of $`g_E`$ as a function of $`\delta K`$, instead of $`K`$, has the advantage that it allows the comparison, on the same scale, of the lineshapes of this function at different cluster energies and for different cluster sizes.

As a general trend, it is found that $`g_E`$ becomes narrower for increasing cluster sizes, indicating a larger relative dispersion of the kinetic energy values for the smaller clusters. This result is expected since it confirms the increment of fluctuations in kinetic energy of a physical system that decrease in size. A common characteristic of $`g_E`$, existing in the six clusters investigated, is the larger broadening of the distribution function when the cluster is at the melting-like transition. At lower (solid-like phase) and higher (liquid-like phase) energies, the width of $`g_E`$ is smaller, whereas at the phase transition the fluctuations in kinetic energy, as expected, should increase.

For the three smaller clusters which do not have negative heat capacity, $`g_E(\delta K)`$ shows a nearly symmetric lineshape independent of the cluster energy (see panels (a), (b), and (c) of Fig. 3 ). In contrast, Na142 and Na147, that show a negative heat capacity, have a distribution function $`g_E(\delta K)`$ with a shoulder towards positive values of $`\delta K`$, at energies in the middle of the melting region (see panels (e) and (f) of Fig. 3 ). Although this difference in the distribution function of the kinetic energy could be a useful feature to determine the existence of a negative heat capacity, the Na135 cluster would be an exception to this rule since its $`g_E(\delta K)`$ does not show a resolved shoulder in its lineshape (see panel (d) in Fig. 3), but it has negative heat capacity.

On the other hand, a characteristic of $`g_E`$ that would be useful to determine the sign of the heat capacity is the width of the distribution function, which can be obtained through its second moment:

$$(\delta K)^2_g=(\delta K)^2g_E(\delta K)d(\delta K).$$
(7)
The second moment of $`g_E`$ corresponds to the second term in the expansion of $`Z_E`$, which can be obtained from the left hand term of Eq. (5) Pearson ; Sugano :

$$Z_E=1+(\delta K)^2_g+\mathrm{}.$$
(8)
By taking terms up to second order in this expansion, assuming that $`\delta K_g<<1`$, $`Z_E`$ can be approximated by:

$$Z_EZ_E^{(2)}=1+(\delta K)^2_g.$$
(9)
Since the $`Z_E^{(2)}`$ values can be calculated using $`g_E`$ and Eq. (7), it is possible to check the validity of this approximation, which can only be applied to systems with a finite number of particles. Fig. 2 shows the values of $`Z_E^{(2)}`$ (stars) as a function of the cluster energy. It can be seen that the difference with $`Z_E`$ is small for the three smaller clusters and negligible for the three larger ones. Then, $`Z_E^{(2)}`$ can be considered as a quantitative measure of the width (second moment) of the distribution function of the kinetic energy, and can be used to determine the sign of the heat capacity. Figure 4 shows the maximum values of $`Z_E^{(2)}`$ (black dots) calculated for energies at the middle of the melting transition, and their comparison with the threshold values $`Z_c`$ as a function of the cluster size (full line). From this figure, it is obvious that although the width of $`g_E`$ is relatively large for the three smaller clusters, the corresponding maximum values of $`Z_E^{(2)}`$ are below the threshold $`Z_c`$, and therefore these clusters do not show a negative heat capacity. On the other hand, the width of $`g_E`$ for Na135, Na142, and Na147 is smaller than for Na13, Na20, and Na55, however, $`Z_c`$ is a faster decreasing function of the cluster size $`N`$, such that the maximum of $`Z_E^{(2)}`$ lies above the threshold values, indicating that the larger clusters have negative heat capacities. Then, the above results suggest that the width of the distribution function of the kinetic energy is a useful property to determine the sign of the heat capacity of clusters. However, this quantity not only depends on the cluster size, but also on the characteristics of the PEL.

To illustrate the importance of the topology of the PEL we have investigated the behavior of the heat capacity of 55-atom clusters using the many-body Gupta potential Gupta , shown in Eqs. (3) and (4), for the different metals listed in Tables I and III of Ref. Cleri, , and in Table II of Ref. Li, . Our results show that the Cd55 cluster (with the following parameter values Cleri : $`A`$=0.0416 eV, $`\xi `$=0.4720 eV, $`p`$=13.639, and $`q`$=3.908) has a negative heat capacity, but none of the other 55-atom clusters show this peculiarity. The upper insets of Fig. 4 show the calculated caloric curve with the corresponding backbending loop and the heat capacity with negative values for a range of total energy values of the Cd55 cluster. The inset at the lower right corner of Fig. 4 shows $`g_E`$ as a function of the normalized mean deviation of the kinetic energy for both, the Na55 and Cd55, clusters. It can be seen that the broadening of $`g_E`$ is larger in Cd55 than in Na55, such that, it generates a maximum value of $`Z_E^{(2)}`$ (this value corresponds to the point represented by a star in Fig. 4) that overcome the threshold $`Z_c`$, and consequently, the Cd55 cluster display a negative heat capacity. This comparison with the Na55 cluster which does not show this peculiarity, indicates that although both clusters have the same size, their dynamical properties defined by their corresponding PEL’s, generate different behavior in their heat capacities.

In order to investigate the influence of the PEL on the different widths of the distribution function of the kinetic energy of the Na55 and Cd55 clusters, further studies are necessary. In this direction, the calculation of short-time averages of the kinetic energy and periodic quenchings of instantaneous configurations during the MD trajectories allow us to obtain the distribution of potential energy minima (isomers) that are accessible at different cluster energies Wales1 . Figure 5 display the normalized distribution of potential energy minima, obtained by periodical quenchings using MD trajectories at a total energy where the cluster is at the middle of the melting transition, for the Na55 and Cd55 clusters. It can be notice that the number of isomers with higher energy relative to the global minimum, are larger for the Cd55 cluster in comparison with the results obtained for Na55. This result can be explained by taking into account that the range of the interatomic forces is shorter in Cd than in Na clusters, mainly due to the higher value of the $`q`$ parameter in the many-body Gupta potential Karo3 . The physical reason for the larger number of minima at short range is the loss of accessible configuration space as the potential wells become narrower, thus producing barriers where there are none at long range Doye .

To show how the distribution of potential energy minima determine the broadening of the distribution function of the kinetic energy, we approximate the complex topology of the PEL by a set of independent harmonic potential wells in the $`3N6`$ dimensional space. Each one of these wells is associated to the different potential energy minima forming the distribution of isomers shown in Fig. 5. For each potential energy minimum denoted by $`l`$, the distribution function of the kinetic energy at a total energy $`E`$, in the harmonic approximation, is given by Labas ; Pearson :

$$g_{E,l}(K)=C_l(E\mathrm{\Delta }_lK)^{\frac{3N7}{2}}K^{\frac{3N7}{2}},$$
(10)
where $`\mathrm{\Delta }_l`$ is the potential energy of the isomer $`l`$, relative to the potential energy value of the lowest-energy isomer, and $`C_l`$ is a normalization constant such that:

$$_0^{E\mathrm{\Delta }_l}g_{E,l}(K)𝑑K=1.$$
(11)
The distribution function of the kinetic energy $`g_{E,har}`$, at a total energy $`E`$, corresponding to the whole PEL can be constructed by adding up the contribution of each harmonic potential well, weighted by the probability, $`\omega _{E,l}`$, of finding a given isomer during the quenching from the MD trajectories. This probability is given by the height of the distribution shown in Fig. 5. Then, $`g_{E,har}`$ is given by:

$$g_{E,har}=\underset{l=1}{\overset{l_{max}}{}}\omega _{E,l}g_{E,l}(K),$$
(12)
with

$$\underset{l=1}{\overset{l_{max}}{}}\omega _{E,l}=1.$$
(13)
By using the data from the whole distribution of isomers in Fig. 5, the distribution function of the kinetic energy $`g_{E,har}`$ was calculated for the Cd55 and Na55 clusters. They are displayed in the insets of Fig. 5 (full lines). A comparison between $`g_{E,har}`$ and the exact $`g_E`$ (obtained from the MD simulation and displayed in the lower right inset of Fig. 4) shows that there is a good agreement between the two distribution functions. This indicates that $`g_E`$ is determined mainly from the number of isomers and the probability to found them (height of the distribution), rather than from other features of the PEL like saddle points. The advantage in introducing $`g_{E,har}`$ in this discussion is related with the fact that it is possible to analyze the broadening of this distribution function of the kinetic energy by considering different subsets of potential energy minima (isomers). This is useful to determine what regions of the PEL are more relevant to increase the width of $`g_{E,har}`$, and investigate the appearance of the negative heat capacity. The insets of Fig. 5 show three partial distribution functions $`g_{E,har}`$, considering different subsets of isomers corresponding to three intervals of low (L), medium (M) and high (H) potential energy values. By analyzing the relative contribution of these subsets to the width of $`g_{E,har}`$ for the Cd55 and Na55 clusters, it is found that the larger broadening in the cadmium cluster is mainly due to the larger contribution of the isomers in the range of high potential energy which are spreaded along a larger interval than those corresponding to the Na55 cluster. The width of $`g_{E,har}`$ for the Na55 cluster is smaller since there are proportionally less isomers with high potential energy, and they are extended over a shorter interval of values. As was mentioned above, the physical reason for this difference in the distribution of isomers between the Cd55 and Na55 clusters is the shorter range of the many-body forces existing in the cadmium cluster as compared with those present in the sodium cluster. A similar result was obtained for 55-atom clusters using a pairwise Morse potential for different values of the range of the interatomic forces Doye . In that case the backbending loop in the caloric curve (negative heat capacity) was obtained using a Morse potential with a range of the interatomic forces that is shorter than the range characteristic of alkali metals which have long-ranged interactions Doye .

Therefore, if a detailed characterization of the distribution of isomers forming the PEL of clusters is performed, the broadening of $`g_E`$ may be estimated, and by the comparison of the corresponding $`Z_E^{(2)}`$ and $`Z_c`$ values, it would be possible to predict the sign of the heat capacity of clusters.

## IV Summary

The microcanonical heat capacity of sodium clusters has been calculated using constant-energy MD simulations and the many-body Gupta potential. Negative values for the heat capacity at energies where the cluster is at the melting-like transition were found for Na135, Na142, and Na147. The smaller sodium clusters NaN, $`N`$=13, 20, and 55, do not show this peculiarity. An analysis of the calculated distribution function of the kinetic energy $`g_E`$ for the six clusters investigated, shows that the width of this distribution function is a useful feature to determine the sign of the heat capacity. It was found, that although the broadening of $`g_E`$ is larger for the smaller clusters, it is not enough to overcome the corresponding threshold value to obtain a negative heat capacity. However, since this threshold is a fast decreasing function of the cluster size, the broadening of $`g_E`$ in the larger clusters is enough to generate a negative heat capacity.

It was also shown that the broadening of $`g_E`$ depends on the distribution of potential energy minima that characterize the PEL of clusters. Specifically, as the range of the many-body interactions is decreased (like in the case of Cd clusters), the number of local minima with higher energy increases generating a larger broadening in $`g_E`$, and consequently a negative heat capacity. The analysis presented in this paper shows how the complex topology of the PEL can be explored to extract the main features that determine the sign of the heat capacity of metal clusters.

###### Acknowledgements.

This work was supported by Conacyt-Mexico under Project No. G32723-E. JARN acknowledges a graduate fellowship from DGEP-UNAM.


text:
# Quantum magneto-oscillations in a two-dimensional Fermi liquid

## Abstract

Quantum magneto-oscillations provide a powerfull tool for quantifying Fermi-liquid parameters of metals. In particular, the quasiparticle effective mass and spin susceptibility are extracted from the experiment using the Lifshitz-Kosevich formula, derived under the assumption that the properties of the system in a non-zero magnetic field are determined uniquely by the zero-field Fermi-liquid state. This assumption is valid in 3D but, generally speaking, erroneous in 2D where the Lifshitz-Kosevich formula may be applied only if the oscillations are strongly damped by thermal smearing and disorder. In this work, the effects of interactions and disorder on the amplitude of magneto-oscillations in 2D are studied. It is found that the effective mass diverges logarithmically with decreasing temperature signaling a deviation from the Fermi-liquid behavior. It is also shown that the quasiparticle lifetime due to inelastic interactions does not enter the oscillation amplitude, although these interactions do renormalize the effective mass. This result provides a generalization of the Fowler-Prange theorem formulated originally for the electron-phonon interaction.

Patterns of quantum magneto-oscillations in thermodynamic (de Haas-van Alphen effect) and transport (Shubnikov-de Haas effect) quantities encode three important parameters of a Fermi-liquid metal. The period of the oscillations gives the area of the extremal cross-section of the Fermi surface, the slope of the temperature dependence of the oscillations amplitude provides the quasiparticle effective mass, and the phase shift between oscillations of spin-up and-down electrons yields the (renormalized) spin susceptibility. Magneto-oscillations studies of the Fermi-liquid (FL) state in two-dimensions (2D) date back to the early 70s, when semiconductor heterostructures first became available . Another surge of the activity in this field, which occurred in mid 90s, was stimulated by the discovery of the metallic state at $`\nu =1/2`$ . Recently, Shubnikov-de Haas oscillations have been used to determine the parameters of the “anomalous” metallic state in Si MOSFETs and other semiconductor heterostructures exhibiting an apparent metal-insulator transition in zero magnetic field . Despite the long and successful history of quantifying FLs in 3D via magneto-oscillations, this method remains controversial in 2D. The primary goal of our paper is to resolve some of the open issues.

The first controversy is related to the applicability of the current theory of magneto-oscillations to the 2D case. The analysis of the experimental data in 2D is often based on the premise that the classic result for magneto-oscillations in a Fermi liquid for the 3D case, known as the “Lifshitz-Kosevich (LK) formula” , is transferrable to 2D upon a trivial change in the electron spectrum. The crucial features of the LK formula, i.e., its validity for arbitrarily strong interactions (without destroying the FL) and the fact that the FL parameters entering the formula are taken at zero magnetic field, survive on this premise. The deviations of the observed oscillation pattern in stronger fields from that predicted by the LK formula are ascribed to oscillations in the effective $`g`$-factor  and the effective mass . On the other hand, there have been warnings that the 3D LK formula is non-transferrable to 2D  for any field strength. Hence the situation needs to be clarified. The second controversy–not specific to 2D–is related to the effect of quasiparticle damping. It is often mentioned in the literature that any scattering of quasiparticles, elastic and inelastic, contributes to the smearing of magneto-oscillations via the effective Dingle temperature (scattering rate) for a given process. Alternatively, Fowler and Prange  showed that the electron-phonon scattering rate does not appear in the oscillations amplitude due to the cancellation of two $`T`$dependent parts of the Matsubara self-energy (cf. also ). To the best of our knowledge, this cancellation has never been discussed for other interactions, including the electron-electron one, which is of a primary importance for 2D electron systems. The last issue to be addressed in this paper (and not discussed previously in the literature) is the effect of interference between electron-impurity and electron-electron scattering on magneto-oscillations, neglected in the LK formula. The theory of interference effects in the ballistic regime , when $`T\tau 1,`$ where $`\tau `$ is the electron-impurity scattering time (we set $`\mathrm{}=k_B=1`$ throughout the paper), offers a plausible explanation of the metallic temperature dependence in the metallic phase of the 2D metal-insulator transition. Unusual (within the LK framework) temperature dependences of the oscillation amplitude are also commonly observed in Si MOSFETs , but the proper theory is currently lacking.

Our answers to these open questions are as follows. i) Although it is true that the LK formula does not work in 2D at $`T=0`$ and in the absence of disorder, it is still applicable to the situation when finite temperature and/or disorder cause the oscillations to be exponentially small. ii) The cancellation of the scattering rate term in the Matsubara self-energy is pertinent to any inelastic interaction, including the electron-electron one. Due to this cancellation, the scattering rate of inelastic processes does not enter the oscillation amplitude. iii) Interference between electron-impurity and electron-electron interactions gives a new, $`T\mathrm{ln}T`$, dependence of the amplitude’s argument, that can be interpreted equivalently either as a “$`T`$dependent” effective mass or Dingle temperature. The functional form of this dependence is the same in the diffusive $`\left(T\tau 1\right)`$ and ballistic $`\left(T\tau 1\right)`$ regimes.

We limit our analysis to the de Haas-van Alphen effect and assume the chemical potential is fixed. The main features of the results for the de Haas-van Alphen effect, in particular, the $`T`$-dependence of the oscillation amplitude, are commonly expected to apply to the Shubnikov-de Haas effect as well, although a rigorous proof of that is currently lacking. Assuming a fixed chemical potential is not essential for the case of small oscillations, which is the focus of this paper (see below). We begin with a brief reminder of how the LK formula is derived in the Luttinger formalism . The key issue here is whether the zero-field FL parameters determine uniquely the oscillation pattern in a finite (and not small) field. The (Matsubara) self-energy (that encodes all FL parameters) has consists of parts $`\mathrm{\Sigma }=\mathrm{\Sigma }_0+\mathrm{\Sigma }_{\text{osc}},`$ where $`\mathrm{\Sigma }_0`$ may contain a monotonic (nonoscillatory) dependence on magnetic field and $`\mathrm{\Sigma }_{\text{osc}}`$ oscillates with the field. For electron-electron interactions in 3D, $`|\mathrm{\Sigma }_{\text{osc}}|/|\mathrm{\Sigma }_0|N^{3/2}`$, where $`N`$ is the number of occupied Landau levels, whereas the leading term in the oscillatory part of the thermodynamic potential $`\mathrm{\Omega }`$ falls off as $`N^{5/2}`$ for $`\mathrm{\Sigma }_{\text{osc}}=0`$. Expanding $`\mathrm{\Omega }`$ in $`\mathrm{\Sigma }_{\text{osc}}`$ up to the second order (the first order term vanishes due to the property $`\delta \mathrm{\Omega }/\delta \mathrm{\Sigma }=0`$), one finds that the oscillatory part of $`\mathrm{\Sigma }`$ can always be neglected in the semiclassical regime ($`N1`$). With this simplification and for a momentum-independent self-energy, the amplitude of the $`k^{\text{th}}`$ harmonic in $`\mathrm{\Omega }`$ is given by

$$A_k=\frac{4\pi ^2kT}{\omega _c}\underset{\epsilon _n>0}{}\mathrm{exp}\left(\frac{2\pi k\left[\epsilon _n+i\mathrm{\Sigma }_0(i\epsilon _n,T)\right]}{\omega _c}\right),$$
(1)
where $`\epsilon _n=\pi (2n+1)T`$ and $`\omega _c=eB/mc`$. Notice that $`i\mathrm{\Sigma }_0`$ is real and does not contain a constant term. The second argument of $`\mathrm{\Sigma }_0`$ emphasizes the fact that the temperature enters $`\mathrm{\Sigma }_0`$ in two ways: via the Matsubara frequency and via the thermal distribution of electrons and other degrees of freedom. For a generic Fermi liquid and in the presence of short-range impurities, $`i\mathrm{\Sigma }_0(i\epsilon _n,T)=\alpha \epsilon _n+\text{sgn}\epsilon _n/2\tau `$, so that the effective mass is defined as $`m^{}=m(1+\alpha )`$. The amplitude then assumes a familiar form

$$A_k=\frac{2\pi ^2kT/\omega _c}{\mathrm{sinh}\left(2\pi ^2kT/\omega _c^{}\right)}\mathrm{exp}\left(\frac{2\pi ^2T_D}{\omega _c}\right),$$
(2)
where $`\omega _c^{}=eB/m^{}c`$ and $`T_D=1/2\pi \tau `$ is the Dingle temperature. Momentum-dependence of $`\mathrm{\Sigma }_0`$ of the form $`\beta v_F\left(pp_F\right)`$ results in a change of the effective mass in (2) to $`m^{}=m\left(1+\alpha \right)(1+\beta )^1`$ and in multiplying (2) by $`Z_sm^{}/m`$, where $`Z_s`$ is the renormalization factor.

In arbitrary dimensionality $`D`$, the estimates for the ratio of oscillatory to monotonic-in-field parts of the self-energy and for the leading oscillatory term in $`\mathrm{\Omega }`$ change to $`N^{D/2}`$ and $`N^{(D+2)/2}`$, respectively. For $`D=2`$, the oscillations in the self-energy are as important as in the thermodynamic potential itself . The Luttinger expansion at $`T=T_D=0`$ breaks down and the LK formula is not, generally speaking, valid . The physical reason is that the ground states of an interacting system at $`B=0`$ and in a finite field are not adiabatically connected in 2D. This fact has been emphasized by recent findings that the ground state of a 2D electron liquid is not a Fermi liquid even for $`N1`$, but rather a charge-ordered state . Nevertheless, an absense of the full LK formula in 2D does not preclude a canonical analysis of magneto-oscillations, if under more restrictive conditions, as the FL-behavior is restored at higher energies.

The power-counting argument for (against) neglecting the oscillatory part of the self-energy in 3D (2D) is valid at $`T=T_D=0`$. If the real and/or Dingle temperatures are sufficiently high, i.e.,

$$2\pi ^2\left(T/\omega _c^{}+T_D/\omega _c\right)1,$$
(3)
the amplitudes of all oscillatory quantities, including the self-energy, are exponentially small. Neglecting the oscillatory part of the self-energy, the amplitude of the first harmonic takes the form

$$A_1=\left(4\pi ^2T/\omega _c\right)\mathrm{exp}\left(2\pi \left[\pi T+i\mathrm{\Sigma }_0(i\pi T,T)\right]/\omega _c\right),$$
(4)
where, due to condition (3), we limited the Matsubara sum by the first term $`\epsilon _0=\pi T`$. The oscillatory part of $`\mathrm{\Sigma }`$ results in a correction to $`A_1`$ which is itself of order $`A_1`$ (with exponential accuracy). The net contribution to $`\mathrm{\Omega }`$ is of order $`A_1^2`$, which is of the same order as $`A_2`$ for $`\mathrm{\Sigma }_{\text{osc}}=0`$. Thus, harmonics with $`k2`$ are affected by the oscillations in $`\mathrm{\Sigma }`$ and a 2D analog of the LK formula, which includes the sum over all $`k`$, can only be derived in a perturbation theory for a specific interaction but not for a generic FL. However, the $`k=1`$ harmonic does not include $`\mathrm{\Sigma }_{\text{osc}}`$ and, as long as (3) is satisfied, the analysis can proceed as in the 3D case. In what follows, we assume that (3) is satisfied and the amplitude of the first (and only important) harmonic is given by (4).

Next we discuss whether the quasiparticle relaxation rate affects the amplitude of magneto-oscillations. We set $`T_D=0`$ temporarily. Suppose that a quasiparticle relaxation rate is measured in a clean Fermi liquid, e.g., via electron heat conductivity, with a result that $`1/\tau _{\text{e-e}}T^2`$. It seems natural to assume that the same rate contributes also to the Dingle temperature of magneto-oscillations. That this is not the case was shown for the electron-phonon interaction by Fowler and Prange . Here we generalize their arguments for the electron-electron interaction in 3D and 2D, and then give a general theorem for an arbitrary interaction. For a generic FL in 3D, the Matsubara self-energy, up to the quadratic in $`\epsilon _n`$ and $`T`$ terms, can be written as

$$i\mathrm{\Sigma }_0(i\epsilon _n,T)=\alpha \epsilon _n+i\beta v_F(pp_F)+\gamma \left[\left(\pi T\right)^2\epsilon _n^2\right].$$
(5)
In addition to a direct calculation, the validity of the quadratic term in (5) is readily established by noticing that upon analytic continuation $`i\epsilon _n\epsilon +i0`$ this term gives the correct form for the imaginary part of the on-shell self-energy: $`Im\mathrm{\Sigma }_0^R\left(\pi T\right)^2+\epsilon ^2`$ . The amplitude of the first harmonic (4) contains $`i\mathrm{\Sigma }_0(i\pi T,T)`$, in which the quadratic term vanishes identically. The $`T^2`$-terms from higher Matsubara frequencies (legitimately considered within this scheme in 3D) increase the amplitude and cannot be interpreted as ”damping.”

In 2D, the integral over momentum transfers diverges logarithmically at the lower limit, changing the behavior of $`Im\mathrm{\Sigma }_0^R`$ to $`E^2\mathrm{ln}E`$, where $`E=\text{max}\{\epsilon ,T\}`$. This change does not alter the principal result. Consider the simplest case of a contact interaction. To the second order in this interaction, the quadratic term in (5) is replaced by

$$i\stackrel{~}{\mathrm{\Sigma }}_0(i\epsilon _n,T)=\frac{U^2m}{\pi ^2v_F^2}T\underset{\omega _m=0}{\overset{\epsilon _n\pi T}{}}\omega _m\mathrm{ln}\left(\epsilon _F/\omega _m\right).$$
(6)
Although the sum in (6) does not have an analytic solution, it obviously vanishes for $`\epsilon _n=\pi T`$.

To analyze a general case of a finite-range and dynamic interaction, including the screened Coulomb one, it is convenient to find the imaginary part of the retarded self-energy first and then continue back to Matsubara frequencies. On the mass-shell, $`Im\mathrm{\Sigma }_0^R`$ is given by

$`Im\mathrm{\Sigma }_0^R\left(\epsilon \right)={\displaystyle \frac{\pi }{2\left(2\pi \right)^D}}{\displaystyle 𝑑\omega F\left(\omega \right)\left[\mathrm{coth}\frac{\omega }{2T}\mathrm{tanh}\frac{\omega \epsilon }{2T}\right]},`$

where $`F\left(\omega \right)=d^Dq\delta \left(\omega 𝐯_F𝐪\right)ImV^R(\omega ,q)`$ and $`V^R(\omega ,q)`$ is the retarded interaction potential. As a function of a complex variable $`z`$, $`f(z)Im\mathrm{\Sigma }^R\left(z\right)`$ has the following properties in the upper half-plane: i) all lines $`Imz=\pi \left(2n+1\right)T`$ are branch cuts on which $`Ref`$ is continuous but $`Imf`$ changes jumpwise; ii) due to the fact that $`\mathrm{tanh}(xi\pi (n+1/2))=\mathrm{coth}x`$, all points $`z=i\pi \left(2n+1\right)T`$ are zeroes of $`f(z).`$ Thus, function $`f(z)`$ is analytic in the band $`0Imz<\pi T`$ including the point $`z=i\pi T`$. Analytic continuation from the real axis into this band is legitimate and at $`z=i\pi T`$ it yields the Matsubara self-energy $`\stackrel{~}{\mathrm{\Sigma }}_0(i\epsilon _0=i\pi T,T)`$, which is equal to zero. Zeroes of $`f\left(z\right)`$ at $`z=i\pi \left(2n+1\right)T`$ with $`n1`$ do not lead to vanishing of $`\stackrel{~}{\mathrm{\Sigma }}_0(i\epsilon _n,T)`$ for $`n1`$ because those zeroes are separated from the real axis by branch cuts and thus are not accessible by analytic continuation. The 2D case is special only in that the $`q`$integration results in the $`\mathrm{ln}\omega `$ factor in $`F(\omega )`$ which does not change the reasoning given above. In particular, for a dynamically screened Coulomb interaction in 2D, $`F\left(\omega \right)\omega \mathrm{ln}|\omega |`$ and still $`\stackrel{~}{\mathrm{\Sigma }}_0(i\pi T,T)=0`$. As this result does not depend on the particular form of the interaction, it can be viewed as a generalization of the Fowler-Prange theorem. Notice that in 3D the Fowler-Prange theorem is of limited applicability because nothing prevents one from considering $`k>1`$ and lower values of $`T+T_D`$, when the effect of $`\mathrm{\Sigma }_0(i\epsilon _{n>0},T)`$ needs to be taken into account. In 2D, one is bound to consider only $`\mathrm{\Sigma }_0(i\epsilon _{n=0},T)`$ within the Luttinger approximation.

Finally, we discuss the effect of interference between electron-electron and electron-impurity scattering on magneto-oscillations, extending the analysis of the interference corrections to the self-energy in 2D from the diffusive ($`T\tau 1`$)  to the ballistic ($`T\tau 1`$) limit. The general form of the interference correction to the Matsubara self-energy is (see Fig. 1)

$`\mathrm{\Sigma }_0^{\text{int}}(i\epsilon _n,𝐩)`$ $`=`$ $`2T{\displaystyle \underset{\epsilon _n\left(\omega _m\epsilon _n\right)>0}{}}{\displaystyle \frac{d^2q}{\left(2\pi \right)^2}V(i\omega _m,q)}`$ (8)
$`\times \mathrm{\Gamma }(i\omega _m,q)G(i\epsilon _ni\omega _m,𝐩𝐪),`$
where $`G(i\epsilon _n,𝐩)=\left(i\epsilon _n\xi _p+i\text{sgn}\epsilon _n/2\tau \right)^1`$ and the effective interaction $`V=V_s+3V_t`$ contains the contributions from both singlet and triplet channels

$`V_s^1(i\omega _m,q)`$ $`=`$ $`(2\pi e^2/q+F_\rho ^0/\nu )^1\mathrm{\Pi }(i\omega _m,q);`$ (9)
$`V_\sigma ^1(i\omega _m,q)`$ $`=`$ $`\nu /F_\sigma ^0\mathrm{\Pi }(i\omega _m,q),`$ (10)
where $`\nu =m/\pi `$. The factor of two in (8) accounts for two possibilities of including the vertex correction in Fig.1a. The Fermi-liquid constants $`F_\rho ^0`$ and $`F_\sigma ^0`$  determine the renormalized charge and spin-susceptibilities, respectively. The vertex $`\mathrm{\Gamma }(i\omega _m,q)=(\sqrt{\left(\left|\omega _m\right|\tau +1\right)^2+\left(qv_F\tau \right)^2}1)^1`$ reduces to $`\mathrm{\Gamma }=\left(D\tau q^2+\left|\omega _m\right|\tau \right)^1`$ and $`\mathrm{\Gamma }=\tau ^1\left(\left|\omega _m\right|^2+\left(qv_F\right)^2\right)^{1/2}`$ in the diffusive and ballistic limits, respectively. The general form of the (small $`q`$) polarizarion operator

$`\mathrm{\Pi }(i\omega _m,q)=\nu \left[1\left|\omega _m\right|\tau \mathrm{\Gamma }(i\omega _m,q)\right]`$ (11)
reduces to $`\mathrm{\Pi }(i\omega _m,q)=\nu Dq^2/\left(Dq^2+\left|\omega _m\right|\right)`$ in the diffusive limit, where $`D=v_F^2\tau /2`$, and to $`\mathrm{\Pi }^0(i\omega _m,q)=\nu \left(1\left|\omega _m\right|/\sqrt{\left(v_Fq\right)^2+\omega _m^2}\right)`$ in the ballistic one. Omitting the details of lengthy but straightforward calculations, we give just the result for the self-energy valid to logarithmic accuracy

$`i\mathrm{\Sigma }_0^{\text{int}}(i\pi T,T)=T\mathrm{ln}\left(\epsilon _F/T\right)Q(T\tau ,F_\sigma ^0)/2\epsilon _F\tau ,`$

where $`Q(T\tau ,F_\sigma ^0)=g_\rho \left(T\tau \right)+\left[3F_\sigma ^0/\left(1+F_\sigma ^0\right)\right]g_\sigma \left(T\tau \right)`$ and $`g_{\rho /\sigma }\left(x\right)`$ are slowly varying functions which interpolate between the diffusive and ballistic regimes. The limiting values of $`g_{\rho /\sigma }\left(x\right)`$ are as follows: $`g_\rho \left(0\right)=1,`$ $`g_\rho \left(x1\right)=3/2,`$ $`g_\sigma \left(0\right)=1,`$ $`g_\sigma \left(x1\right)=1/2`$. Apart from the numerical coefficients, the $`T`$dependence of $`i\mathrm{\Sigma }_0^{\text{int}}(\pi T,T)`$ is the same in the diffusive and ballistic regimes. In that sense, the behavior of the self-energy is sumilar to that of the tunneling density of states . Notice that the interference correction to the scattering rate in the ballistic regime $`\left|Im\mathrm{\Sigma }_0^R\right|_{\text{int}}\left(T/\epsilon _F\tau \right)\mathrm{ln}\left(\epsilon _F/T\right)`$ is smaller than the scattering rate in a clean FL, $`\left|Im\mathrm{\Sigma }_0^R\right|_{\text{int}}\left(T^2/\epsilon _F\right)\mathrm{ln}\left(\epsilon _F/T\right)`$, in the parameter $`\left(T\tau \right)^11`$. However, due to the cancellation of the $`T^2\mathrm{ln}T`$ term in $`i\mathrm{\Sigma }_0(i\pi T,T)`$ discussed earlier in this paper, the interference correction is the main nonlinear $`T`$dependent term in the Matsubara self-energy, leading to a modification of the LK-formula. The $`T\mathrm{ln}T`$dependence of the self-energy can be interpreted as a logarithmic $`T`$dependence of the effective mass. Following this interpretation, the interference effect leads to a replacement of the effective mass in the argument of the exponential in $`A_1`$ by

$$m^{}\left(T\right)=m^{}\left(1\frac{m}{m^{}}\frac{\mathrm{ln}\left(\epsilon _F/T\right)}{2\pi \epsilon _F\tau }Q(T\tau ,F_\sigma ^0)\right).$$
(12)
Notice that the effective mass is reduced by the singlet-channel interaction but enhanced by the ferromagnetic $`\left(F_\sigma ^0<0\right)`$ interaction in the triplet channel. The $`\mathrm{ln}T`$dependence of the effective mass is a characteristic feature of the “marginal Fermi liquid” model .

Equivalently, the nonlinear $`T`$dependence of $`i\mathrm{\Sigma }_0(i\pi T,T)`$ may be interpreted as a $`T`$dependent Dingle temperature:

$$T_D\left(T\right)=T_D\left(1\left(T/\epsilon _F\right)\mathrm{ln}\left(\epsilon _F/T\right)Q(T\tau ,F_\sigma ^0)\right).$$
(13)
One of the empirical procedures used in  to account for the observed deviations from the LK formula was to assume that the effective Dingle temperature has the same $`T`$dependence as the zero-field resistivity. Comparing (13) with the result for the interference correction to the resistivity , we see that, although this recipe is not precise, it has some theoretical justification: the general structure of the $`T_D\left(T\right)`$ and $`\rho \left(T\right)`$ depedences is similar in the ballistic regime, except for a factor of $`\mathrm{ln}\left(\epsilon _F/T\right)`$ present in $`T_D\left(T\right)`$ but not in $`\rho \left(T\right)`$.

This work was supported by NSF DMR-0077825. We are grateful to I. L. Aleiner, A. V. Chubukov, J. Klauder, and B. Z. Spivak for valuable comments and to M. E. Gershenson and V. M. Pudalov for extensive discussions of the experimental data.