text:
# Discovery and quantitative spectral analysis of an Ofpe/WN9 (WN11) star in the Sculptor spiral galaxy NGC 3001footnote 11footnote 1Based on observations obtained at the ESO Very Large Telescope [65.N-0389(B) and 67.D-0006(A)].

## 1 Introduction

With the new telescopes of the 8-10 meter class stellar astronomy is branching out beyond the Local Group. Ideal targets for our understanding of young stellar populations in distant galaxies are hot massive stars. These objects have strong stellar winds producing broad and easily detectable spectral features distributed over the whole wavelength range from the UV to the IR, and providing unique information on chemical composition, galactic evolution and extragalactic distances. With these perspectives in mind we have recently begun a systematic spectroscopic study of luminous blue stars in galaxies beyond the Local Group, and presented spectral classification and first quantitative results for A supergiants in NGC 3621 (6.7 Mpc, Bresolin et al. 2001) and NGC 300 (2.0 Mpc, Bresolin et al. 2002). Here we report on the discovery and detailed quantitative analysis of the first example of an Ofpe/WN9 star outside of the Local Group. We will present a detailed chemical abundance pattern – the first in a galaxy beyond the Local Group – together with stellar parameters and a determination of the stellar wind properties.

The Ofpe/WN9 class was introduced to include objects which show in their spectra high excitation emission lines from He ii and N iii, typical of Of stars, together with low excitation lines from He i and N ii, seen in late WN stars (Walborn 1982, Bohannan & Walborn 1989). Objects of this class have so far been identified in the Galaxy (possibly several stars in the Galactic center: Allen et al. 1990, Najarro et al. 1997, Figer et al. 1999), the LMC (ten stars: Bohannan & Walborn 1989), M33 (seven stars: Massey et al. 1996, Crowther et al. 1997) and M31 (one star: Massey & Johnson 1998). The importance of Ofpe/WN9 stars as objects in a transitional stage of evolution between O and W-R stars has been recognized in the last decade, and a connection to the LBV class has been suggested, Ofpe/WN9 stars being observed during a quiescent or post-LBV phase (Crowther & Smith 1997, Pasquali et al. 1997). Indeed, in at least a couple of instances Ofpe/WN9 stars have been observed to turn into LBVs (R127: Stahl et al. 1983; HDE 269582: Bohannan 1989). The LBV AG Car is also known to show an Ofpe/WN9-like spectrum during its hot phase at visual minimum (Stahl 1986). The discovery of ejected circumstellar nebulae, in some cases measured in a state of expansion, associated with some of the LMC Ofpe/WN9 stars by Nota et al. (1996) and Pasquali et al. (1999) brings forward strong evidence for the occurence of violent episodes of mass loss in these stars, similar to the shell-producing, eruptive outbursts of LBVs.

Smith et al. (1994) revised and extended the classification of late WN stars to include lower excitation objects. In their scheme, stars showing spectra like AG Car at minimum, where no N iii is detected, are reclassified as WN11. Given the spectral resemblance of the NGC 300 star we analyze here to AG Car at visual minimum (see § 4) we will adopt the WN11 classification for the remainder of this Letter. We describe the observational data in § 2, and the photometry in § 3. In § 4 we present our VLT spectra, together with the stellar and wind properties derived from a quantitative spectral analysis.

## 2 Observations

The data presented here are part of a spectroscopic survey of photometrically selected blue supergiants in the Sculptor spiral galaxy NGC 300 obtained at the VLT with the FORS multiobject spectrograph, and described in detail by Bresolin et al. (2002). This latter paper presents spectra obtained in September 2000 in the blue spectral region ($``$ 4000-5000 Å) and used for the spectral classification of about 70 supergiants. The emission-line star analyzed here corresponds to star B-16 of the spectral catalog (see Table 2 and finding charts in the aforementioned paper), a rather isolated bright star apparently not associated with any prominent OB association, star cluster or nebulosity. The nearest OB association (from the catalog of Pietrzyński et al. 2001) lies at a deprojected distance of $``$450 pc. Given the estimated age of 3.8 Myr (§4.1) a speed on the order of 100 km/s would be required for the star to reach its current position from this association. This is much larger than the typical stellar random velocities, and comparable to that of Galactic runaway stars. To the best of our knowledge the star is not included in published listings of emission-line or W-R stars in NGC 300 (Schild & Testor 1992, Breysacher et al. 1997, and references therein), and we will refer to it as star B-16.

In a more recent observing run (September 2001), in order to measure the mass-loss rates of blue supergiants in NGC 300 from the H$`\alpha `$ line profile, we secured MOS spectra in the red with a 600 lines mm-1 grating and 1″ slitlets, which, together with the older spectra, provides us with complete coverage at $`R1,000`$ resolution of the 3800-7200 Å wavelength range, except for a 50 Å-wide gap centered at 5025 Å. The total exposure time was 13,500 s for both the blue and the red spectra, while better seeing conditions favored the blue observations. The mass-loss rates we derive for the whole blue supergiant sample will be discussed in a forthcoming paper, and we concentrate here on the single B-16 star.

## 3 Photometry

As part of a multi-epoch photometric study of NGC 300 carried out for the discovery and monitoring of Cepheid variables (Pietrzyński et al. 2001, see also Pietrzyński et al. 2002), $`BV`$ photometry from ESO/MPI 2.2m telescope WFI images is available for B-16 covering nearly a one-half year period. The resulting $`B`$ and $`V`$ light curves are shown in Fig. 1. An apparently irregular variability, on the order of 0.1 mag (peak-to-peak), is detected in both bands, at a significant level above the typical photometric uncertainty of 0.015 mag. The color index remains approximately constant around $`BV=0.07`$. While bona fide LBV variables show large amplitude variability on timescales of years (Stahl et al. 2001), partly characterizing the outburst phenomena in these objects, smaller-amplitude variability is observed in Ofpe/WN9 stars or LBVs during quiescence (e.g., Stahl et al. 1984, Sterken et al. 1997). Assuming the intrinsic color index from the stellar atmosphere models ($`E_{BV}=0.06`$), $`R=A_V/E_{BV}=3.1`$ and an average magnitude $`V=19.00`$, we obtain from the adopted distance modulus for NGC 300, $`mM=26.53`$ (Freedman et al. 2001), an absolute magnitude $`M_V7.72\pm 0.13`$.

## 4 Spectral analysis

The blue and red portions of the rectified VLT optical spectrum are shown in Fig. 2. In the top panel we have also superimposed the spectrum of AG Car during minimum, taken from the Walborn & Fitzpatrick (2000) atlas, in order to illustrate the remarkable similarity between the two spectra. Line identification is provided for most of the recognizable features, which characterize B-16 as a WN11 star. These include the hydrogen Balmer series lines (pure emission), He i lines (together with He ii$`\lambda `$ 4686) and mostly low-excitation metal lines such as N ii, Fe iii and Si ii. P-Cygni profiles are seen in most of the He i lines, as well as in Si iii$`\lambda \lambda `$ 4552-4575, N ii$`\lambda \lambda `$ 4601-4643 and Fe iii$`\lambda \lambda `$ 5127,5156. The absence of N iii is a discriminant factor against an earlier (WN9-10) classification. The WN11 classification is further confirmed by the observed equivalent widths of He i$`\lambda `$ 5876 (32 Å) and He ii$`\lambda `$ 4686 (1 Å), compared to those given by Crowther & Smith (1997).

### 4.1 Atmosphere models

For the quantitative analysis of the stellar spectrum we have used the iterative, non-LTE line blanketing method presented by Hillier & Miller (1998). The code solves the radiative transfer equation in the co-moving frame for the expanding atmospheres of early-type stars in spherical geometry, subject to the constraints of statistical and radiative equilibrium. Steady state is assumed, and the density structure is set by the mass-loss rate and the velocity field via the equation of continuity. The velocity field (Hillier 1989) is characterized by an isothermal effective scale height in the inner atmosphere, and becomes a $`\beta `$ law in the wind (e.g., Lamers et al. 1996). Better fits are obtained if instead of the standard $`\beta `$ law a combined 2-$`\beta `$ law is used, where each value of $`\beta `$ (given in Table 1) characterizes the shape of the velocity field in the inner and outer parts of the wind, respectively. We allow for the presence of clumping via a clumping law characterized by a volume filling factor $`f(r)`$, so that the ‘smooth’ mass-loss rate, $`\dot{M}_S`$, is related to the ‘clumped’ mass-loss rate, $`\dot{M}_C`$, through $`\dot{M}_S=\dot{M}_C/f^{1/2}`$. The model is then prescribed by the stellar radius $`R_{}`$, the stellar luminosity $`L_{}`$, the mass-loss rate $`\dot{M}`$, the velocity field $`v(r)`$, the volume filling factor $`f`$, and the abundances of the elements considered. The reader is referred to Hillier & Miller (1998, 1999) for a detailed discussion of the code.

