text:
# Nanofriction behavior of cluster-assembled carbon films

## I INTRODUCTION

The understanding and control of friction, adhesion, lubrication and wear in nanostructured systems is an essential requisite to validate the use of nanomaterials for highly demanding structural applications Bhushan (1995); Carpick and Salmeron (1997). The rapidly growing number of applications based on microelectromechanical systems (MEMS) and the new perspectives opened by the production of nano electromechanical systems (NEMS) make necessary the development of an entirely novel class of protective and lubricant coatings with improved mechanical properties at the nanoscale mrs (2001). Carbon-based materials have emerged as a promising class of materials Sullivan et al. (2001). Diamond and the various types of amorphous carbon show interesting tribological properties such as high elastic moduli, good lubrication properties, low stiction surfaces, etc. Sullivan et al. (2001). These properties are controlled by a combination of physico-chemical properties of the surfaces such as sp2/sp3 ratio, surface roughness and porosity and by the presence of contaminant layers Sundararajan and Bhushan (1999).

Cluster-assembled carbon represents a novel nanostructured material obtained by the deposition of clusters produced in a supersonic expansion Milani et al. (2001). Cluster-assembled carbon films can be obtained starting from cluster mass distributions that contain clusters with fullerene-like structure and/or linear and planar structures Donadio et al. (1999). Low-energy deposition (fractions of eV per atom) drastically reduces fragmentation on the substrate allowing the films to be structured at the nano- and mesoscale by keeping memory of the original cluster distribution Barborini et al. (1999a).The mechanical properties of cluster-assembled carbon films have been studied by Brillouin light scattering showing that these systems have elastic properties similar to graphite (shear modulus) and Young modulus and Poisson’s ratio typical of a soft very porous material (note that the Poisson’s ratio is very close to zero or, in certain cases, even negative) Casari et al. (2001).

The Atomic Force-Friction Force Microscope (AFM-FFM) has emerged as a powerful tool for the characterization of the tribological properties of materials from the micrometer down to the atomic scale Carpick and Salmeron (1997); Dedkov (2000). However, the use of a nanometer-sized probe in nano-friction experiments carried out on corrugated samples in humid environments causes the experimental conditions to be usually different from both those encountered in typical ultra high-vacuum experiments, where flat crystalline surfaces are investigated in a humidity and contaminants-free environment, leading to a single-asperity contact, and those of macroscopic tribology, where the contact regime is always multi-asperity like Adams and Nosonovsky (2000). A number of parameters such as adhesion, surface and tip micro and nano-roughness, load range, as well as tip radius and shape, influence the contact-friction regime. Moreover, local corrugation on a scale larger than that of the probe can affect the friction measurement Koinkar and Bhushan (1997); Sundararajan and Bushan (2000). The case of nanostructured materials is somehow peculiar since the typical size of the probe is close to or even comparable with cluster size. In this case peculiar tip-sample interactions and topographic effects should be expected.

In this paper we present the results of an AFM-FFM characterization of the frictional properties of cluster-assembled carbon films. We have studied the dependence of frictional parameters on both the relative humidity and the structural composition of films deposited with different cluster mass distributions.

The experimental data are discussed on the basis of a modified Amonton’s law for friction Bowden and Tabor (1950), stating a linear dependence of friction on load plus an adhesive offset accounting for a finite friction force in the limit of null total applied load. A new procedure for the correction of the lateral force maps from the contributions of the local tilt of the surface (the topographic correction) was applied in order to extract intrinsic values of the frictional parameters independent on surface roughness on a scale larger than that of the tip-sample contact. In order to validate the use of the Amonton’s law for the interpretation of the experimental data, we have simulated the AFM tip-nanostructured carbon interaction via molecular dynamics (MD) simulations. MD results of the nanofriction experiment support the validity of the friction model used for this system.

## II EXPERIMENTAL DETAILS

### II.1 Deposition of nanostructured carbon

The use of supersonic beams of clusters for deposition of thin films has attracted a large interest from more than two decadesMilani and Iannotta (1999). This technique consists in preparing clusters in the gas phase diluted in a carrier light inert gas, typically Helium, and letting the mixture expand through a nozzle in high vacuum so that a very collimated, intense supersonic beam is produced. The cluster beam is intercepted by a suitable substrate in order to deposit a thin film. In the case of carbon, the kinetic energy per atom in the clusters (below 0.4 eV/atom) is smaller than the binding energy per atom avoiding massive fragmentation: the resulting film has a structure at the nanoscale keeping the memory of the nanometer-sized building blocks used for the assembling.

Nanostructured carbon films have been deposited from a supersonic cluster beam produced by a pulsed microplasma cluster source (PMCS) as described in detail in Ref. Barborini et al., 1999b. With normal PMCS operation conditions, the cluster beam is characterized by a log-normal cluster mass distribution peaked at about 500 atoms/cluster and extending to several thousands atoms per cluster. We have controlled and varied the cluster mass distribution and deposition rates by exploiting aerodynamic focusing effects Piseri et al. (2001). Using the standard cluster mass distribution described above, we have deposited films with thicknesses of several hundreds of nanometers on silicon substrates. The deposition rate was 4-5 nm/min and the density of the films was 0.8-0.9 g/cm3. With a suitable nozzle configuration we have produced a beam depleted from clusters with diameters roughly larger than 2 nm Piseri et al. (2001). These conditions produce films with densities of 1.2-1.3 g/cm3 at a deposition rate of 5 nm/sec.

Raman analysis have shown that in films grown with large clusters the graphitic $`sp2`$ bonding, due to the large number of cage-like particles, is more pronounced than in films grown with small clusters, whose structure resembles more that of the amorphous carbon. In the following, we shall refer to films grown with the two different nozzles as to the films grown with small and large clusters, accordingly.

### II.2 AFM-FFM

The atomic force microscope is a Nanoscope Multimode IIIa from Digital Instruments with phase extender and Signal Acquisition Module (SAM). We scanned different points of the samples, with scan size 500 nm, sliding velocity typically 1 $`\mu m/s`$. We used rectangular cantilevers 450 $`\mu m`$ long with silicon tip with radius 5-40 nm. The AFM can be housed in a sealed chamber connected to a humidifier to work in controlled humidity and atmosphere (typically in dry and wet nitrogen). A second PC is used to remotely control the external applied load in friction measurements. The application of the external load is synchronized via the reading of the end-line and end-frame triggers from the microscope. We can acquire a complete friction vs. load curve in a single AFM scan of typically 512 lines, 512 points per line. This allows performing many measurements in each session, fast and reliably. The experimental data are extracted in ASCII format and processed via dedicated software.

## III MOLECULAR DYNAMICS SIMULATIONS

The theoretical investigation has been carried out by classical MD simulations based on the Tersoff potential Tersoff (1988a). This computational framework has proved to be reliable for the investigation of the structural properties of carbon based materials Tersoff (1988b) and particularly of nanostructured cluster assembled carbon films. The growth of nanostructured carbon films by supersonic cluster beam deposition on a (001) diamond substrate has been simulated by a protocol described in Ref. Donadio et al., 1999 which accounts well for the experimental deposition process. The simulation cell is a slab with periodic boundary conditions applied in the plane orthogonal to the growth direction, the substrate is four layers thick and its bottom layer is rigid, while the second and third layers are thermostat. The dynamics of the atoms of the forth layer is newtonian as that of the impinging clusters. Two samples have been produced with the same growth conditions (cluster kinetic energy, vibrational temperature, substrate temperature) but different size distribution of the precursors: sample (A) has been obtained from a cluster beam mostly containing small precursors (1 to 23 atoms per cluster) and sample (B) has been grown mainly from larger precursors (46 to 120 atoms per cluster). The structural properties of the two sample, which are strongly influenced by the features of the precursors, are discussed elsewhere Donadio et al. (1999).

Fig. 1 illustrates the molecular dynamics model for the FFM experiment: a small crystalline diamond tip slides on the surface of the nanostructured carbon. The tip consists of 178 carbon atoms arranged in a truncated pyramid shape, obtained by cutting an ideal diamond structure crystal along (111) and (001) directions and letting it relax at room temperature. Tip radius is about 8 Å. Although silicon tips are used in experimental measurements we have chosen to perform nanofriction simulations with a stiffer carbon tip, in order to avoid the complete wearing of the tip even at relatively low normal loads, due to its small size. Thus it has been possible to perform FFM simulations with loads up to 30 nN.

Since the thickness of the samples exploited for nano-friction simulations is about 60 Å, in order to reduce computational costs, a 20 Å thick layer of the film is fixed, the atoms belonging to the 2 Å thick layer above it are thermostat and the dynamics of the other atoms of the film is Newtonian. The two top layers of the tip are rigid and they are subject to the external forces that mimic the action of the AFM cantilever. The overall external force acting on the tip has been decomposed into three independent components: (i) an elastic dragging force that produces the sliding motion of the tip, (ii) a constant load perpendicular to the sliding plane and (iii) an elastic torsion force that keeps the tip vertical with respect to the sliding plane. The frictional force is calculated at every timestep as the reaction to the elastic dragging force and it is averaged over six paths $`20`$ Å long.

Since the formation of chemical bonds between the atoms of the AFM tip and those of the film is quite unlikely in such AFM-FFM experiments, the film–tip interaction is ruled by Van der Waals forces. Therefore the film–tip interaction is modelled by a two-body modified Morse potential given by:

$$V(r_{ij})=D(e^{2\alpha (r_{ij}dr_0)}2e^{\alpha (r_{ij}r_0)}),$$
(1)
where $`r_{ij}`$ is the distance between two atoms. As for $`D`$ and $`\alpha `$ parameters we have chosen the values fitted by Cheong et al. Cheong and Zhang (2000): namely $`D=0.435`$ eV and $`\alpha =46.487`$ . The parameter $`d`$ has been introduced in order to increase the equilibrium distance and to reduce the binding energy between the tip and the film. This reproduces better the actual experimental conditions: as the FFM measurements are usually performed in air, we have to assume that the surface-tip interaction is modified by hydrogen passivation of the dangling bonds and by the presence of surface lubricants or contaminants. Hence, in principle, by tuning this parameter one may study the effects of lubricants and contaminants without explicitly introduce them. We have chosen $`d=0.011`$ nm, by performing sliding friction simulations on the (001) non reconstructed surface of diamond and comparing the calculated lateral force to those measured in Ref. Mate et al., 1987.

## IV Results and discussion

### IV.1 AFM-FFM measurements

Our characterization protocol assumes a modified Amonton’s law for friction, i.e. a linear dependence on total load plus an offset representing the zero-load friction force in analogy with the single-asperity JKR contact modelJohnson et al. (1971):

$$f=\mu N+c$$
(2)
Here $`N`$ represents the total applied load in the direction perpendicular to the surface, including the contribution of adhesion $`A`$, and $`c`$ accounts for a zero-load friction force.

We apply a procedure, described in details elsewhere Podestà (2002); Podestà et al. (2002), to correct lateral force maps from the spurious contributions due to the presence of a local tilt of the surface. Actually, in the case of a locally tilted surface, the measured forces in the directions parallel and perpendicular to the AFM reference plane do not necessarily coincide with the forces acting parallel and perpendicularly to the sample surface, which actually define the friction coefficient and the friction vs. load characteristics of the interface under investigation. The topographic correction is necessary in order to extract quantitative and accurate information from corrugated sample and to compare results from different samples.

We have investigated frictional properties of films grown with large and small clusters in dry and humid nitrogen environment.

In Fig. 2 are shown two AFM pictures at different magnifications of a 250 nm thick ns-carbon film. Typical film morphology consists in a fine raster of grains with typical diameter of 10-20 nm. Films grown with large clusters are usually rougher than films grown with small clusters with the same thickness. However, our protocol automatically corrects the friction maps for the topographic contributions, allowing the comparison of the results obtained on different samples.

In Fig. 3 is shown a typical lateral force vs. applied load curve measured on ns-carbon. This curve is obtained extracting from the original lateral force-load dispersion the subset of all the lateral force-load pairs corresponding to the same slope in the topographic map. Each of such curves is processed separately in order to apply the topographic correction.

We show in Fig. 4 the measured dynamic friction coefficients and friction offsets for the film grown with large clusters and grown with small clusters accordingly, measured in ambient condition (RH$`40\%`$)111One should keep in mind that we have measured the friction coefficient and adhesive offset of the pair (ns-carbon, Silicon AFM tip). Using AFM tips with different composition can lead in principle to different friction coefficients and adhesive offsets..

The average values for both $`\mu `$ and $`c`$ are:

$`\mu _{large}^{humid}=0.086\pm 0.013,`$ $`c_{large}^{humid}=0.44\pm 0.29nN`$
$`\mu _{small}^{humid}=0.085\pm 0.007,`$ $`c_{small}^{humid}=0.26\pm 0.15nN`$ (3)
We also studied the frictional behavior of ns-carbon film in dry environment (RH$`2\%`$). The results are shown in Fig. 5. The averaged values are:

$`\mu _{large}^{dry}=0.081\pm 0.012,`$ $`c_{large}^{dry}=0.69\pm 0.28nN`$
$`\mu _{small}^{dry}=0.082\pm 0.008,`$ $`c_{small}^{dry}=0.47\pm 0.20nN`$ (4)
The errors shown in Eqs. 3 and 4 are calculated as: $`\sigma =1/N\sqrt{\sigma _i^2}`$, where $`\sigma _i`$ are the errors associated to each measurements.

### IV.2 MD results

Exploiting the model setup described in Section III, we have performed several nanofriction simulations, calculating the values of lateral friction force as a function of the applied normal load, both for sample (A) and for sample (B). Each simulated measurement at a fixed load, consists of six independent scannings 25 Å long. The friction force for each scansion is obtained by the average of the lateral force calculated at every time-step of the simulation and the values reported in Fig. 6 are calculated by averaging over the six scannings. This procedure provides a sufficient amount of data so to devise a linear trend in the lateral force as a function of applied load and to calculate the slope, i.e. the friction coefficient, and the offset of the curves according to the modified Amonton’s law discussed above.

The friction coefficient and offset values thus obtained are:

$`\mu ^A=0.47\pm 0.110,`$ $`c^A=1.81\pm 1.43,nN`$
$`\mu ^B=0.42\pm 0.053,`$ $`c^B=0.70\pm 0.59,nN`$ (5)
Simulation have been performed with loads up to 22.6 nN for sample (A) and 19.2 nN for sample (B). Beyond these threshold values, combined wearing of the tip and of the film occurs and lateral force rapidly increases above the linear behavior observed so far.

No topographic correction to the theoretically calculated friction vs. load curves has been applied, since the method used to calculate each single point of the curves averages over the scanned profile and this should already account for the topographic effects. Moreover, the average tilt angle of the MD profiles is in general smaller than 15 degrees, making the topographic correction negligible.

### IV.3 Discussion

The observed lateral force vs. applied load curves for ns-carbon are linear, except for very low loads, close to the pull-off limit, and for loads larger than about 30 nN (Fig. 3). Moreover, in the limit of the experimental error, the measured values of the friction offset are definitely different from zero. These observations confirm that the friction law of a silicon tip on cluster-assembled carbon is well represented by an Amonton’s like equation, and in addition that a zero-load offset must be included in the friction model.

MD results further support our model based on a modified Amonton’s law. Simulations show a linear trend of the friction vs. load curve with a finite offset, different from zero in the limit of the experimental error. The quantitative discrepancy between theoretical and experimental results is not surprising and it can be attributed to the parameterization of the interaction potentials and in the choice of using carbon tips in the simulations of the FFM experiments.

For loads larger than about 30 nN we have evidence of the onset of non-linear trends in friction vs. load curves, probably related to wear and indentation during scanning. The phenomenon of wearing of the ns-carbon film has been directly observed in MD simulations at loads larger than 30 nN.

It has been reported that the values of the friction coefficients measured in successive scans are likely to be different Carpick and Salmeron (1997); Koinkar and Bhushan (1996). This is currently attributed to modification of the tip-sample interface after repeated sliding due to material removal, tip wear and blunting and even structural modifications. Blunting of the tip is a very unlike event in this case because Silicon is much harder than cluster-assembled carbon. The observed fluctuations in the measured parameters (see Figs. 4 and 5) could be attributed to dynamic modifications of the contact area during sliding because of the contamination of the AFM tip, which may pick up carbon during scanning, particularly at higher loads, and even indent the film, which is much softer than the tip.

The strongest fluctuations are observed for the experimental friction offset $`c`$. This is not surprising, because this parameter is expected to be highly dependent not only upon environmental conditions (such as relative humidity and temperature), but also upon tip shape and contamination, these being the parameters directly influencing adhesion. This fact is clearly represented by the large errors in the measured adhesive offset. While the relative error for the friction coefficient are below 15$`\%`$, which is basically the limit set by the accuracy of the calibration procedure for the force constants of the cantileverGibson et al. (1996)222The statistical error in each measurement is in general very small, because with our procedure we may acquire up to 512x512 friction-load pairs in each image and many images fast and reliably., those of the adhesive offset are as large as 60$`\%`$. In the absence of an accurate characterization of tip shape and size however it is not possible to definitely individuate the causes of the observed changes in the friction offset $`c`$. The investigation of adhesive properties of cluster-assembled carbon should be more reliable through the analysis of force vs. distance curves, or even via the study of thermal fluctuations of the cantilever close to the nanostructured surface, provided a reliable model of the tip-sample interaction and/or contact.