The main stellar parameters of our best-fitting model are summarized in Table 1, whereas fits to some of the most important line diagnostics are illustrated in Fig. 3. As can be seen, the majority of the spectral features are well reproduced by the adopted model. The Fe iii$`\lambda \lambda `$ 5127,5156 lines, together with He ii$`\lambda `$ 4686, provide excellent constraints for the effective temperature. Using the strength of the electron scattering wings of the Balmer lines we derive a clumping factor $`f=0.25`$. Thus the ”smooth” mass loss rate would be twice the value displayed in Table 1. The parameters describing the stellar wind, i.e., the terminal velocity $`v_{\mathrm{}}=325`$ km s-1, the mass-loss rate $`\dot{M_C}=4.6\times 10^5`$ $`M_{}`$ yr-1, and the wind performance number $`\eta =\dot{M}v_{\mathrm{}}/(L_{}/c)=0.84`$, lie in the range found for Local Group WN11 stars (Crowther et al. 1997). We estimate the uncertainty in $`L`$ and $`\dot{M_C}/f^{1/2}`$ to be approximately 10% and 30%, respectively. As a further comparison, Fig. 4 shows the location of star B-16 in the H-R diagram, together with additional WN9-11 stars in the Local Group (Crowther, Hillier & Smith 1995; Smith, Crowther & Willis 1995; Crowther et al. 1997), as well as LBVs at minimum in the Galaxy (P Cyg: Pauldrach & Puls 1990; AG Car: Smith et al. 1994) and in the LMC (R71: Lennon et al. 1994).

Table 2 summarizes the model fractional mass abundances. Uncertainties in the latter range between 0.2 dex (He, N, Si and Fe) and 0.3 dex (C, O and Al). The chemistry of B-16 resembles that of other known Local Group WN9-11 stars. The severe depletion of hydrogen inferred from our analysis (mass fraction $`X=0.27`$), and the correspondingly large helium surface abundance (H/He = 1.5 by number), are fairly typical for late WN stars (H/He = 0.8-3 by number), and are indicative of heavy mass-loss stripping of the stellar outer layers during the post-main sequence evolution of massive stars. We find a total CNO abundance close to that of Galactic main sequence B stars (Gies & Lambert 1992, Kilian 1992), consistent with a half solar chemical composition in these elements. This agrees with the empirical O abundance derived for an H ii region located 30″ away (region 5 observed by Pagel et al. 1979). On the other hand, Fe has a roughly solar abundance, suggesting an $`\alpha `$/Fe ratio different than in the Galaxy. In comparison with the Galactic B star abundance pattern, N is overabundant (by mass) by a factor of 10, while C and O are reduced by factors of $``$2 and 6, respectively, although the latter abundances are highly uncertain.

The current evolutionary status of B-16 can be estimated by means of the recent Geneva stellar tracks including rotation (Meynet & Maeder 2000) with an initial rotational velocity $`v_{ini}=300`$ km s-1 and Z = 0.02 (Fig. 4). The star’s position in the H-R diagram approaches the 60 $`M_{}`$ track, and an initial mass of $``$55 $`M_{}`$ can also be derived from the relation between stellar luminosity and initial mass for WNL stars given by Schaerer & Maeder (1992). From the 60 $`M_{}`$ stellar model more closely matching the observed H and He surface abundances we infer an age of 3.8 Myr and a present mass of $``$36 $`M_{}`$. At this stage, the predicted N overabundance with respect to the ZAMS abundance is 12 (mostly occuring during the MS phase, as a consequence of rotational mixing), in good agreement with our finding.

According to its He- and N-enriched surface chemistry, B-16 might be in a dormant or post-LBV phase of evolution. This is also indicated by the position in the H-R diagram, well above the Humphreys-Davidson limit (Humphreys & Davidson 1979), and intermediate between the LBVs AG Car and P Cygni, which have comparable mass-loss rates to B-16, as well as similar wind performance numbers. The NGC 300 star shows higher He enrichment at the surface (nHe/n$`{}_{H}{}^{}0.7`$, versus 0.4-0.5 for the LBVs), and a somewhat faster wind (325 versus 185-250 km s-1; Najarro, Hillier & Stahl 1997, Langer et al. 1994), so that a more advanced (post-LBV) state seems more likely. This picture is in agreement with the conclusions reached by Crowther, Hillier & Smith (1995) and Crowther et al. (1997) regarding the close connection between WN9-11 stars and LBVs.

An unexpected reward out of our spectroscopic survey of blue supergiants in NGC 300, the WN11 star B-16 studied here is the first star in this galaxy (and beyond the Local Group) for which we have attempted a detailed quantitative analysis. In the near future the systematic study of the massive stellar population in NGC 300 and similar galaxies within a few Mpc, well within current observational capabilities, will certainly provide new insights into massive stellar evolution and stellar abundances.

We thank J. Hillier for providing his code, and G. Meynet for the Geneva stellar evolutionary tracks including rotation. WG gratefully acknowledges support for this research from the Chilean Center for Astrophysics FONDAP No. 15010003. FN acknowlegdes Spanish MCYT PANAYA2000-1784 and Ramon y Cajal grants.


format: text

text:
# Spectral densities of Kondo impurities in nanoscopic systems

## I Introduction

The physics at the nanoscale has emerged as one of the most active areas in condensed matter physics. This new field includes the study of small metallic and superconducting islands, quantum dots and nanostructurated semiconductors, nanotubes and quantum corrals, nanoelectronics and nano-characterization among other things.nanoscience

The advances in nanotechnologies revived the interest in the Kondo effect,Hewson1993  one of the paradigms of strongly correlated systems. On one hand, Scanning Tunneling Microscopy (STM) allowed the direct measurement of local spectroscopic properties of Kondo impurities on noble metal surfacesmadhavan98 ; nagaoka2002  and in nanoscopic systems.Manoharan2000  On the other hand, it has been shown that single electron transistors and single walled carbon nanotubes weakly coupled to contacts may behave as Kondo impurities generating new alternatives to study the phenomena.Meir1991

In its simplest form, the Kondo problem is a single quantum spin interacting with an ideal electron gas. An antiferromagnetic coupling between the impurity and the free electron spins gives rise to an anomalous scattering at the Fermi energy leading to a large impurity contribution to the resistivity. Simultaneously the impurity spin is screened by the conduction electrons and the magnetic susceptibility saturates at low temperature. There is a characteristic temperature $`T_K`$ that separates the low temperature from the high temperature regimes. At $`TT_K`$, the impurity spin is essentially free and the problem can be treated by perturbations in a dimensionless coupling constant $`\lambda `$. At $`TT_K`$ the impurity spin is screened forming a singlet complex with the conduction electrons and the system is described by an infinite effective coupling. The crossover regime with $`TT_K`$ is more difficult to describe and the best treatment corresponds to the numerical renormalization group (NRG) as done by Wilson.Wilson1975

The characteristic Kondo temperature is given by $`T_K`$ $`D\sqrt{\lambda }e^{1/\lambda }`$ where $`2D`$ is the free electron bandwidth. Associated to this energy scale there is a characteristic length scale known as the Kondo screening length $`\xi _K=\mathrm{}v_F/`$ $`k_BT_K`$ where $`v_F`$ is the Fermi velocity. The physical meaning of the screening length is that, in the low temperature regime where the impurity spin is screened, the antiferromagnetic correlations between the impurity and conduction electron spins extend up to a distance of the order of $`\xi _K`$.

The problem of Kondo impurities in nanoscopic systems for different experimental realizations has been the subject of many recent theoretical and experimental works.Thimm1999  When a Kondo system, either an atomic impurity or an artificial atom like a quantum dot (QD), is embedded in a small system of volume $`L^d`$ where $`d`$ is the spacial dimension, the length-scale $`L`$ should be compared with the characteristic Kondo length $`\xi _K`$. For $`L<\xi _K`$ finite-size effects are expected to be important. The condition $`L\xi _K`$ is equivalent to $`k_BT_K\mathrm{\Delta }`$ where the energy $`\mathrm{\Delta }`$ gives the average level spacing of the finite system. For a finite system the characteristic energy $`\mathrm{\Delta }`$ acts as a low energy cutoff for the charge and spin excitations and consequently it modifies the low temperature behavior of the system.

The ground state properties of a Kondo impurity in a small system have been addressed by a number of authors using different approximations.Thimm1999 ; Cornaglia2002 ; Affleck2001  The thermodynamic properties and the effect of coupling the system to a macroscopic reservoir have recently been studied using the Wilson’s renormalization group.Cornaglia2002  In this work we extend renormalization group calculations to evaluate the impurity spectral density and analyze how the new energy or length scales introduced by the finite size affects the Kondo resonance at the Fermi energy. We also study the temperature dependence of the low energy spectral density.

The rest of the paper is organized as follows: in section II we present the model and describe how the numerical renormalization group is adapted to our case. We then recap the most relevant thermodynamic properties. Section III contains the spectral properties of the impurity for frequencies close to the Fermi energy. After presenting the general formulation we show results for finite systems, systems in contact with a reservoir and the temperature dependence of the Kondo resonances. Finally section IV includes a summary and discussion.

## II Model and Thermodynamic properties

In this section we present the model for an impurity in a nanoscopic system coupled with a macroscopic reservoir and briefly discuss how Wilson’s renormalization group is adapted to this problem. Then we review the thermodynamic properties of the system.

### II.1 Model Hamiltonian and Wilson’s Renormalization Group

Our starting point is the Anderson model for magnetic impurities with a Hamiltonian which in the usual notation reads:Anderson1961