As confirmed also by MD simulations, in the limit of the experimental error, the nanostructure of the films seems to have very weak or no influence on the friction coefficient. In analogy to the macroscale, one could expect that films grown with large clusters and hence containing a larger number of well-organized graphitic regions should show a lower friction coefficient Grill (1997). This is not observed.

The friction coefficient is also found to be independent on relative humidity. Null or weak dependence of friction coefficient and adhesion on the relative humidity was reported in literature for carbon-based materialsBinggeli and Mate (1994); Koinkar and Bhushan (1996); Bhushan and Dandavate (2000). One possible explanation for the behavior of our films in dry environment is the hydrophilic nature of cluster-assembled carbon Ostrovskaya et al. (2002). The low density and large porosity of these systems in conjunction with a non-negligible wettability may favor the formation of a water layer lasting even in dry nitrogen with the formation of a water meniscus between the tip and the surface

Friction coefficients measured on carbon-based materials span from 0.009 (for HOPG in air) to 1 for DLC in UHV Grill (1997). A numerical comparison with nanostructured carbon would not be of great interest because of the large scattering of the results and the arbitrariness of the calibration procedures adopted. What is in general observed is that all carbon-based materials are better lubricants in the presence of surface contamination, while in very clean and dry environments the covalent interaction between sample and probe (usually Silicon or Diamond) leads to higher friction coefficients. In the case of cluster-assembled carbon we do not observe an increase of the friction coefficient upon removal of the environmental humidity.

The adhesive offset $`c`$ is found experimentally to decrease in humid environment. The decrease of the adhesion force with increasing relative humidity was observed in other systems Binggeli and Mate (1994). This is explained by a decrease of the capillary attractive force exerted by the water meniscus in the tip-surface gap, with higher relative humidity, upon a certain threshold, depending on the material. This could explain also the larger adhesion coefficient observed in the films grown from large clusters. This system with a more pronounced graphitic character may have a smaller wettability compared to that grown by small clusters Ostrovskaya et al. (2002).

In the simulations the effect of humidity is poorly reproduced by the introduction of the parameter $`d`$ in the empirical potential describing the tip-film interaction since it is impossible to establish a clear relation between the relative humidity and the value of $`d`$. Anyway, we have observed that in completely dry conditions, corresponding to $`d=0`$ Å, strong adhesion of the tip on the surface occurs, resulting in the complete wearing of the tip and in large deformation of the scanned surface, even for loads as small as 5 nN.

Recently Buzio et al. reported a friction characterization on cluster-assembled carbon films produced by our technique R. Buzio and Valbusa (2002). They have found a non-linear dependence of friction on load, well fitted by the Hertzian-plus-offset model Schwarz et al. (1997), using a silicon tip of $`30`$ nm and scan length of 50 nm. On the other hand, a linear dependence of friction on load is found by scanning on a larger area (1x1 $`\mu `$m2) Gnecco et al. (2000). The authors attributed the non-linear dependence to a possible passivation of the sample surface and to the presence of a water layer, which would induce a transition from a multi-asperity to a single-asperity contact (i.e. from a linear to a power-law dependence) as explained in the framework of the composite-tip model Putman and Igarashi (1995). For larger scan length and relatively larger scan velocity, the tip-sample interface is broken and re-formed rapidly and the smoothing action of water and impurities is less effective. This hypothesis is confirmed by our measurements apart from differences observed in the numerical values of the friction coefficients.

## V CONCLUSIONS

We have characterized by FFM-AFM the tribological properties of nanostructured carbon films grown by supersonic cluster beam deposition. We have found that the friction behavior at the nanoscale can be described by a modified Amonton’s law demonstrating that a nanometer-sized contact not necessarily leads to a single-asperity friction regime. Molecular Dynamics simulations support our interpretation and show a qualitative agreement with the experiments. Our results show that the presence of water and the high reactivity of $`sp2`$-$`sp3`$ material plays an important role in determining a multi-asperity like contact also at these very small scales. Although numerical comparison with other carbon-based systems must be considered with high caution, cluster-assembled carbon shows a low friction coefficient even in a dry environment. The nanoscale structure of the films seem to have a negligible influence on the friction at the observed scales. This suggests that the tribological behavior of a nanostructured solid is ”scale-sensitive” and one should always consider the scale factor before making comparison or extrapolations.

###### Acknowledgements.

Partial financial support from Italian Space Agency is acknowledged.

## References

## Figure captions

1. Representation of the model used for AFM-FFM molecular dynamics simulations. Only a $``$50 Å thick top layer of nanostructured carbon film is shown. Color map reflects the coordination of atoms in the ns-carbon film: Blue, 1-fold coordinated; yellow, 2; green, 3; red, 4.
2. AFM pictures of a 250 nm thick ns-carbon film at different magnifications. The scan size is 2 $`\mu `$m (top) and 500 nm (bottom). The vertical color scale is 100 nm. The higher magnification picture clearly shows the surface granularity of such nanostructured films.
3. Lateral force vs. (external) applied load curve measured on ns-carbon.
4. Measured friction coefficients $`\mu `$ and friction offset $`c`$ for films grown with large clusters and grown with small clusters in humid nitrogen (RH$`40\%`$).
5. Measured friction coefficients $`\mu `$ and friction offset $`c`$ for films grown with large clusters and grown with small clusters in dry nitrogen (RH$`2\%`$).
6. Lateral force vs. applied load values and fitting from MD simulations, for sample (A) (circles and solid line) and for sample (B) (squares and dashed line). Error bars represent the standard deviation of the lateral force calculated for each scansion at a given load.


format: text

text:
# Nonlinear screening in two-dimensional electron gases

## Abstract

We have performed self-consistent calculations of the nonlinear screening of a point charge $`Z`$ in a two-dimensional electron gas using a density functional theory method. We find that the screened potential for a $`Z=1`$ charge supports a bound state even in the high density limit where one might expect perturbation theory to apply. To explain this behaviour, we prove a theorem to show that the results of linear response theory are in fact correct even though bound states exist.

PACS Numbers: 73.20.-r, 73.20.Hb

Screening is a fundamental property of an electron gas in arbitrary dimensions. The example of two dimensions is of particular interest because of the possible realization of quasi-two-dimensional systems in a variety of contexts: semiconductor heterostructures, image or band-gap surface states at metal surfaces, electrons on the surface of liquid helium, and layered materials. In all of these cases, the interaction of external charges with the two-dimensional electron gas (2DEG) is a problem of both fundamental and practical interest. For example, the transport of electrons in a 2DEG is often limited by charged impurity scattering and a detailed knowledge of the scattering potential is needed for an accurate determination of the electron mobility. Scanning tunnelling microscopy offers an even more direct means of determining the screening response of a quasi-2DEG through the observation of adsorbate-induced Friedel oscillations. Still another class of problems involves the interaction with moving charges as might arise in low-energy electron scattering or tunnelling experiments. In this case, the dynamic response of the 2DEG is important in that electronic excitation, and hence energy loss, will occur.

A charged impurity or projectile typically represents a strong perturbation and a nonlinear screening theory is in general needed to account for the modifications of the local electronic structure. However, in certain situations the screened impurity potential may be relatively weak and therefore amenable to a perturbative treatment. This is usually the method adopted to deal with donor impurities that are spatially removed from the 2DEG within a heterostructure, although it is rare to find quantitative agreement between theory and experimentally measured mobilities. The situation of acceptor impurities within the gas is a much more severe perturbation and quite dramatic effects can arise as a result of the modified electronic structure. In such situations, the screening response has to be determined nonlinearly.

One of our objectives in this Letter is to provide a fully self-consistent description of the nonlinear screening in an ideal 2DEG within the context of density functional theory. A second objective is to use these calculations to establish the range of validity of linear response theory. Somewhat surprisingly, this latter objective is more subtle than anticipated as a result of a peculiarity of potential scattering in 2D, namely the fact that any purely attractive potential always has at least one bound state, in marked contrast to the situation in 3D. For a positively charged impurity we in fact encounter a situation in which the screened potential supports a bound state even in the high density limit where one would intuitively expect a perturbative treatment to be valid. To resolve this apparent paradox, we prove what will be referred to as the high density screening theorem which states that the screening charge can indeed be calculated perturbatively even when bound states exist. In this way we are able to justify the use of perturbation theory when at first sight it would seem inapplicable.

The problem we address is the nonlinear screening of a stationary point charge, $`Z`$, located in the plane of a 2DEG. The latter is treated as ideal in the sense that the electrons are confined to the plane. Of course in real applications, such as a heterostructure, a more accurate treatment of the electronic states is required. We ignore these complications in order to focus on those aspects of nonlinear screening which are expected to be independent of these details. To stabilize the system, the electrons move in the presence of a uniform neutralizing positive background. In addition, we assume the electrons to have an isotropic effective mass $`m^{}`$ and to be immersed in an extended dielectric with permittivity $`\epsilon `$. We use the effective Bohr radius, $`a_0=\epsilon \mathrm{}^2/m^{}e^2`$, as the unit of length and the effective Hartree, $`H=e^2/\epsilon a_0`$, as the unit of energy. The density of the gas, $`n_0`$, is characterized by the density parameter $`r_s=1/\sqrt{\pi n_0}`$.

The static screening response of the 2DEG is determined by solving self-consistently the two-dimensional Kohn-Sham equations

$$\frac{1}{2}^2\psi _i(𝐫)+\mathrm{\Delta }v_{\mathrm{eff}}(𝐫)\psi _i(𝐫)=E_i\psi _i(𝐫)$$
(1)
where the effective potential is given by

$$\mathrm{\Delta }v_{\mathrm{eff}}(𝐫)=v_{\mathrm{ext}}(𝐫)+\mathrm{\Delta }v_H(𝐫)+\mathrm{\Delta }v_{xc}(𝐫).$$
(2)
Here, $`v_{\mathrm{ext}}=Z/r`$ is the external potential and $`\mathrm{\Delta }v_H`$ is the Hartree potential due to the electronic screening density, $`\mathrm{\Delta }n(𝐫)=n(𝐫)n_0`$. The change in the exchange-correlation potential, $`\mathrm{\Delta }v_{xc}(𝐫)=v_{xc}[n(𝐫)]v_{xc}[n_0]`$, is defined in the local density approximation (LDA) using the parameterization of the 2D exchange-correlation energy given in Ref. .

The total screening charge is given by

$$\mathrm{\Delta }n(𝐫)=\underset{b}{}|\psi _b(𝐫)|^2+\underset{i}{}\left[|\psi _i(𝐫)|^2|\psi _i^0(𝐫)|^2\right]$$
(3)
where the first sum extends over all bound states of the effective potential, and the second extends over all occupied continuum states up to the Fermi level $`E_F`$. We assume that each spatial orbital is doubly-occupied for spin. The scattering states $`\psi _i(𝐫)`$ are shifted asymptotically relative to the free-particle solutions $`\psi _i^0(𝐫)`$ by a phase $`\eta _m(E)`$. These scattering phase shifts are related to the total screening charge according to the 2D Friedel sum rule (FSR)

$$Z_{FSR}=\frac{2}{\pi }\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\eta _m(E_F).$$
(4)
Details of the self-consistent solution of Eqs. (1) and (2) will be presented elsewhere.

We begin by considering the case of a negatively charged impurity ($`Z=1`$), such as an antiproton or acceptor state in a semiconductor. Fig. 1 presents the self-consistent effective potential $`\mathrm{\Delta }v_{\mathrm{eff}}`$ for $`r_s=4`$ (solid curve) as a function of the distance from the impurity. This potential repels electrons almost completely from the impurity’s vicinity, leaving exposed a positively charged disc of radius $`Rr_s`$ which neutralizes the impurity charge. As is the case in 3D, this behaviour cannot be reproduced in a linear theory. The screened $`Z=1`$ potential does not support bound states for any $`r_s`$ value, as might have been expected. We mention this since it has been claimed that the introduction of a negative test charge into a 2D gas can give rise to a potential which is sufficiently strong to bind an electron.

To make contact with this earlier work we compare in Fig. 1 the nonlinearly screened potentials with those obtained on the basis of linear response theory. The chain curve shows the linear response Hartree potential ($`v_{\mathrm{ext}}+\mathrm{\Delta }v_H`$) as obtained when local field corrections (LFC) are included in the determination of the electron screening density. This potential has a large attractive region in real space and supports a bound state for a unit negative test charge of one electron mass. This observation led to the suggestion of a possible pairing mechanism that could be responsible for a correlation-induced instability at low densities. Such a conclusion, however, is invalid on two counts. First, the screening of the impurity is strongly nonlinear. The dashed curve in Fig. 1 shows the corresponding Hartree potential when the nonlinear screening density is used. It has a much shallower attractive region. But more importantly, an electron, as opposed to a negative test charge, also feels the effect of the induced xc potential. With this contribution included in the full nonlinear potential $`\mathrm{\Delta }v_{\mathrm{eff}}`$ (solid curve), there is no tendency for bound state formation, as confirmed numerically.

The results for a positive impurity are quite different in that the attractive screened potential supports bound states for all densities of physical interest (including $`r_s0`$). For $`Z=1`$ there is one $`m=0`$ bound state that is doubly occupied. Since the total screening density integrates to unity to satisfy the FSR, the continuum screening density must itself contribute a total charge of +1 in order to compensate for the overscreening provided by the bound states. In other words, the $`Z=1`$ impurity with 2 bound electrons can be viewed as an H<sup>-</sup> ion which acts as a $`Z=1`$ impurity. This was confirmed by comparing the $`Z=1`$ and $`Z=1`$ continuum screening densities. For $`r_s=10`$, the two are virtually the same.

In Fig. 2 we show the $`Z=1`$ bound state energy as a function of $`r_s`$. The behaviour seen is surprising in view of the corresponding behaviour in 3D. There, the bound state energy increases with decreasing $`r_s`$ since the impurity potential is screened more effectively with increasing density and, as a result, the bound state eventually ceases to exist. Beyond this point, the accuracy of perturbation theory improves with increasing density. The contrary behaviour exhibited in Fig. 2 calls into question the applicability of perturbation theory in the 2D case.

To address this question, we now prove a theorem regarding electronic screening in the high density limit. We consider the introduction of an external potential $`\lambda V(𝐫)`$ into a uniform noninteracting Fermi gas for arbitrary dimension $`D`$. The parameter $`\lambda `$ is a coupling constant whose physical value is unity. The problem is to determine the induced density $`\mathrm{\Delta }n(𝐫)`$ due to the introduction of $`\lambda V(𝐫)`$; the effect of interactions will be dealt with afterwards. By making use of Dyson’s equation, $`G_{\lambda +\delta \lambda }=G_\lambda +\delta \lambda G_\lambda VG_{\lambda +\delta \lambda }`$, for the single particle Green function $`G_\lambda (𝐫,𝐫^{},z)`$, one can show that the screening density satisfies

$`{\displaystyle \frac{\mathrm{\Delta }n(𝐫;\lambda )}{\lambda }}=`$ $``$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle 𝑑𝐫^{}V(𝐫^{})_{\mathrm{}}^{E_F}𝑑E}`$ (5)
$`\times `$ $`\text{Im}\left[G_\lambda (𝐫,𝐫^{},E+iϵ)G_\lambda (𝐫^{},𝐫,E+iϵ)\right].`$ (6)
The crucial next step in the argument is to use the analytic properties of the Green functions to change the energy integral from $`(_{\mathrm{}}^{E_F})`$ to $`(_{E_F}^{\mathrm{}})`$. In the $`E_F\mathrm{}`$ limit it is then permissible to replace $`G_\lambda `$ by the free-particle Green function $`G_0`$ since the energy $`E`$ can now be assumed to be much larger than the strength of the (bounded) potential $`|V(𝐫)|`$. Once this is done, the energy integration can be changed back to its original range, and after integrating with respect to the coupling constant, we obtain

$`\mathrm{\Delta }n(𝐫)`$ $``$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle 𝑑𝐫^{}V(𝐫^{})\underset{E_F\mathrm{}}{lim}_{\mathrm{}}^{E_F}𝑑E}`$ (7)
$`\times `$ $`\text{Im}\left[G_0(𝐫,𝐫^{},E+iϵ)G_0(𝐫^{},𝐫,E+iϵ)\right].`$ (8)
This asymptotic result is valid for any $`D`$ and applies even when the potential $`V(𝐫)`$ supports bound states.

An alternative expression for (6) is

$$\mathrm{\Delta }n(𝐫)\underset{E_F\mathrm{}}{lim}\frac{1}{L^D}\underset{𝐪}{}\chi _0(𝐪)V(𝐪)e^{i𝐪𝐫}$$
(9)
which is the result of linear response theory. Here, $`\chi _0(𝐪)`$ is the noninteracting static density response function of the system. In the high density limit, $`qk_F`$ for all wavevectors for which $`V(𝐪)`$ is finite, and we have $`\mathrm{\Delta }n(𝐫)\chi _0(0)V(𝐫)`$. In 2D, $`\chi _0(0)=\pi ^1`$ and