$`H_{AM}`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}\epsilon _dd_\sigma ^{}d_\sigma +Ud_{}^{}d_{}d_{}^{}d_{}+{\displaystyle \underset{\nu ,\sigma }{}}\epsilon _\nu c_{\nu \sigma }^{}c_{\nu \sigma }`$ (1)
$`+{\displaystyle \underset{\nu ,\sigma }{}}(V_\nu ^{}c_{\nu \sigma }^{}d_\sigma +V_\nu d_\sigma ^{}c_{\nu \sigma })\mu _iBS_{iz},`$
where the operator $`d_\sigma ^{}`$ creates an electron with spin $`\sigma `$ at the impurity orbital with energy $`\epsilon _d`$ and Coulomb repulsion $`U`$, and $`c_{\nu \sigma }^{}`$ creates an electron in an extended state with quantum numbers $`\nu `$ and $`\sigma `$ and energy $`\epsilon _\nu `$. The last term represents the effect of an external magnetic field along the $`z`$-direction coupled to the impurity spin $`𝐒_i`$. Hereafter we will use $`D=1`$ as our unit of energy.

In this notation, the nanostructure of the system is hidden in the structure of the one-electron extended states with energies $`\epsilon _\nu `$ and wavefunctions $`\psi _\nu (𝐫)`$. In equation (1) the hybridization matrix elements are taken proportional to the extended state wavefunctions at the impurity position, i.e. $`V_\nu =V_0\psi _\nu (0)`$ where the impurity position is defined as the origin of coordinates. We consider a simple structure consisting of a spherical metallic cluster of radius $`R_c`$ with the impurity at the center. The cluster is embedded in a bulk material with which it is weakly coupled through a large surface barrier. The Hamiltonian can be put in the form $`H_{AM}=H_{AM}^0+W(r,R_c)`$ where the first term is the Anderson model Hamiltonian for an impurity in an infinite homogeneous host and the last term is a spherically symmetric potential barrier placed at a distance $`R_c`$ from the impurity. Figure 1(a) illustrates the configuration described by the model. For an infinite impenetrable barrier, the central cluster is decoupled from the macroscopic reservoir and the model describes an impurity in a small sample. In this situation, the extended states that are coupled to the impurity, are confined in the central cluster and their energy spectrum is a discrete spectrum with a mean energy level separation given by the characteristic energy $`\mathrm{\Delta }.`$ For a finite barrier these states are hybridized with the continuum acquiring a finite lifetime, the local density of states inside the central cluster then presents resonances separated by the energy $`\mathrm{\Delta }`$ and with widths $`\gamma `$ that are determined by the barrier. The model then incorporates the new energy scales $`\mathrm{\Delta }`$ and $`\gamma `$, that may drastically change the impurity behavior.

This geometrical structure is particularly appropriate to use the numerical renormalization group approach developed by Wilson. Wilson introduced a logarithmic discretization of the energy of the conduction electrons, dividing the band in a series of energy intervals with exponentially decreasing width. By means of a canonical transformation the resulting Hamiltonian can be mapped into a linear chain with variable hoppings. Each site representing an orbital surrounding the impurity with an associated energy (or length) scale. Wilson proposed to solve the linear chain by iterative perturbation. A truncated chain with $`N`$ sites, described by an effective Hamiltonian $`H_N`$, gives the correct physics on an energy scale $`\stackrel{~}{\omega }_N`$. A renormalization group transformation corresponds to adding a site to the chain and relates the Hamiltonians describing successive lower energy scales. This leads to a systematic way of calculating the thermodynamic properties at successive lower temperatures and the spectroscopic properties at successive lower frequencies.

In our case, the Hamiltonian $`H_{AM}^0`$ is rewritten in Wilson’s basis as schematically shown in Fig.1(b). The potential barrier is described including a higher diagonal energy to the orbital centered around $`R_c`$. Due to the structure of Wilson’s basis wavefunctions, a given potential barrier leads to a diagonal energy that decreases as $`R_c`$ increases. Alternatively the barrier can be simulated with a smaller hopping matrix element in the region where $`W(r,R_c)`$ is different from zero. We adopted this last description reducing one hopping term by a factor $`\alpha `$. In what follows it is assumed that the band of extended states is half filled - the Fermi energy is set equal to zero - and $`ϵ_d=U/2`$. This guarantees that the electron-hole symmetry of the problem is preserved.

The properties of a Kondo impurity depend on the local density of states at the impurity coordinate that in Wilson’s representation is the local density of states at site 1. We end this section with a brief discussion on the effect of the confining potential $`W(r,R_c)`$ on the local density of states close to the Fermi energy. For a given $`N`$ the spectrum of $`H_N`$ depends on the parity of $`N`$. In the absence of impurity, a one particle state at zero energy exists only for odd $`N`$ as schematically shown in Fig. 2. The spectrum, as a function of $`N`$, alternates between those of Figs. 2(a) and 2(b). The mean energy separation, between the one-electron states, is given by the characteristic energy $`\mathrm{\Delta }`$ that decreases as $`N`$ increases. This is the spectrum of a finite system ($`N=N_c`$) described by an infinite barrier. The other one electron states, not shown in Fig. 2, are not at energies $`n\mathrm{\Delta }`$ or $`(n+1/2)\mathrm{\Delta }`$ with integer $`n`$. This is due to the logarithmic discretization of the band. How the NRG can be used to analyze finite size effects was shown by NozieresNozieres . To show that the logarithmic discretization correctly describes the low energy spectrum of a impurity in finite systems, in the next section, we compare the exacts results on a finite system with those obtained with the NRG. The alternating spectrum with a fixed Fermi energy corresponds to an alternating parity in the number of electrons, in fact for a fixed electron density, as the radius $`R_c`$ of the cluster increases, the number of electrons alternates between even and odd. We recall that due to the symmetry of the system only the sector of the Hilbert space with s-wave electrons is considered here.

With a finite barrier ($`\alpha 0`$) the discrete states of the cluster are hybridized with the continuum of the host metal and become resonances. The local density of states then becomes a continuum but retains some structure characterized by the energy $`\mathrm{\Delta }`$. In the NRG approach, Hamiltonians $`H_N`$ with $`N>N_c`$ accumulate states at low energies representing the broadening of the central peak for odd $`N_c`$ or the tails of the states at $`\pm \mathrm{\Delta }/2`$ for even $`N_c`$. The local densities of states as obtained with the NRG are also shown in Fig. 2. For $`\alpha 0`$ the number of electrons in the central cluster is no longer a good quantum number and in what follow we refer to the situations of Figs. 2(a) and 2(b) as the Fermi energy being at a resonance (at-resonance) or between two resonances (off-resonance) respectively.

As we show below, this structure of the local density of states determines the thermodynamic and the spectral properties of the impurity.

### II.2 Thermodynamic Properties

Here we briefly recap the thermodynamic properties of the model. As stated above, for $`\alpha =0`$ the at-resonance situation corresponds to an odd number of extended electrons. An impurity in the Kondo limit contributes with an extra electron and the ground state is a singlet indicated as $`|0.`$ The expectation value of the impurity spin $`0|S_{iz}|0`$ is zero reflecting the complete screening of the impurity spin and the zero temperature susceptibility is finite. For the off-resonance situation the ground state of the isolated cluster with a Kondo impurity is a Kramers spin-$`1/2`$ doublet, $`|`$ and $`|`$. The expectation value $`|S_{iz}|=|S_{iz}|`$ is different from zero and, in the low temperature limit the impurity susceptibility diverges as $`\chi =\mu _i^2|S_{iz}|^2/k_BT`$.

For a finite barrier the susceptibility always saturates, however the at-resonance and off-resonance situations give rise to very different temperature dependences. We calculated the impurity magnetic susceptibility given by Desc

$`k_BT\chi `$ $`=`$ $`\mu _i^2({\displaystyle \underset{\nu }{}}P_\nu |\nu |S_{iz}|\nu |^2+`$ (2)
$`2k_BT{\displaystyle \underset{\nu \xi }{}}|\nu |S_{iz}|\xi |^2{\displaystyle \frac{P_\nu }{E_\xi E_\nu }}),`$
where the summation is done over the low energy states $`|\nu `$ with energies $`E_\nu `$ and $`P_\nu =\mathrm{exp}(E_\nu /k_BT)/Z`$. The matrix elements $`\nu |S_{iz}|\xi `$ have to be evaluated in a recursive way at each renormalization step. The susceptibility reflects the thermodynamic properties of the system and we use the effective magnetic moment $`\mu ^2=k_BT\chi `$ as an indicator of the degree of screening of the impurity spin.

In absence of barrier - corresponding to the infinite homogeneous system - the characteristic energy scale is the Kondo temperature indicated as $`k_BT_K^{\mathrm{}}`$. A finite barrier at $`R_c`$ \- enclosing $`N_c`$ shells- introduces the new energy scale $`\mathrm{\Delta }D\mathrm{\Lambda }^{N_c/2}`$, where $`\mathrm{\Lambda }=2`$ is Wilson’s discretization parameter. For $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ the fine structure (on the scale of $`T_K^{\mathrm{}}`$) of the density of states does not change the properties of the system. Conversely, for $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ new confinement induced regimes are observed: for the system at-resonance, as the temperature is lowered, there is a rapid decrease in the magnetic moment $`\mu ^2`$ when $`k_BT\mathrm{\Delta }`$; for the system off-resonance, as $`k_BT`$ approaches $`\mathrm{\Delta }`$ the magnetic moment saturates leading to a plateau in the temperature dependence of $`\mu ^2`$, only at lower temperatures the screening is completed (see Fig. 3). This plateau can be interpreted as the behavior of an isolated cluster. Only at very low temperatures the coupling with the host states becomes relevant and the complete screening can occur. The condition for the existence of these new regimes $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ can be put in the form $`\xi _K^{\mathrm{}}`$ $`R_c`$ where $`\xi _K^{\mathrm{}}=\mathrm{}v_F/k_BT_K^{\mathrm{}}`$ is the Kondo screening length of the infinite system. In other words, only if the Kondo screening length is of the order or larger than the system size, the confinement effects become evident in the thermodynamic properties.

Although the behavior of the system is not universal, at very low temperatures the susceptibility can be used to define an effective Kondo temperature. In fact the low temperature tail of $`k_BT\chi `$ can be scaled to define the effective energy scale $`T_K^{eff}`$ for the low energy spin excitations. We have done this using two different criteria, one is simply to use Wilson’s result:

$$\chi (T=0)=0.103\frac{\mu _i^2}{k_BT_K},$$
(3)
to relate the low temperature susceptibility to an effective $`T_K^{eff}`$. The other alternative is to plot $`\mu ^2`$ vs. $`T/T_K^{eff}`$ fitting $`T_K^{eff}`$ to have a good scaling at low temperatures. The effective Kondo temperatures describing the low temperature behavior are shown in Fig. 4 as a function of the barrier height for the at-resonance and off-resonance situations. For $`\alpha >0.2`$ the two criteria give the same results, for smaller values of $`\alpha `$ the universal behavior is obtained only at extremely low temperatures where numerical errors become important.

## III Spectral Densities of Impurities in Nanoscopic Systems

Using the standard definitions and notation, the impurity Green’s function can be written, using the Lehmann representation, as:

$$G_\sigma (\omega ,T)=\frac{1}{Z(T)}\underset{\lambda ,\lambda ^{}}{}|\lambda |d_\sigma |\lambda ^{}|^2\frac{e^{E_\lambda /k_BT}+e^{E_\lambda ^{}/k_BT}}{\omega (E_\lambda ^{}E_\lambda )},$$
(4)
where $`Z(T)`$ is the grand partition function

$$Z(T)=\underset{\lambda }{}e^{E_\lambda /k_BT}.$$
(5)
The corresponding impurity spectral density is

$`\rho _{i\sigma }(\omega ,T)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}[G_\sigma (\omega +i0^+,T)G_\sigma (\omega i0^+,T)]`$ (6)
$`=`$ $`{\displaystyle \frac{1}{Z(T)}}{\displaystyle \underset{\lambda ,\lambda ^{}}{}}|M_{\lambda ,\lambda ^{}}|^2(e^{E_\lambda /k_BT}+e^{E_\lambda ^{}/k_BT})\times `$
$`\delta (\omega (E_\lambda ^{}E_\lambda )),`$
with $`M_{\lambda ,\lambda ^{}}=\lambda |d_\sigma |\lambda ^{}`$. In the absence of an external magnetic field, the impurity spectral density is spin independent and from here on we drop the spin index in $`\rho _{i\sigma }(\omega ,T)`$.

We now briefly discuss the application of the NRG to evaluate these quantities. The zero temperature limit of expression (6) can be put in the form:

$`\rho _i(\omega ,T=0)`$ $`=`$ $`{\displaystyle \frac{1}{Z(0)}}{\displaystyle \underset{\lambda ,0}{}}|M_{\lambda ,0}|^2\delta (\omega +E_\lambda )`$ (7)
$`+|M_{0,\lambda }|^2\delta (\omega E_\lambda ),`$
where the subindex $`0`$ indicates the ground state, the summation is over all excitations $`\lambda `$ with energies $`E_\lambda `$ and on the components of the ground state in the case of a degenerate ground state, $`Z(0)`$ is the zero temperature partition function that gives the ground state degeneracy. We recall that in the NRG the chemical potential is set at zero and all many body energies $`E_\lambda `$ are measured from the ground state energy.

In order to evaluate the impurity spectral density of the real system using the results of the NRG, consider the spectral density $`\rho _i^N(\omega ,T)`$ corresponding to Hamiltonian $`H_N.`$ As mentioned above, the RG strategy of an iterative diagonalization of a sequence of Hamiltonians $`H_N`$ with $`N=0,1,\mathrm{}`$ is based on the fact that, on an energy scale $`\stackrel{~}{\omega }_N`$, the spectrum of $`H_N`$ is representative of the spectrum of the infinite Hamiltonian $`H`$. Hence, for $`\omega \stackrel{~}{\omega }_N`$ it is a good approximation to takeHofstetter

$$\rho _i(\omega ,T=0)=\rho _i^N(\omega ,T=0).$$

A typical choice is to consider for the $`N^{th}`$ iteration the energy range $`\mathrm{\Lambda }^{N/2}<\omega <2\mathrm{\Lambda }^{N/2}`$ , that is to consider all the states of $`H_N`$ with energies $`E_\lambda ^N`$ in this range and all the corresponding matrix elements $`M_{\lambda ,0}^N`$ to compute the spectral density as

$`\rho _i(\omega ,T`$ $`=`$ $`0)={\displaystyle \frac{1}{Z_N(0)}}{\displaystyle }_{\lambda ,0}|M_{\lambda ,0}^N|^2\delta (\omega +E_\lambda ^N)`$ (8)
$`+|M_{0,\lambda }^N|^2\delta (\omega E_\lambda ^N).`$
Since for the infinite systems the energy spectrum consists of a continuum, at each energy scale $`\stackrel{~}{\omega }_N`$, the discrete spectra are usually smoothed by replacing the delta function $`\delta (\omega E_\lambda ^N)`$ by a smooth distribution $`P_N(\omega E_\lambda ^N)`$. At each $`N`$ we use a Gaussian logarithmic distribution with a width that decreases as the characteristic energy scale of $`H_N`$. It is then clear that at high energies, say the atomic energy of the impurity level $`\epsilon _d`$, the NRG cannot resolve any detailed structure, only at low frequencies there is enough resolution to see, for example, confinement effects. We analyze in detail the impurity spectral density around the Fermi energy which in fact contains all the physics at an energy scale relevant for the Kondo effect and for the electron confinement in a nanoscopic sample.

Before presenting the renormalization group results we show, as a reference calculation, some exact results in finite systems.

### III.1 Exact results in finite systems

Here we present the zero temperature results for an impurity in a finite system consisting of a linear chain of $`N`$ equivalent sites with the impurity at one end. In this case the energy spectrum of extended states $`\epsilon _\nu `$ and the hybridization matrix elements $`V_\nu `$ of Hamiltonian (1) are given by:

$$\epsilon _\nu =2t\mathrm{cos}(\nu \pi /(N+1)),$$

and

$$V_\nu =V_0\mathrm{sin}(\nu \pi /(N+1)),$$

where $`t`$ is the hopping matrix element in the chain and $`\nu =1,2,\mathrm{}.,N`$.

Using a Lanczos algorithm we first calculate the ground state energy and wavevector. In the present case it is convenient to use a second Lanczos algorithmGagliano89  to evaluate directly the impurity spectral function, rather than evaluating excited states, their corresponding energies and matrix elements $`M_{\lambda ,\lambda ^{}}`$. For a better comparison of the different cases, we always measure the frequency from the Fermi energy $`E_F`$ defined as $`E_F=[E_0(N_e+1)E_0(N_e1)]/2`$ with $`E_0(N_e)`$ the ground state energy of the system with $`N_e`$ particles. The spectral density around the Fermi energy shows a series of peaks separated by the characteristic energy $`\mathrm{\Delta }`$. Each peak may be composed by a one or a few delta lines as shown in Figs. 5(a) and 5(b). In what follows we discuss the low energy structure ($`\omega \mathrm{\Delta }`$) of the impurity spectral function for the case of even and odd number of particles corresponding to the at-resonance and off-resonance situations respectively. Note that the results of Figs. 5(a) and 5(b) are not for a half filed system and as a consequence the electron-hole symmetry is not exactly preserved.

Even number of particles: The ground state of the system is a singlet and the impurity spectral density shows a central peak at $`\omega 0`$ and two satellites at $`\omega \pm \mathrm{\Delta }`$. The central peak consists of two lines, one for $`\omega >0`$ and one for $`\omega <0`$ corresponding to adding and subtracting a particle respectively. The splitting between these two lines, that decreases as $`V_0`$ decreases, is not given directly by the confinement energy $`\mathrm{\Delta }`$ but by the energy gained by forming the singlet, i.e. is given by the Kondo temperature of the finite system. The peaks at $`\omega \pm \mathrm{\Delta }`$ are also made of a couple of delta lines each. The justification of this description of a central line and satellites at $`\pm \mathrm{\Delta }`$ is given by the fact that for smaller hybridizations the two central lines approach $`\omega 0`$ and the impurity spectral density reproduces the structure of the unperturbed local density of states (the at-resonance situation of Fig. 2(a)). Here we used these parameters for a better comparison; results with a smaller hybridization calculated with the NRG are shown below.

Odd number of particles: The ground state of the system is a spin-1/2 doublet. The spectral densities have a low frequency gap with peaks at $`\omega \pm \mathrm{\Delta }/2`$ . Each of these peaks at $`\omega \pm \mathrm{\Delta }/2`$ consist of two delta lines that correspond to final states with different total spin: when an electron is created or destroyed on the spin-1/2 ground state, the final state may have total spin zero or one. The energy of the final state depends on its total spin. Again in this case the general structure is that of the underlying local density of states.

These results should serve us as a guide for the numerical renormalization group calculation. Since the NRG involves some approximations and will be extended to more realistic cases, as is the case of a system with a finite barrier, it is important to have this exact result as a reference calculation.

### III.2 Numerical Renormalization Group: Zero Temperature Results

For a finite system, the renormalization procedure is truncated at an iteration $`N_c`$. In Figs. 5(c) and 5(d) we present the NRG results for the impurity spectral density of a finite system with an even and an odd number of particles respectively. As for the calculation of the thermodynamic properties, the system with an even (odd) number of particles is evaluated with an odd (even) number of shells $`N_c`$. The spectrum consists of a discrete collection of delta lines and, as in the previous section, the smoothing is done only for practical purposes in the presentation of the data. Our low frequency NRG results compare very well with the exact results of a finite system \[Figs. 5(a) and 5(b)\]. Again, for an even number of particles the impurity spectral density consists of a central peak, composed by two delta functions, and the two satellites at $`\omega \pm \mathrm{\Delta }`$. For an odd number of particles the approximate NRG results clearly show the central pseudogap and the two satellites at $`\omega \pm \mathrm{\Delta }/2`$ with a structure that comes from the spin-dependent energy of the final state.

For the remainder, we focus on the more interesting case of a system in contact with a macroscopic reservoir. This is described with a non-zero parameter $`\alpha `$ representing a finite wall.

In Fig. 6(a) the impurity spectral density for the at-resonance situation is shown for different values of $`\alpha `$. For $`\alpha =0`$ the isolated cluster results are reproduced with a peak at low frequencies. This peak is well separated from the Fermi energy due to a large hybridization $`V`$ used in the calculation for practical purposes. In the figure only a detail of the low frequency structure is shown and the satellite at $`\omega \mathrm{\Delta }`$ is not observed. As $`\alpha `$ approaches one, corresponding to an infinite homogeneous system, the two structures merge into a single Kondo peak centered at the Fermi energy.

The impurity spectral density when the Fermi energy is off-resonant is shown in Fig. 6(b). For a small $`\alpha `$ the excitations around $`\omega \mathrm{\Delta }/2`$ are clearly observed. As $`\alpha `$ increases the $`\omega \mathrm{\Delta }/2`$ structure is washed out and simultaneously a narrow Kondo resonance develops at the Fermi energy.

For more realistic parameters (a smaller hybridization and consequently a smaller Kondo temperature) the at-resonance and off-resonance impurity spectral densities are shown in Fig. 7. These results show that, for the general case of a Kondo impurity in a nanoscopic system with an intermediate barrier and the Fermi energy at a resonant state, we should expect the low frequency spectral density $`\rho _i(\omega ,T=0)`$ to present a broad Kondo resonance with some structure. In fact, for large barriers (small $`\alpha `$), the Kondo resonance has a minimum at the Fermi energy. For the off-resonant case, we expect the spectral density to present the structures at $`\omega \pm \mathrm{\Delta }/2`$ separated by a pseudogap and - at zero temperature - a central Kondo peak.

As we show in the next section, a small temperature completely destroys this central Kondo peak while the broad structures survive up to much higher temperatures.

### III.3 Phenomenological approach

Here we briefly discuss the structure of the Kondo resonance in terms of what would be obtained with a simple slave boson theory in the saddle point approximation. In this theory ($`U=\mathrm{}`$), the low energy impurity Green’s function is given by Hewson1993

$$G\left(\omega \right)=\frac{b^2}{\omega \epsilon _\lambda b^2V^2g^0\left(\omega \right)}$$
(9)
where $`b^2`$ is the square of the mean value of the boson field, $`\epsilon _\lambda 0`$ gives the position of the Kondo resonance and $`g^0(\omega )`$ is the bare local propagator of conduction electrons. Here as in the NRG calculation we take the Fermi energy equal to zero. In an homogeneous system it is a good approximation to take $`g^0(\omega )=i\pi \rho `$ where $`\rho `$ is the frequency independent local density of states and defining the Kondo temperature $`T_K`$ as the width of the Kondo resonance the local propagator can be put as $`G(\omega )b^2/(\omega +iT_K)`$ with $`b^2=T_K/\pi V^2\rho `$ in agreement with the Friedel’s sum rule. For our case of a central grain weakly coupled with a reservoir, the propagator $`g^0(\omega )`$ is taken as a sum of poles separated by the characteristic energy $`\mathrm{\Delta }`$ and widths $`\gamma `$

$$g^0(\omega )=\underset{l}{}\frac{1}{\omega \mathrm{\Delta }_l+i\gamma }$$
(10)
with $`\mathrm{\Delta }_l=l\mathrm{\Delta }`$ or $`\mathrm{\Delta }_l=(l+1/2)\mathrm{\Delta }`$ for the at-resonance and off-resonance cases respectively. The structure of the Kondo resonance so obtained is shown in Fig 8. The results are in good qualitative agreement with the NRG results. In the at-resonance case, depending on the value of the parameters, the low energy spectrum consists of two peaks with a deep at the Fermi level or a single peak with a maximum at the Fermi level (not shown). In the off-resonance case, a central peak at the Fermi level is obtained. Note however that the area of the central peak relative to the area of the satellite peaks obtained is this approximation is different from that obtained with the NRG.

### III.4 Finite Temperature Results

The NRG calculation of the finite temperature spectral density $`\rho _i(\omega ,T)`$ relies on the same approximations of the $`T=0`$ case. The spectral density at a fixed temperature $`T`$ is evaluated as above - using Eq. (6) instead of Eq. (7) - if $`\omega >k_BT`$. To calculate the spectral density at $`\omega k_BT`$ a Hamiltonian $`H_N`$ with $`N`$ such that $`\stackrel{~}{\omega }_Nk_BT`$ is used.

When $`T_K^{\mathrm{}}\mathrm{\Delta }`$ and the Fermi energy lies between resonances, the off-resonance case, two energy scales are clearly observed in the temperature dependence of the susceptibility. In fact as the temperature is lowered they determine the onset of the plateau in $`k_BT\chi `$, that is given by $`\mathrm{\Delta }`$, and the offset at $`T_K^{eff}`$ that is determined by $`\alpha `$. For the at-resonance case the static magnetic properties are dominated by the larger scale $`\mathrm{\Delta }`$ since for $`k_BT<\mathrm{\Delta }`$ the magnetic moment is completely screened. The spectral densities are also sensitive to the temperature and as we will see, there are changes in $`\rho _i(\omega ,T)`$ each time the temperature approaches each one of the characteristic energy scales.

In Fig. 9(a) the spectral density $`\rho _i(\omega ,T)`$ is shown at different temperatures for a system with the parameters as in Fig. 6. The Kondo peak develops and reaches its zero temperature value as soon as the temperature goes below $`\mathrm{\Delta }`$. The parameters are the same as in Fig. 3 with $`\alpha =0.35`$.

For the off-resonant situation the temperature evolution is shown in Fig. 9(b). In this last case the two characteristic temperatures are also reflected in the magnetic susceptibility from which we can extract their values as the onset and offset of the plateau in $`\mu ^2`$. As can be seen from the figure, as the temperature is lowered, the broad satellite at $`\omega \mathrm{\Delta }/2`$ develops at a temperature corresponding to the onset of the plateau while a narrow Kondo structure develops at a temperature corresponding to the offset of the plateau. The first and second temperatures at which the narrow Kondo peak and the broad structure disappear indicate the decoupling of the impurity spin with the electron spin density at large ($`r>R_c`$) and short ($`r<R_c`$) distances respectively.

## IV Summary and Discussion

We have presented a model for a Kondo impurity in a nanoscopic sample or grain coupled with a reservoir through a large barrier. The natural energy - or length - scale associated to the Kondo effect in a bulk material is to be compared with the new scale introduced by the confinement of electrons. These new scales modify the local density of states at the impurity site and consequently their thermodynamic and spectral properties.

The local density of states of the host material consists of resonant states separated by a characteristic energy $`\mathrm{\Delta }`$ and with a width $`\gamma `$. While $`\mathrm{\Delta }`$ is determined by the size ($`R_c`$) of the grain, the width of the resonances is determined by the transparency ($`\alpha `$) of the barrier. We have considered cases in which the Fermi energy lies at a resonance or between two of them.

The main results of the paper are the frequency and temperature dependence of the impurity spectral density $`\rho _i(\omega ,T)`$. We have performed a detailed analysis of the structure of $`\rho _i(\omega ,T)`$ for $`\omega `$ around the Fermi energy for the case $`T_K^{\mathrm{}}\mathrm{\Delta }`$. This regime with $`T_K^{\mathrm{}}\mathrm{\Delta }`$ is characterized by a interplay between electron-electron correlation and confinement effects. The width of the Kondo resonance ($`T_K^{\mathrm{}}`$) and the characteristic energy ($`\mathrm{\Delta }`$) that defines the structure of the local density of states are of the same order and as a consequence the Kondo resonance shows a superstructure. When the Fermi energy lies at a resonance, the low temperature spectral density consists of a broad structure around the Fermi level. For large barriers, the spectrum presents a two peaked structure with a deep at the Fermi level, where the distance between peaks is given by the Kondo energy of the finite system. As the height of the barrier is decreased, the peaks merge into a single broader peak.

For the off-resonance case $`\rho _i(\omega ,T=0)`$ presents structures at $`\omega \pm \mathrm{\Delta }/2`$ and a central Kondo-like peak at $`\omega 0`$. In temperature dependence two energy scales can be distinguished, at which the two types of structures disappear. These two energy scales $`\mathrm{\Delta }`$ and $`T_K^{eff}`$ are also reflected in the magnetic susceptibility. Phenomenological results based on the slave boson mean field theory for the impurity spectrum at zero temperature give results in qualitative agreement with the ones obtained with NRG.

Different experiments like the STM conductance for impurities in quantum corrals or small terraces or the transport through QD in nanoscopic rings could be used to test the theory.

This work was partially supported by the CONICET and ANPCYT, Grants N. 02151 and 99 3-6343.


text:
# 1 Introduction

## 1 Introduction

For thousands of years mathematicians have enjoyed competing with one another to compute ever more digits of the number $`\pi `$. Among modern physicists, a close analogy is computation of the ground state energy of the Helium atom, begun 75 years ago by E. A. Hylleraas .

Many have contributed incremental steps in this endeavor, flexing their computational muscle and ingenuity, often trying to use mathematical insight for advantage. The strongest line of theory has been to focus on the analytic properties of the wavefunction, especially following the 1954 work of V. A. Fock  which showed a weak logarithmic singularity at the three-particle coalescence.

The recent work of V. I. Korobov  stands out for its simplicity and its success. His trial functions use the three scalar coordinates packaged as $`e^{\alpha r_1\beta r_2\gamma r_{12}}`$, with many sets of the complex nonlinear parameters $`\alpha ,\beta ,\gamma `$ selected in a quasirandom manner within specified intervals $`A_i\alpha _iB_i`$, etc. With a linear combination of N=2200 terms of this type, grouped into four sets $`[A_i,B_i]`$, Korobov surpassed the best previous work by three decimals of accuracy in the Helium ground state energy; and more recently he went on to N=5200 and added four more decimal places.

What struck me as surprising in Korobov’s work was the fact that it seemed to ignore that earlier “wisdom” about analytic properties of wavefunctions. His basis functions are, analytically, no different from the original Hylleraas basis; yet his computational results appear to converge so much more rapidly. This perception motivated the investigations reported below.

Section 2 presents the new experimental data - systematic variational calculations using a variety of different basis functions that have been suggested over the years: these include negative powers, fractional powers and logarithms of the coordinates. Vastly different rates of convergence are observed, as illustrated in Figure 1.

Section 3 offers a qualitative discussion and attempts to interpret this wealth of new data. Section 4 is a review, and an attempt to apply, the theoretical approach for understanding, at least semi-quantitatively, the observed convergence rates. We conclude that this state of the theory is far from satisfactory.

## 2 Experiments - Data

Several different sets of basis functions were used in the standard variational calculations for the ground state energy of the Hamiltonian (with Z=2),

$$H=\frac{1}{2}[_1^2+_2^2]Z/r_1Z/r_2+1/r_{12}$$
(2.1)
and they are detailed below in terms of the Hylleraas coordinates,

$$s=r_1+r_2,t=r_1r_2,u=r_{12}=|\stackrel{}{x}_1\stackrel{}{x}_2|.$$
(2.2)
Table I presents summary results for the primary bases studied. The Accuracy of any energy value E is defined as Log\[E\*/(E\*-E)\] where E\* is our best estimate of the exact value. Thus, the numerical value of Accuracy tells how many correct decimal places there are in the calculated result. Generally, we see that convergence was found to be more rapid as one progressed through this series, A,B,C,D,E,F. More detailed discussion of the results will be deferred to the following sections.

$``$ Basis A is just the original Hylleraas set:

$`\psi ={\displaystyle C_{l,m,n}e^{ks/2}s^lu^mt^n}`$ (2.3a)
$`l,m=0,1,2,3,\mathrm{},n=0,2,4,6,\mathrm{}`$ (2.3b)
and we use the order parameter $`\omega =l+m+n`$ to group the terms. We also designate a calculation of order $`\omega `$ to mean a basis set including all terms with $`l+m+n\omega `$. The total number of terms, N, grows asymptotically as $`\omega ^3/12`$. This same organizational scheme is used for all the experiments listed below.

$``$ Basis B allows negative powers of s, which were introduced by Kinoshita :

$$e^{ks/2}s^l(u/s)^m(t/s)^n;$$
(2.4)
and we can rearrange the terms as

$$e^{ks/2}s^l(u/s)^{m+n}P_n(t/u),$$
(2.5)
using the Legendre polynomials to get the selection rule $`\mathrm{\Delta }n=0,\pm 2`$. This technique follows the work of Goldman  and leads to more efficient use of computer space and time.

$``$ Basis C allows fractional powers of s, which were first introduced by H. M. Schwartz :

$$(1,s^{1/2})e^{ks/2}s^lu^mt^n$$
(2.6)
which doubles the size N of the basis for each order $`\omega `$. (The $`s^{1/2}`$ term is omitted for $`\omega =0`$.)

$``$ Basis D uses the logarithm of s, first introduced by Frankowski and Pekeris :

$$(1,ln(s))e^{ks/2}s^lu^mt^n$$
(2.7)
(The ln(s) term is omitted for $`\omega =0,1`$. The values of N shown in Tables I and II are two units off for bases D and F.)

$``$ Basis E has both negative powers and fractional powers of s.

$$(1,s^{1/2})e^{ks/2}s^l(u/s)^m(t/s)^n.$$
(2.8)
$``$ Basis F has both the logarithm and negative powers of s.

$$(1,ln(s))e^{ks/2}s^l(u/s)^m(t/s)^n.$$
(2.9)
For each of the basis sets described above, the scale parameter k was varied somewhat to find the lowest energy at each order. For the Table I data on Basis A, k varied from 5.0 to 8.2; for Basis B, from 3.8 to 5.9; for Basis C, from 4.7 to 6.6; and for Basis D, from 4.6 to 6.6. For Bases E and F, it was found that the optimum value of k stayed close to 2.0 for the mid-size and larger orders, so k was fixed at this value for all the data shown.

In two cases I replaced the set of functions $`e^{ks/2}s^l`$ by the set $`e^{\alpha _ls}`$, using Korobov’s quasi-random method for selecting the (real) nonlinear parameters $`\alpha `$ in a single group. The results for these experiments are shown as Bases B’ and E’ in Table I.

Several variants of these basis sets were also explored briefly but discarded when they appeared less effective, as functions of N, than their counterparts above. Among these were the use of:

$`\text{Negative powers of s and u }(u/s)^m(t/u)^n`$ (2.10a)
Third roots of s instead of the square root (2.10b)
$`\text{Two or more powers of }ln(s)`$ (2.10c)
$`\text{The coordinate }r=\sqrt{r_1^2+r_2^2}`$ (2.10d)
$`\text{The coordinate }R=|\stackrel{}{x}_1+\stackrel{}{x}_2|\text{ (See the Appendix.)}`$ (2.10e)
| Table I. Energy Accuracy – Summary data | | | | |
| --- | --- | --- | --- | --- |
| $`\omega `$ | N | A-Basis | B-Basis | B’ -Basis |
| 9 | 125 | 7.9 | 8.7 | 9.4 |
| 12 | 252 | 8.7 | 10.2 | 10.7 |
| 16 | 525 | 10.4 | 11.7 | 12.3 |
| 21 | 1078 | 11.6 | 13.2 | 13.9 |
| 27 | 2135 | 12.8 |  |  |
| $`\omega `$ | N | C-Basis | D-Basis |  |
| 7 | 139 | 9.5 | 9.6 |  |
| 9 | 249 | 11.2 | 11.3 |  |
| 12 | 503 | 13.4 | 13.5 |  |
| 16 | 1049 | 15.8 | 16.0 |  |
| 21 | 2155 | 18.3 | 18.5 |  |
| $`\omega `$ | N | E-Basis | E’ -Basis | F-Basis |
| 7 | 139 | 9.5 | 9.8 | 9.4 |
| 9 | 249 | 11.7 | 11.6 | 11.5 |
| 12 | 503 | 14.5 | 14.5 | 14.3 |
| 16 | 1049 | 18.2 | 18.3 | 18.1 |
| 21 | 2155 | 22.5 |  | 22.4 |
| 27 | 4269 | 27.2 |  | 27.6 |
| 34 | 8093 | 29.7 |  | 33.0 |

Accuracy = Number of correct decimals

In Table II are the detailed results (in a compact format) for the two bases - E and F - that show the most rapid convergence. One quantitative measure of the rate of convergence is the Ratio of successive differences and this is also shown in the table.

Technical notes. For these computations, I wrote a set of subroutines for multiple-precision arithmetic (in C), eschewing more professional packages which are available. The last row of data in Table II used 101 decimals of precision and took one week running on an otherwise idle desktop computer equipped with a 300 MHz processor and 320 MB of memory.

| Table II. Calculated results with the two best bases | | | | | |
| --- | --- | --- | --- | --- | --- |
| $`\omega `$ | N | E-Basis Energies | Ratios | F-Basis Energies | Ratios |
| 4 | 43 | -2.90372 183 |  | -2.90371 941 |  |
| 5 | 67 | -2.90372 42300 | 17.9 | -2.90372 415 | 22.5 |
| 6 | 99 | ” 43643 | 11.3 | -2.90372 43610 | 14.2 |
| 7 | 139 | 43762 2 | 15.8 | ” 43758 7 | 14.0 |
| 8 | 189 | 43769 759 | 14.3 | 43769 382 | 12.2 |
| 9 | 249 | 43770 28348 | 10.5 | 43770 25283 | 11.2 |
| 10 | 321 | 43770 33352 | 7.32 | 43770 33068 | 8.35 |
| 11 | 405 | ” 34036 | 9.37 | ” 34000 4 | 8.80 |
| 12 | 503 | 34109 6 | 8.30 | 34106 294 | 8.94 |
| 13 | 615 | 34118 444 | 8.63 | 34118 13633 | 9.02 |
| 14 | 743 | 34119 46685 | 8.79 | 34119 44847 | 10.1 |
| 15 | 887 | 34119 58229 | 8.06 | 34119 57846 | 7.33 |
| 16 | 1049 | ” 59667 | 10.4 | ” 59620 | 10.0 |
| 17 | 1229 | 59806 | 6.32 | 59797 | 5.90 |
| 18 | 1429 | 59828 5 | 10.5 | 59827 3 | 9.38 |
| 19 | 1649 | 59830 654 | 5.12 | 59830 456 | 5.50 |
| 20 | 1891 | 59831 06419 | 9.69 | 59831 03831 | 8.89 |
| 21 | 2155 | 59831 10650 | 5.08 | 59831 10381 | 6.18 |
| 22 | 2443 | ” 11482 | 9.14 | ” 11442 | 8.22 |
| 23 | 2755 | 11573 8 | 5.62 | 11571 6 | 7.25 |
| 24 | 3093 | 11589 997 | 7.94 | 11589 408 | 6.66 |
| 25 | 3457 | 11592 03902 | 6.18 | 11592 08081 | 8.87 |
| 26 | 3849 | 11592 36947 | 5.29 | 11592 38154 | 4.79 |
| 27 | 4269 | ” 43186 | 4.80 | ” 44444 | 11.0 |
| 28 | 4719 | 44484 | 2.97 | 45017 4 | 3.57 |
| 29 | 5199 | 44922 | 2.88 | 45177 752 | 13.1 |
| 30 | 5711 | 45074 35 | 2.32 | 45189 95689 | 3.00 |
| 31 | 6255 | 45139 97 | 2.29 | 45194 02040 | 14.8 |
| 32 | 6833 | 45168 66 | 2.16 | 45194 29566 | 2.76 |
| 33 | 7445 | 45181 95 | 2.11 | ” 39521 | 15.1 |
| 34 | 8093 | 45188 24 | 2.06 | 40182 | 2.78 |
| 35 | 8777 | 45191 29 | 2.02 | 40420 7 | 13.5 |
| 36 | 9499 | 45192 80 |  | 40438 342 | 3.11 |
| 37 | 10259 |  |  | 40444 00495 |  |

extrapolate: E\* = -2.90372 43770 34119 59831 11592 45194 40444 6

## 3 Experiments - Discussion

Figure 1 provides a visual comparison of the convergence rates for the different variational basis sets (A - F), plotting Accuracy versus the Log of N, the number of basis functions used. I found a number of surprises in these results.

Surprise 1. Basis B (negative powers of s) shows a significant improvement over basis A (the original Hylleraas basis).

Surprise 2. Bases E and F do a great deal better than any of the others. Somehow, the benefits of B and C (or B and D) are cumulative.

Surprise 3. The performance of C and D are nearly identical, as with E and F (until we reach very high orders.)

Surprise 4. The performance of basis E drops off dramatically after $`\omega =27`$; but basis F keeps up its rapid convergence, although with marked oscillation, as seen from the Ratios in Table II.

The surprising performance of Korobov’s basis has already been noted: his published results are shown by the line labelled with the letter “K” in Figure 1.

The fact that Basis C performs a lot better than Basis A was not a surprise, since earlier work  had already shown that. Ditto for Basis D . Also, the smallness of the gain shown by basis B’ over B (and by E’ over E) is as expected, based upon the analytic equivalence of exponentials and power series.

While I cannot explain the surprising results, I can readily offer suggestions on how one might interpret them. The virtue of Bases C and D lies in providing more flexibility to the “radial” behavior of the wavefunction (coordinate s); while that of the negative powers lies in providing more flexibility in the “angular” behavior (ratios u/s and t/s). The similarity between C and D (and between E and F) indicates that the precise analytic behavior in the “radial” coordinate is not important - any flexibility will do - until one gets to the very high orders.

This appeal to “flexibility” is just armwaving; it lacks any mathematical foundation. Such appeal to flexibility is also the best way I know to understand the success of Korobov’s calculations: his work seems akin to the “adaptive” techniques used in numerical integration, where one puts additional mesh points into any region that shows a slower rate of convergence.

In varying the scale parameter k, I most always found a simple minimum in the energy. However, in one case - Basis E at $`\omega =21`$ \- a more complex behavior emerged: see Table III. While the variation shown here is not very great, this does raise the general question of how effectively one may search for the minimum of a complicated function of many nonlinear parameters. This is a possible source of worry in using Korobov’s technique, especially when it comes to error estimation and extrapolation.

| Table III. Double minimum in E(k) for Basis E at $`\omega `$ = 21 | | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| k | 1.7 | 1.8 | 1.9 | 2.0 | 2.1 | 2.2 | 2.3 | 2.4 |
| E | 10567 | 10697 | 10691 | 10650 | 10635 | 10640 | 10636 | 10603 |
| These numbers for the energy E follow the first 20 decimal places. | | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |

## 4 Theories

The first lesson in analysis of atomic wavefunctions concerns the two-particle cusps: linear behavior as any one of the coordinates $`r_1,r_2`$, or $`r_{12}`$ goes to zero. All of the basis functions studied in this paper are correct in that regard; we are concerned here with what comes next.

Take the Hylleraas expansion (2.3a) and put it into the Schrodinger equation (2.1). Then collect the coefficients of each monomial in s, u, and t and set that equal to zero. Early in this infinite set of algebraic equations for the expansion coefficients $`C_{l,m,n}`$ one finds the following *inconsistent* equations

$`C_{1,0,0}+ZC_{0,0,0}=0`$ (4.1a)
$`2C_{1,1,0}+C_{1,0,0}=0`$ (4.1b)
$`4C_{1,1,0}C_{1,0,0}=0.`$ (4.1c)
This contradiction in the Hylleraas basis was discovered by Bartlett, Gibbons and Dunn  in 1935 and it led them to consider an alternative expansion: one that involved logarithms of the hyperradius $`r=\sqrt{r_1^2+r_2^2}`$. Later, Fock  independently developed a systematic expansion of the wavefunction with such terms.

In 1962 this author  developed a general theory about the convergence rate of variational calculations, based upon analogy with least-squares fitting of functions and one-dimensional model problems. This line of analysis was expanded by others  and in some cases given a more rigorous mathematical basis.

That theoretical work led directly to the idea that the convergence rate in Hylleraas-type calculations for the Helium ground state was controlled by the Fock logarithmic singularity; and the semi-quantitative analysis seemed to fit the available data. It also led to the successful exploitation of the fractional power basis C . Shortly thereafter, Frankowski and Pekeris  took logarithmic terms explicitly into their trial functions and this also seemed to confirm the importance of the Fock behavior.

However, when Kinoshita  considered the systematic use of negative powers - in the form (2.10a), not (2.4) - he found that there was no contradiction of the type noted above. I have confirmed that this conclusion holds also for the basis (2.4) used in the current work.

Thus, one might be drawn to believe that the logarithmic singularity is not an intrinsic property of the He atom wavefunction itself but rather results from a particular analysis that relies on the six-dimensional coordinate r. This idea may be dispelled by reading the work of Myers et al , where they take one for a walk around the three-particle coalescence and show that the (finite) discontinuity in the local energy disappears when one includes the full set of terms that accompany the Fock logarithm.

This approach helps us understand some other experimental results. We noted earlier that inclusion of negative powers $`t/u`$ gave poorer results than $`t/s`$. Each of these ratios shows a (finite) discontinuity when one walks around the place where the denominator vanishes - something that the correct wavefunction should not allow. In the case of (t/u) this error occurs along a line, while in the case of (t/s) this error occurs only at a point. A similar situation appears in the work of Goldman , who holds the previous record for basis size (N=8066). His use of the coordinates $`r_<,r_>`$ led to a very efficient computational scheme, but demonstrates relatively slow convergence. (See the point “G” in Figure 1.) This may be attributed to the discontinuity of his basis functions on the 2-dimensional surface $`r_1=r_2`$.

### 4.1 Fitting the Data

The goal of a good theoretical understanding should be the ability to predict or to explain, at least semi-quantitatively, the observed rates of convergence for systematic variational calculations with different basis sets. In my earlier work , the attempt to do this was based upon analogies with one-dimensional model problems, doing least-squares fit with appropriate orthogonal bases to represent functions with various types of singularities:

$`Minimize`$ $`{\displaystyle \rho (x)𝑑x[f(x)\underset{i=0}{\overset{n1}{}}C_iu_i(x)]^2}`$ (4.2c)
$`C_i={\displaystyle \rho (x)𝑑xf(x)u_i(x)}`$
$`Error(C_n)^2.`$
For one example we find:

$`f(x)=x^\nu lnx,`$ $`\rho (x)=x^\mu \text{on the interval (0,1)}`$
$`C_n`$ $``$ $`1/n^{\mu +2\nu +3/2};`$ (4.3)
and an alternative example is:

$`f(x)=x^\nu lnx,`$ $`\rho (x)=x^\mu e^x\text{on the interval (0,}\mathrm{}\text{)}`$
$`C_n`$ $``$ $`1/n^{\mu /2+\nu +1}.`$ (4.4)
The difference in convergence rates for these two examples may be understood qualitatively as follows. The basis functions $`x^ne^x`$ peak at $`x=n`$. Therefore, at higher n these basis functions on the interval $`(0,\mathrm{})`$ get farther and farther away from the singularity, which is at $`x=0`$. One may improve the situation by using basis functions $`x^ne^{kx}`$, where k is a scale parameter that may grow as one proceeds to higher orders. I do not have a quantitative theory for this result but it is qualitatively relevant to the current study.

In my 1962 work I applied this simple modeling to the He atom problem, identifying the Fock term $`r^2lnr`$ as the dominant singularity which is neglected in conventional Hylleraas coordinates. This led me to predict a convergence rate formula,

$$E(\omega )E(\omega 1)const./\omega ^p$$
(4.5)
and I estimated that p should be between 5.5 and 10, due to uncertainties in replacing the real 3-dimensional problem with the one-dimensional model. The then best results with Hylleraas variables (work of Pekeris, 1959, up to order 21, using a cleverly orthogonalized basis) fit the convergence rate formula (4.5) with a value of p between 7 and 8. This was good confirmation of the theory. The extended computations reported here (Basis A data in Table I) fit the convergence rate formula (4.5) with a value of p which varies from 7, at the lower orders, to a value about 12 at the higher orders. This improvement is probably due to my allowing the scale parameter k to vary, which was not done in the earlier work.

Also, in 1962, I introduced the half-powers of coordinate s, explicitly for the purpose of increasing the convergence rate, following this theory. That was successful, with the observed value of p approximately doubled to 14 or 15 at $`\omega 8`$. The extended computations reported here (Basis C data in Table I) are fit to values of p which vary from about 16 to 21. Again, this is fairly good confirmation of the theory; and again we acknowledge some improvement by allowing the scale parameter to vary.

Following that earlier theory one would certainly not expect Basis D to converge at the same rate as Basis C – but this is exactly the behavior we have found in the present experiments.

What can I say about the observed convergence rate of Basis B, introducing negative powers into the Hylleraas functions? The data in Table I are fit with a value of the exponent p around 13. I do not understand this but will only offer a guess that it may have to do with fitting the complex “angular” behavior around the Fock singularity, which was described in Ref. . Maybe this is connected with the difference in convergence rates noted above, in Eqs. (4.3, 4.4), for the model problems on (0,1) and on (0,$`\mathrm{}`$ ).

Finally, look at the results for Bases E and F. The data in Tables I and II are fit with values of the exponent p which grow from the 20’s to the 40’s in the middle range of $`\omega `$; at the top end, the data for Basis E drop to around p=25, while the data for Basis F climb to about p = 65. I am at a loss to explain these large exponents following the former analysis.

An alternative to the power law convergence rate formula (4.5) is the exponential rate formula

$$E(\omega )E(\omega 1)const.(\sigma )^\omega $$
(4.6)
which one could expect from a model fitting problem that involved no singularities at all. For example, expanding $`e^{ax}`$ in terms of $`x^ne^{bx}`$ would yield the formula (4.6) with $`\sigma =(\frac{ab}{a+b})^2`$. If one plots the data for Basis F (Log of increments in E vs. $`\omega `$), it does look close to a straight line; and the smoothed data in Table II may be fitted with a value of $`\sigma `$ in the range 0.13 - 0.16 for $`\omega >16`$. If one looks at the asymptotic behavior of the He wavefunction as $`r_1`$ goes to infinity, the behavior in $`r_2`$ should be as $`e^{Zr_2}`$ with $`Z=2`$. The trial functions I used for this basis have the exponential envelope $`e^{k(r_1+r_2)/2}`$ with $`k=2.`$ Using the formula quoted above, this model gives us the parameter $`\sigma `$ as $`(\frac{21}{2+1})^2=0.11`$. This looks like a fairly good fit to the data; but accepting this explanation would lead us to doubt the relevance of the Fock singularities for the He wavefunction.

Also, I know of no published theoretical attempts to explain the excellent convergence found by Korobov with his highly nonlinear fitting of the trial wavefunction. John Morgan has suggested (in private communication) that Korobov’s approach may be likened to the work of fitting the Hydrogen radial wavefunction with a set of gaussians, using “floating exponents”. This sounds plausible, but at present it is just more handwaving about “flexibility.”.

I conclude that theoretical understanding of the convergence of variational calculations on the two-electron atom is far outstripped by the raw computing power of available machinery.

Some may ask if any of this is really relevant to current issues in physics. One response is to point to high accuracy measurements performed on atomic systems which may check the current theories of fundamental particles and interactions. A recent paper by Pachucki and Sapirstein  aims to determine the fine structure constant to a few parts-per-billion. This is based upon measurements of the $`2^3P_J`$ states in Helium and detailed calculations that rely upon a Korobov-type representation of the atomic wavefunction.

Then, again, all this may be nothing more than an expression of $`\pi `$-envy.

ACKNOWLEDGEMENT

I am grateful to John D. Morgan III for several very helpful discussions.

Appendix: Integrals

Integrals of the following type were needed in the calculations reported here:

$$\mathbf{}_\mathrm{𝟎}^{\mathbf{}}𝒅𝒔𝒆^\mathbf{}𝒔𝒔^𝒑\mathbf{(}𝒍𝒏\mathbf{(}𝒔\mathbf{)}\mathbf{)}^𝒒\mathbf{.}$$
(A.1)
There is a simple recursion formula on the index p; and for the minimum values of p I used a particular technique of numerical integration. (See Ref. .) First change variables, $`𝒔\mathbf{=}𝒆𝒙𝒑\mathbf{(}𝒚\mathbf{)}`$; then use the simple rule,

$$\mathbf{}_{\mathbf{}\mathbf{}}^{\mathbf{}}𝒇\mathbf{(}𝒚\mathbf{)}𝒅𝒚\mathbf{}𝒉\underset{𝒏\mathbf{=}\mathbf{}\mathbf{}}{\overset{\mathbf{}}{\mathbf{}}}𝒇\mathbf{(}𝒏𝒉\mathbf{)}\mathbf{.}$$
(A.2)
The summation is truncated when terms are smaller than the desired accuracy; and the answer converges exponentially as the interval $`𝒉`$ is decreased.

For the two-electron atom, one can evaluate the most conventional integrals from the formula,

$$\mathbf{}\frac{𝒅^\mathrm{𝟑}𝒙_\mathrm{𝟏}}{\mathrm{𝟒}𝝅}\mathbf{}\frac{𝒅^\mathrm{𝟑}𝒙_\mathrm{𝟐}}{\mathrm{𝟒}𝝅}\frac{𝒆^{\mathbf{}𝒂𝒓_\mathrm{𝟏}}}{𝒓_\mathrm{𝟏}}\frac{𝒆^{\mathbf{}𝒃𝒓_\mathrm{𝟐}}}{𝒓_\mathrm{𝟐}}\frac{𝒆^{\mathbf{}𝒄𝒓_{\mathrm{𝟏𝟐}}}}{𝒓_{\mathrm{𝟏𝟐}}}\mathbf{=}\frac{\mathrm{𝟏}}{\mathbf{(}𝒂\mathbf{+}𝒃\mathbf{)}\mathbf{(}𝒃\mathbf{+}𝒄\mathbf{)}\mathbf{(}𝒄\mathbf{+}𝒂\mathbf{)}}\mathbf{,}$$
(A.3)
and derivatives of this simple result with respect to the parameters a,b,c.

In exploring more complicated functions, I was able to find another simple formula for the following integral, which involves $`𝑹\mathbf{=}\mathbf{|}\stackrel{\mathbf{}}{𝒙}_\mathrm{𝟏}\mathbf{+}\stackrel{\mathbf{}}{𝒙}_\mathrm{𝟐}\mathbf{|}`$.

$`{\displaystyle \mathbf{}\frac{𝒅^\mathrm{𝟑}𝒙_\mathrm{𝟏}}{\mathrm{𝟒}𝝅}\mathbf{}\frac{𝒅^\mathrm{𝟑}𝒙_\mathrm{𝟐}}{\mathrm{𝟒}𝝅}}`$ $`{\displaystyle \frac{𝒆^{\mathbf{}𝒂𝒓_\mathrm{𝟏}}}{𝒓_\mathrm{𝟏}}}{\displaystyle \frac{𝒆^{\mathbf{}𝒃𝒓_\mathrm{𝟐}}}{𝒓_\mathrm{𝟐}}}{\displaystyle \frac{𝒆^{\mathbf{}𝒄𝒓_{\mathrm{𝟏𝟐}}}}{𝒓_{\mathrm{𝟏𝟐}}}}{\displaystyle \frac{𝒆^{\mathbf{}𝒅𝑹}}{𝑹}}`$ (A.4)
$`\mathbf{=}`$ $`{\displaystyle \frac{\mathrm{𝟏}}{\mathbf{(}𝒂^\mathrm{𝟐}\mathbf{+}𝒃^\mathrm{𝟐}\mathbf{}\mathrm{𝟐}𝒄^\mathrm{𝟐}\mathbf{}\mathrm{𝟐}𝒅^\mathrm{𝟐}\mathbf{)}}}𝒍𝒏{\displaystyle \frac{\mathbf{(}𝒂\mathbf{+}𝒃\mathbf{+}\mathrm{𝟐}𝒄\mathbf{)}\mathbf{(}𝒂\mathbf{+}𝒃\mathbf{+}\mathrm{𝟐}𝒅\mathbf{)}}{\mathrm{𝟐}\mathbf{(}𝒂\mathbf{+}𝒄\mathbf{+}𝒅\mathbf{)}\mathbf{(}𝒃\mathbf{+}𝒄\mathbf{+}𝒅\mathbf{)}}}\mathbf{.}`$
To derive this, insert the Laplacian operators into the middle of the integral and let them work both ways. It appears that one could almost deduce these results (A.3) and (A.4) purely by arguments of analyticity and symmetry. Consider, for example, how the integral behaves as $`𝒓_\mathrm{𝟏}\mathbf{}\mathbf{}`$: by counting powers one sees the nature of the singularity as (a+c), or (a+c+d), goes to zero.

As noted earlier, using this variable R in the He trial wavefunction did not produce good results - as one might expect since it introduces a spurious cusp when the two electrons are on opposite sides of the nucleus. I have, nevertheless, recorded the above information here in case it might be useful to others.

The result (A.4) can be generalized with $`𝑹\mathbf{=}\mathbf{|}𝝁\stackrel{\mathbf{}}{𝒙}_\mathrm{𝟏}\mathbf{+}𝝂\stackrel{\mathbf{}}{𝒙}_\mathrm{𝟐}\mathbf{|}`$.