$$\mathrm{\Delta }n(𝐫)\frac{1}{\pi }V(𝐫).$$
(10)
Thus we arrive at the interesting conclusion that in 2D the screening density takes on a density-independent form in the high density limit, and is simply proportional to the perturbing potential.

This result can easily be checked numerically. In Fig. 3 we show $`\mathrm{\Delta }n(𝐫)`$ for a 2D gas with $`r_s=0.5`$ for a model potential $`V(r)=V_0\mathrm{sin}^2(2\pi r/r_0)\theta (r_0r)`$ which is an axially symmetric, double-well potential. With $`r_0=5`$ a.u. and $`V_0=0.125`$ H the potential has a single $`m=0`$ bound state while for $`V_0=0.25`$ H there are two $`m=0`$ bound states. In both cases we see that the total screening density is well-approximated by the asymptotic result in Eq. (10).

To include the effect of interactions, we identify $`\mathrm{\Delta }v_{\mathrm{eff}}`$ in Eq. (2) with $`V(r)`$ and make use of Eq. (10). In the high-density limit we then find

$$\mathrm{\Delta }v_{\mathrm{eff}}(q)=\frac{2\pi Z}{q+2},$$
(11)
which in real space gives the Thomas-Fermi potential. This potential is purely attractive and has a bound state eigenvalue of $`E_0=0.2862`$ H which is the $`r_s0`$ limit of the curves in Fig. 2. This explains why a bound state exists in the high density limit and why, in spite of this, the result is not in conflict with the applicability of linear response theory. It should be emphasized that the same argument in 3D leads to the conclusion that no bound state can exist in this case.

As a final practical application we present results of the calculation of the energy loss per unit length (or stopping power, $`S`$) for a projectile of charge $`Z`$ moving with velocity $`v`$ in the plane of the 2DEG. Within the so-called kinetic theory framework, the stopping power is given by the expression

$$S=n_0vv_F\sigma _{tr}(E_F),$$
(12)
where $`\sigma _{tr}(E_F)`$ is the momentum-transfer cross section defined in terms of the scattering phase shifts by

$$\sigma _{tr}(E_F)=\frac{4}{v_F}\underset{m=0}{\overset{\mathrm{}}{}}\mathrm{sin}^2[\eta _m(E_F)\eta _{m+1}(E_F)].$$
(13)
To leading order in the velocity, it is sufficient to determine the scattering phase shifts using the static nonlinearly screened potentials calculated in the present paper.

In Fig. 4 we show the stopping power as a function of the projectile charge $`Z`$. For small $`Z`$, $`S`$ has the expansion $`S=S_1Z^2+S_2Z^3+\mathrm{}`$ where the first two terms are the linear and quadratic response results, respectively. To emphasize the deviations from linear response, we present the results in the form $`S/(vZ^2)`$. In this representation, the stopping power including the quadratic response correction appears as a straight line with slope $`S_2/v`$. This correction was previously calculated within the quadratic random phase approximation and is shown in Fig. 4 as the straight line. We can see that corrections beyond quadratic response theory are large, especially for positive charges. Furthermore, the inclusion of xc is seen to enhance the stopping power considerably in the range $`1Z1`$, even to a greater extent than found in 3D. Finally, comparison with earlier calculations demonstrates that important differences arise when the self-consistently determined nonlinear screening potentials are used to evaluate the momentum scattering cross-sections.

In summary, we have performed self-consistent calculations of the nonlinear screening of a point charge in a 2DEG using density functional theory. We have also proved a screening theorem which clarifies the behaviour of the screening in the high density limit. These results find application in a variety of problems, such as charged impurity scattering and the stopping power of charged projectiles.

The work of I.N. has been supported by the OTKA (Grant Nos. T025019 and T029813), and that of E.Z. by a grant from the Natural Sciences and Engineering Research Council of Canada. P.M.E. thanks the University of the Basque Country, the Basque Hezkuntza, Unibertsitate eta Ikerketa Saila and the Spanish MCyT for support.


text:
# Magneto- and Electrotransport Properties of Nd0.67Sr0.33Mn1-xFexO3 (𝑥 = 0.00, 0.05)

## Abstract

Simultaneous studies on the magnetic and electrotransport properties of Nd0.67Sr0.33Mn1-xFexO3 polycrystalline bulk and epitaxial thin films ($`x`$ = 0.00 and 0.05) have been carried out. A magnetoresistance (MR) as high as $`33\%`$ is observed for Nd0.67Sr0.33MnO3 bulk at the metal-insulator transition (MIT) temperature ($`T_p`$) of 273 K in a magnetic field (H) of 10 kOe. Fe substitution at Mn sites reduces the ferromagnetic (FM) ordering temperature ($`T_c`$) and leads to an overall increase in MR. MR as high as $`90\%`$ is observed for thin film ($`x`$ = 0.05) at H = 10 kOe with $`T_p`$ = 100 K while the corresponding bulk has a MR of only $`43\%`$. Low field MR in epitaxial thin films also display a pronounced anisotropy effect compared to bulk materials.


text:
# Electronic correlation in cyclic polyenes. Behavior of approximate coupled-pair theories for large rings

## I Introduction

The problem of electronic correlation in the $`\pi `$-electron systems corresponding to small cyclic polyenes (annulenes) has been studied extensively within the Pariser–Parr-Pople (PPP) model by Paldus et al. paldus:84a ; paldus:84b ; paldus:82 ; takahashi:83 ; paldus:84c ; takahashi:85 ; piecuch:90a ; piecuch:90b ; piecuch:90c ; piecuch:91 ; paldus:92 ; piecuch:92 . They have found that the coupled-cluster doubles method (CCD) \[in annulenes equivalent to the coupled-cluster singles and doubles (CCSD) one\] breaks down when the correlation effects become sufficiently strong. They devised the approximate coupled-pair method corrected for connected quadruple excited clusters (ACPQ) paldus:84a ; paldus:84b . The ACPQ method, together with the similar ACP-D45 method (or ACP, in short) jankowski:80 , was shown to perform very well for the small annulenes, being convergent and giving correlation energies close to the full configuration interaction (FCI) results.

Recently, small annulenes have been studied with the inclusion of the doubles (D), triples (T), and quadruples (Q) in the cluster operator podeszwa:02a . It turned out that although the inclusion of quadruple excitations considerably improves the results in the weakly and moderately correlated regime (which correspond to large and moderate absolute values of the PPP resonance integral $`\beta `$, respectively), the CCDT, CCDQ, and CCDTQ methods break down in the strongly correlated regime. It was also found that in the strongly correlated regime the $`t_2`$ amplitudes of the ACP and ACPQ methods differ considerably from the FCI $`t_2`$ amplitudes.

Even though the ACP and ACPQ methods are satisfactory for small annulenes, it has been unknown how the methods perform for large annulenes, where the quasi-degeneracy problems may be too strong. The strength of the correlation (and the level of quasi-degeneracy) can be adjusted by decreasing the absolute value of $`\beta `$ for small annulenes, but it has been unknown whether this procedure is equivalent to the increasing of the size of the ring and keeping the resonance integral constant. The behavior of the methods for large systems is important if one wants to extrapolate the results and study the limit of the one-dimensional metal. That motivated us to investigate the performance of the ACP and ACPQ methods for large annulenes.

## II Theory

In this Letter, we provide only a brief description of the PPP model of annulenes; more details may be found in Ref. podeszwa:02a, .

We shall consider annulenes of formula $`\mathrm{C}_M\mathrm{H}_M`$, where $`M=4m_0+2`$, $`m_0=1,2,\mathrm{}`$, and the number of electrons $`N=M`$. The C atoms form a regular polygon, and the C–C bonds are assumed to be of the length $`R^0=1.4`$ Å. The Fock-space Hamiltonian $`\widehat{H}`$, built according to the prescriptions of the PPP model, is given in Eq. (2) of Ref. podeszwa:02a, . The following semiempirical parameters are used: the Coulomb integral $`\alpha =0`$, the resonance integral $`\beta =2.5`$ eV (the so-called spectroscopic value), and the two-center two-electron integrals $`\gamma _{mn}`$ are parametrized with the Mataga–Nishimoto formula mataga:57 , $`\gamma (R)=e^2[R+e^2(\gamma ^0)^1]^1`$, where $`e`$ is the electron charge and $`\gamma ^0=\gamma (0)=10.84`$ eV. Within the PPP computer code, the atomic units are used: $`1\text{ bohr}=0.529177\text{ Å}`$, $`1\text{ hartree}=27.2116\text{ eV}`$ (the conversion factors are the same as employed in Refs. stolarczyk:88,  and podeszwa:02a, ).

The ACP and ACPQ methods are approximations to the CCD method. In this method we use exponential Ansatz and expand the ground-state wave function $`\mathrm{\Psi }`$ by using the cluster operator $`\widehat{T}=\widehat{T}_2`$,

$$\mathrm{\Psi }=\mathrm{exp}(\widehat{T_2})\mathrm{\Phi },$$
(1)
where $`\mathrm{\Phi }`$ is a single-determinantal reference configuration, usually the RHF wave function. The $`\widehat{T_2}`$ operator depends linearly on some parameters, the so-called $`t_2`$ amplitudes, that can be obtained by solving the set of nonlinear CCD equations given by

$$\mathrm{\Phi }^{}|\mathrm{exp}(\widehat{T_2})\widehat{H}\mathrm{exp}(\widehat{T}_2)\mathrm{\Phi }=0,$$
(2)
where $`\mathrm{\Phi }^{}`$ represent all the doubly excited configurations. The correlation energy for the ground state is completely determined by the $`t_2`$ amplitudes and can be calculated as

$$E_{\mathrm{corr}}=\mathrm{\Phi }|\widehat{H}\widehat{T}_2\mathrm{\Phi }E_{\mathrm{HF}}.$$
(3)
The cyclic symmetry of our PPP model imposes some extra conditions on the CC model. The $`t_1`$ amplitudes (singles) must vanish, and therefore $`\mathrm{\Phi }`$ is also the Brueckner determinantal function paldus:73 ; stolarczyk:84 . Thus, the CCD method becomes equivalent to the CCSD one. In the presence of the cyclic symmetry, the amplitudes depend on 3 indices (instead of 4), and the number of amplitudes scales as $`M^3`$ (instead of $`M^4`$). Moreover, the computational cost for the CCD iteration scales as $`M^4`$ instead of $`M^6`$. This makes enormous savings in the CPU time for large annulenes and makes them computationally accessible. The ACP and ACPQ methods violate the alternancy symmetry koutecky:85 , and therefore we have not used this symmetry in the present calculations. The details of the implementation of the cyclic symmetry are given in Ref. podeszwa:02a, .

The explicit form of the CCD equations is most conveniently presented in the form of diagrams (see, e.g., Ref. kucharski:86, ). The left hand side of Eq. (2) is then a sum of terms represented by these diagrams. In the ACP method, some of the quadratic terms of the CCD equations are omitted (these are the diagrams that do not factorize with respect to the hole lines). It was shown paldus:84b  that these terms are approximately canceled by the terms corresponding to the connected quadruple excitations. Thus, the omission of these terms may improve the results and was shown to be very effective for the annulenes in the strongly correlated regime paldus:84b . The only quadratic terms that are present in the ACP method correspond to the diagrams in Fig. 1 (the labels below are the same as in Ref. paldus:84a, , where the complete set of the diagrams may be found).

For closed-shell calculations, the CC equations may be put in a spin-adapted form. We use the nonorthogonally spin-adapted formalism cizek:66 , which for the CCD method and its variants is equivalent to the orthogonally spin-adapted formalism used by Paldus et al. paldus:84b . With the latter formalism they devised an improvement to the ACP method called ACPQ. In this method the diagram (5) of Fig. 1 is multiplied by 9 in the triplet coupled equations of Ref. paldus:84b, . The ACPQ method better approximates the quadruples and was shown to provide the exact solution in the limit of $`\beta 0`$ paldus:84b .

The ACPQ method can also be translated to the nonorthogonally spin-adapted formalism. Let $`g_{ij}^{ab}`$ represent the nonorthogonally spin-adapted term corresponding to diagram (5) of Fig. 1 in the ACP method ($`i,j`$ and $`a,b`$ stand for the occupied and unoccupied orbitals, respectively). In the ACPQ method, the term should be modified in the following way:

$$(\text{ACP})g_{ij}^{ab}5g_{ij}^{ab}4g_{ji}^{ab}(\text{ACPQ})$$
(4)
## III Results

We performed the ACP and ACPQ calculations for the annulenes, $`\mathrm{C}_M\mathrm{H}_M`$, for $`M`$ up to $`442`$. In Fig. 2, the ACP and ACPQ results for small annulenes are compared with the results of ACPQ+T(ACPQ) paldus:92 , FCI results bendazzoli:94 , and the density matrix renormalization group (DMRG) calculations fano:98 . These are the methods that were shown to provide converged results for the small annulenes. It can be seen that the correlation energy corrected for the triples in the ACPQ+T(ACPQ) method is much better than in the ACPQ method. The method depends, however, on the $`t_2`$ amplitudes from the ACPQ method, and will work only if the latter method converges. For small annulenes the DMRG results are the most accurate and the deviation from the FCI is so small that it cannot be seen on the figure.

Since the ACP and ACPQ correlation energies do not change much with the increasing size of the ring around $`M=50`$, it may suggest that the results are close to the saturation. However, it is not true for the larger annulenes, as shown in Fig. 3. Some numerical values are also shown in Table 1. Surprisingly, the correlation energy is not monotonic with the increasing size of the ring and reaches a minimum. It means that the extrapolation of the results for smaller rings would lead to a completely wrong result.

We found that the ACP method is convergent up to $`M=442`$. For $`M>270`$, the use of the direct inversion in the iterative subspace method (DIIS) pulay:82 ; scuseria:86  was necessary to provide converged results. The DIIS method was also very useful in accelerating the convergence for the smaller annulenes.

The results of the ACPQ method are even more complicated. For $`M<174`$ there was one solution found, but between 174 and 194 we found two solutions. One set of solutions, later referred to as ACPQ1, is a continuation of the results for $`M<174`$. It was found by using the ACP result for a given $`M`$ as a starting point for the iterative process with the DIIS method. The other set of solutions, of lower energy, was found by applying the DIIS iterative scheme to the MP2 starting point. For $`M<174`$, both starting points resulted in the same solution (ACPQ1).

None of the ACPQ1 and ACPQ2 converged for $`M198`$. In Fig. 4, the results of the ACPQ method for $`M`$ from $`150`$ to $`198`$ are shown. It can be seen that the correlation energies of the two solutions approach each other and the states become quasi-degenerate for $`M198`$. It is interesting to check whether the amplitudes of the states also approach each other. To test the similarity of the $`t_2`$ amplitudes of the two solutions, we calculated the parameters $`\theta `$ and $`\eta `$ defined in Eq. (22) of Ref. podeszwa:02a, . Here $`\theta `$ measures the angle between the vectors formed of the $`t_2`$ amplitudes, and $`\eta `$—the ratio of the vector lengths. For $`M=174`$, we found $`\theta =2.69^{}`$ and $`\eta =0.9989`$, while for $`M=194`$, $`\theta =0.74^{}`$ and $`\eta =0.9999`$. This shows that the $`t_2`$ amplitudes of the ACPQ1 and the ACPQ2 solutions are quite similar for $`M=174`$, and are almost identical for $`M=194`$.

The existence of multiple solutions of the CC equations has been observed previously for the H4 model meissner:93 ; jankowski:94 , and the complete set of solutions has been obtained with the homotopy method kowalski:98a . In those studies, one of the solutions was easily identified as the ground-state solution. In contrast, for the ACPQ method of our annulene, such identification is ambiguous. The ACPQ1 solution, which exists for small annulenes, and which may be considered as the standard solution, has the energy higher than the ACPQ2 solution. Some exotic solutions that have the energy lower than the ground state have been observed for the PPP model of benzene podeszwa:02b . Our preliminary results for $`\mathrm{C}_{10}\mathrm{H}_{10}`$ also show the existence of such solutions podeszwa:02d .

## IV Conclusions

The ACP an ACPQ methods, while performing satisfactorily for small annulenes, have convergence problems for large annulenes. The ACP and ACPQ correlation energies are not monotonic with the increasing size of the ring. It shows that an extreme caution must be taken if one tries to extrapolate the results of the small annulenes to the infinite limit. In the case of the ACPQ method, two solutions are found for $`174M194`$, behaving as if they coalesce into a single solution between $`M=194`$ and $`M=198`$. However, no converged real solution of the ACPQ equations has been found beyond $`M=194`$, which may indicate that the solutions become complex for $`M198`$. The above observations indicate that the ACP and ACPQ methods cannot be applied to the one-dimensional metallic case.

Among all the methods used in the calculations of the PPP annulenes, the DMRG method seems the most promising. However, it should be noted that the method was not tested in a broader range of the resonance integral $`\beta `$, especially in the strongly correlated regime of $`\beta 0`$, where the various CC methods break down. It is also unclear whether the method can be applied to the large or infinite annulenes, since the method may not be size-extensive for the metallic case chan:02 .

The existence of a method that can describe the large annulenes and the one-dimensional metal in the PPP model is still an open problem. It is unlikely that the method that cannot describe the strongly correlated regime of a small annulene would be able to describe the infinite case. It seems, however, that the proper behavior of a method in the strongly correlated regime of small annulenes still does not guarantee the correct description of the larger systems.

I would like to express my gratitude to L. Z. Stolarczyk for encouragement and helpful discussions. I would also like to thank B. Jeziorski and L. Z. Stolarczyk for the critical reading of the manuscript. The work was supported by the Committee for Scientific Research (KBN) through Grant No. 7 T09A 019 20.


text:
# Hot White Dwarf Donors in Ultracompact X-Ray Binaries

## 1. INTRODUCTION

The recent discovery by the Rossi X-Ray Timing Explorer of two accreting millisecond pulsar transients XTE J1751-305 ($`\nu _s435\mathrm{Hz}`$; Markwardt et al. 2002) and XTE J0929-314 ($`\nu _s185\mathrm{Hz}`$; Galloway et al. 2002) has allowed us to learn about the donors in these ultracompact binaries. At $`P_{\mathrm{orb}}43`$ minute orbital periods, H-rich donors are ruled out (Nelson, Rappaport & Joss 1986), and He-rich stars (see Podsiadlowski, Rappaport & Pfahl 2002 for an updated discussion of the range of H/He ratios) or WDs are filling the Roche lobe (RL). If a degenerate WD, then the orbital period increases as mass is transferred at the rate set by angular momentum losses from gravity waves, $`\dot{J}_{\mathrm{GR}}`$ (see Verbunt 1993).

The measured pulsar orbital parameters yield the RL filling companion’s mass, $`M_\mathrm{c}`$, and the minimum mass transfer rate, $`\dot{M}_{\mathrm{GR}}=3M_\mathrm{c}\dot{J}_{\mathrm{GR}}/2J`$, which are shown in Figure 1 for neutron stars of $`M_x=1.42.0M_{}`$. The low $`\dot{M}`$’s are the likely cause for the transient behavior, as the steady-state outer disk temperature is below these element’s ionization temperature (Tsugawa & Osaki 1997; Menou et al. 2002), even if X-ray heating is included at the rate inferred in other X-ray binaries (Dubus, Hameury & Lasota 2001). If identical WDs are the donors in both of these binaries, then their identical orbital periods would require the same RL filling solutions. This would constrain the inclination for XTE J0929-314 to less than 37 degrees and $`M_c>0.013M_{}`$ for both XTE J0929-314 and XTE J1751-305.

The dotted lines in Figure 2 show the mass-radius relation for RL filling donors of XTE J1751-305 (lower line) and XTE J0929-314 (upper line). The solid and dashed lines (which exhibit a maximum radius due to the onset of Coulomb physics) are the cold ($`T=0`$) WDs of Zapolsky & Salpeter (1969, hereafter ZS) for pure He and C, respectively. A cold He WD will fill the RL for XTE J0929-314 (Galloway et al. 2002), whereas there are no cold WD solutions that fill the RL for XTE J1751-305 (Markwardt et al. 2002).

This evidence for finite $`T_c`$ WDs motivated my calculations shown in Figure 2 by the solid and dashed lines that diverge at low $`M_c`$. These models (see §2) are for He at $`T_c=10^5\mathrm{K}`$ and $`10^6\mathrm{K}`$ and C at $`T_c=10^6\mathrm{K}`$ and $`T_c=3\times 10^6\mathrm{K}`$. These provide RL filling solutions, but don’t a priori differentiate between He or C or any other element that might dominate the donor star. I consider C WDs throughout this paper, as Schulz et al. (2001) have measured a high Ne to O ratio in the matter transferred onto NSs in ultracompact X-ray binaries (confirmed by Juett, Psaltis & Chakrabarty 2002; Homer et al. 2002; Juett & Chakrabarty 2002). These measurements led Schulz et al. (2001) and Juett et al. (2001) to suggest that the donors in these binaries are the cores of previously crystallized C/O WDs.

In §2 I justify (on evolutionary grounds) that hot WDs are expected, describe the construction of their $`R_c(M_c)`$ relations, and point out the existence of a minimum $`M_c`$ solution for a given $`T_c`$. I put the models in the context of mass donors in §3 and calculate

$$n_{\mathrm{Ad}}=\frac{d\mathrm{ln}R_c}{d\mathrm{ln}M_c},$$
(1)
the adiabatic exponent needed to evaluate the stability of mass transfer. Though uncertainties remain regarding the ability of the accretion disk to store angular momentum on long timescales (Verbunt & Rappaport 1988), my initial work finds that the expansion of a He WD donor under mass transfer can exceed that of the RL, possibly allowing for the Ruderman & Shaham (1983) instability.

## 2. FINITE ENTROPY WHITE DWARFS AND EVAPORATION

The material that fills the RL in these systems today was once the deep interior of an initially more massive WD, so I begin by discussing the extent to which prior evolution (before RL filling) can affect the temperature of the current low-mass WDs. There are a few evolutionary scenarios that place a He or C/O WD in a tight enough orbit (following a common-envelope induced spiral-in) about a NS that gravity wave emission will place it in contact within 5-10 Gyr (Rasio, Pfahl & Rappaport 2000; Yungelson, Nelemans & van den Heuvel 2002). In this case the WD has had some time to cool before initiating RL overflow. Nelemans et al. (2001) discuss a similar scenario for the origin of AM CVn binaries.

The rate at which they cool differentiates He and C/O WDs. The larger specific heat of a He WD slows its cooling (Althaus & Benvenuto 1997; Hansen & Phinney 1998). For example, a $`0.20M_{}`$ He WD would have $`T_c=10(3.3)\times 10^6\mathrm{K}`$ at $`1.0(4.0)`$ Gyr (Althaus & Benvenuto 1997) and be liquid when RL filling occurs. The smaller specific heat of a $`0.6M_{}`$ C/O WD allows it to cool to $`2.5\times 10^6\mathrm{K}`$ in 4 Gyr and begin crystallization (Salaris et al. 2000). Hence, most C/O WDs would be solid prior to RL filling (this is needed for the enhancement of $`{}_{}{}^{22}\mathrm{Ne}`$ by fractionation; Schulz et al. 2001) and might only melt if heated during the GW inspiral (Iben, Tutukov & Fedorova 1998). The other differentiation between such WDs comes from their mass, as Yungelson et al. (2002) argue that the mass transfer at the onset of RL filling will be unstable if $`M_c>0.44M_{}`$, excluding traditional C/O WDs from ever reaching short orbital periods. However, Yungelson et al. (2002) note that “hybrid” WDs with C/O cores can have $`M_c<0.44M_{}`$ and thus evolve to lower mass systems filling the RL at 20-40 minutes.

The inability for the donor WDs to cool on the $``$ Gyr timescale of the mass transfer phase (Rappaport et al. 1987) makes it clear that the initial entropy is the minimum value attainable. All of the entropy is in the liquid ions, which have an adiabatic scaling, $`T_c\rho _c^{0.50.6}`$ (Hernanz et al. 1988), less steep than for an ideal gas ($`T\rho ^{2/3}`$). Since $`\rho _cM_c^2`$, a $`0.013M_{}`$ He WD made by adiabatically expanding the deep interior of a $`0.20M_{}`$ WD would be a factor of $`15`$ cooler, or $`T_c6.6(2.2)\times 10^5\mathrm{K}`$ for the initial range of $`T`$’s discussed earlier. Of course, any tidal heating will increase $`T_c`$ (see §3).

So I have constructed $`M_c<0.03M_{}`$ WDs of finite $`T_c`$ that contains the Coulomb physics (from Chabrier & Potekhin 1988) that causes the turnover to constant density “rocks” in the mass-radius relation at $`10^3M_{}`$ seen in the ZS models in Figure 2. I have imposed an arbitrary (but convectively stable) temperature profile of $`T=T_c(P/P_c)^{1/5}`$ while integrating hydrostatic balance and mass conservation. I used the electron equation of state of Paczynski (1983) and halt integrations at the point where the pressure has fallen to $`10^5`$ of the central value, which avoids the need for an envelope model. The results are shown in Figure 2 and make clear that the temperatures expected from adiabatic evolution of a He WD are adequate to provide the slight radius expansion needed to fill the RL for these accreting millisecond pulsars. This is easier for He WDs because of the extra number of ions per electron compared to C/O WDs. White dwarfs made of C/O will need some tidal heating to fill these RLs.

I also found that there is a minimum mass solution for a WD of fixed $`T_c`$. This is related to the electron Fermi energy of the nearly constant density solutions at low masses and is reflected in the radius divergence of the models in Figure 2 at low $`M_c`$. My semi-analytic modeling finds that these finite $`T_c`$ solutions eventually “turnover” to follow a track of $`R_cM_c`$, as expected for an ideal gas polytrope (where $`T_cM_c/R_c`$). Since there is no way to reach these solutions of higher $`M_c`$ while mass is being lost, these “evaporative” endpoints where $`d\mathrm{ln}R_c/d\mathrm{ln}M_c`$ diverges can only be reached if tidal heating is adequate to keep these models hot while the central density is dropping. More work remains to actually show that this will eliminate the companion and lead to an isolated millisecond radio pulsar.

## 3. MASS TRANSFER INSTABILITIES AND QUIESCENT EMISSION

Barring tidal heating adequate to reach the “evaporative” endpoint alluded to in §2, I will now discuss a possible mass transfer instability. Matter leaving the WD takes angular momentum with it as it settles into the accretion disk. This angular momentum is returned to the WD via tidal torques once there has been enough time for viscosity to move the material outwards from the splash point. The timescale for this angular momentum loop to be closed has been considered to be long enough (or uncertain enough) that Ruderman & Shaham (1983) raised the possibility that a mass transfer instability could occur when the expansion of a low-mass donor ($`M_cM_x`$) due to mass loss (measured by $`n_{\mathrm{Ad}}`$ of eq. ) exceeds that of the RL, measured by

$$n_\mathrm{R}=\frac{d\mathrm{ln}R_R}{d\mathrm{ln}M_c}=\frac{d\mathrm{ln}a}{d\mathrm{ln}M_c}+\frac{1}{3},$$
(2)
where $`a`$ is the orbital separation. The value of $`n_\mathrm{R}`$ presuming the angular momentum “sink” of the disk (using the fitting formula of Verbunt & Rappaport 1988) is shown in Figure 3 by the dotted lines for (from left to right) $`M_x=0.6,1.0,1.4,`$ & $`1.8M_{}`$.

Ruderman & Shaham (1983) presumed $`n_{\mathrm{Ad}}=1/3`$, which was shown to be an overestimate by Hut & Paczynski (1984) and Bonsema & van den Heuvel (1985), who found $`n`$ from the ZS relation (shown in Figure 3 by the solid\[dashed\] lines for He\[C\]). These authors noted that a ZS model moved the crossing ($`n_{\mathrm{Ad}}<n_\mathrm{R}`$) to such low $`M_c`$ that the binary might not reach it in a Hubble time. However, for finite $`T`$ He WDs, I show that this conclusion is altered. The solid circles are my evaluation of $`n_{\mathrm{AD}}`$ for the $`T_c=10^5\mathrm{K}`$ (set closest to the solid line) and $`T_c=10^6\mathrm{K}`$ (bottom-most set) He WDs of Figure 2 and show that the finite entropy model sits between the $`T=0`$ value and the “naive” $`n=1/3`$ guess for a perfect gas. For a $`M_x=1.4M_{}`$ accretor, the instability question is first raised (since $`n_{\mathrm{Ad}}<n_\mathrm{R}`$) when a $`10^6\mathrm{K}`$ He WD has $`M_c0.01M_{}`$. This will occur about a Gyr after mass transfer has started presuming just gravity wave emission. My evaluations of $`n_{\mathrm{Ad}}`$ for C WDs at the $`T_c`$’s of Figure 2 never found such a crossing.

The instability question is thus raised for He WDs and requires some extra entropy from tidal heating. The amount of heating is small, as for liquid He to reach $`T_c10^6\mathrm{K}`$ from about 1/2 that value (which it had initially) requires $`3k_B(0.5T_c)/4m_p3\times 10^{13}\mathrm{erg}\mathrm{g}^1`$ (the specific heat in the liquid state is $`3k_BT`$ per ion). This energy is comparable to the current rotational energy per gram for a tidally locked WD at the 40 minute orbital period. Hence, only a fraction of the higher rotational energy per gram from the tidally locked WD at a shorter orbital period needs to be deposited to yield $`T_c10^6\mathrm{K}`$ today.

The remaining criticism against the mass transfer instability is one of relative timescales (Verbunt & Rappaport 1988; Priedhorsky & Verbunt 1988). The mass-transfer instability most likely grows on a timescale $`\tau _gM_cH/\dot{M}_cR_c`$ (Verbunt & Rappaport 1988), where $`H`$ is the scale height in the WD atmosphere, fixed by the X-ray heating in quiescence from the hot NS (Bildsten & Chakrabarty 2001). The binary spends most of the time undergoing mass transfer to the outer disk, but little to no accretion onto the NS. Bildsten & Chakrabarty (2001) showed that the X-ray emission always detected from NS binaries in quiescence (see Bildsten & Rutledge 2001 for an overview) heat the companion on the side facing the NS. Presuming the NS thermal emission is at the level predicted by Brown, Bildsten & Rutledge (1998) (giving $`L_{X,q}5\times 10^{32}\mathrm{erg}\mathrm{s}^1`$) then the WDs in XTE J1751-305 and XTE J0929-314 will have $`T_{\mathrm{Eff}}50006500\mathrm{K}`$.111These high $`T_{\mathrm{Eff}}`$’s should aid detection of the quiescent counterpart, as the luminosity (presuming $`R_c=0.045R_{}`$) will be $`L_c(0.51.0)\times 10^3L_{}`$ just from reprocessing the NS’s thermal emission. Presuming $`T_{\mathrm{Eff}}6000\mathrm{K}`$ and using the WD radius implied by RL filling, I get $`M_V12`$ for the heated WD face. For the high latitude source XTE J0929-314 at $`d57`$ kpc (Galloway et al. 2002), this gives $`V26`$.  The scale height is $`H/R_c3\times 10^4`$, yielding $`\tau _g10^6\mathrm{yr}`$. This is still much longer than the time for viscosity to finally play a role in moving material outwards in the outer accretion disk, most likely allowing for the angular momentum to get back to the donor and move it out, so that $`n_\mathrm{R}=5/3`$. The adiabatic index is never less than $`5/3`$. However, the isothermal WD response is $`<5/3`$ and mass transfer would be unstable for $`M_c<4(9.5)\times 10^3M_{}`$ in a He WD of $`T_c=4(10)\times 10^5\mathrm{K}`$ as long as internal heating can keep the WD isothermal as mass is lost. These masses are near the values I noted earlier for the “evaporation”, namely where the mass-radius relation turns over.

## 4. CONCLUSIONS

The timing of two accreting millisecond pulsars (Markwardt et al. 2002; Galloway et al. 2002) in ultracompact binaries has probed the WD donor properties to new levels and shown that they are of finite entropy. This motivated my calculations of low-mass WDs of finite $`T_c`$ that allow for $`T_c`$ to be constrained. For He WDs, the implied $`T_c`$ are nearly that expected just from adiabatic expansion of the initially hot WD that filled the RL. Only a small amount of tidal heating is needed. More tidal heating is needed to make a C/O WD fill the RL.

These finite $`T`$ solutions allowed for a re-evaluation of Ruderman & Shaham’s (1983) scenario for making isolated millisecond radio pulsars via a mass transfer instability. I find that the adiabatic mass transfer instability can occur for a hot He WD as long as the angular momentum leaving the RL filling star is not returned to the orbit. I have thus eliminated one criticism of their model, though the question of angular momentum elimination remains a serious one. The physics of hot, low-mass WDs also yields a minimum mass WD solution for a fixed $`T_c`$ so that a mass transfer instability can occur if the donor remains isothermal under mass loss. I thus speculate that an evaporative or mass transfer instability endpoint might occur as long as tidal (or other) heating persists at 40-80 minute orbital periods.

This work also impacts AM CVn binaries, where a low mass He star donates material to a more massive WD (see Solheim 1995 for a review). The larger WD radii lead to more GW emission and a higher $`\dot{M}_c`$ than expected for a given $`P_{\mathrm{orb}}`$. Hence, models which track the WD entropy will fall between the degenerate and non-degenerate models in Nelemans et al. (2001) and depend on both the age of the system when RL filling occurs (as this fixes the initial WD entropy) and any tidal heating that occurs during the mass transfer. If either of the instabilities discussed above occur, then the endpoint of AM CVn’s could well be a DB WD (e.g. Tutukov & Yungelson 1996).

I thank Deepto Chakrabarty for alerting me to the discovery of these transients and for many conversations. Ira Wasserman provided great physics insights during the progress of this work, which was supported by NASA via grant NAG 5-8658 and by the NSF under grants PHY99-07949 and AST01-9642. L. B. is a Cottrell Scholar of the Research Corporation.