text:
# Quasi-Particle Spectra, Charge-Density-Wave, Superconductivity and Electron-Phonon Coupling in 2H-NbSe2

Permanent address: \]Advanced Light Source, LBNL, Berkeley, CA 94720

## Abstract

High-resolution photoemission has been used to study the electronic structure of the charge density wave and superconducting dichalcogenide, 2H- NbSe2. From the extracted self-energies, important components of the quasiparticle interactions have been identified. In contrast to previously studied TaSe2, the CDW transition does not affect the electronic properties significantly. The electron-phonon coupling is identified as a dominant contribution to the QP self-energy and is shown to be very anisotropic ($`k`$\- dependent) and much stronger than in TaSe2.

The family of layered 2H dichalcogenides represents an interesting system in which charge-density wave (CDW) order co-exists with the superconductivity (SC) Wilson ; Castro . The fact that the CDW transition temperature decreases, while the superconducting critical temperature ($`T_C`$) increases from TaSe2 through TaS2 and NbSe2 to NbS2 suggests that these two order parameters represent competing ground states. Indeed, it has been found that in TaS2 and NbSe2, $`T_C`$ increases under pressure while $`T_{CDW}`$ decreases Smith ; Moline . After CDW order disappears, $`T_C`$ remains approximately constant. In NbS2, the system without CDW order, $`T_C`$ is insensitive to pressure. Although various anomalies, including an apparent anisotropy of the superconducting gap, have been observed Garoche ; Graebner ; Corcoran ; Yokoya , it is generally believed that superconductivity in the dichalcogenides is of conventional BCS character, mediated by strong electron-phonon coupling Corcoran . However, consensus on the exact mechanism that drives the system into the CDW state has still not been reached. Some authors Wilson ; Doran ; Wilson1 ; Castro  argue, in analogy with a Pierls transition in one- dimensional systems, that the CDW transition is driven by a FS instability (nesting), where some portions of the FS are spanned by a CDW vector- $`q_{CDW}`$. In another scenario, the CDW instability is induced by the nesting of van Hove singularities (saddle points in the band structure with high density of states) if they are within a few $`k_BT_{CDW}`$ of the Fermi energy Rice . The Fermi surface of the 2H dichalcogenides is rather complicated, being dominated by several open (2D-like) sheets and with one small 3D S(Se)-derived pancake-like FS Corcoran . In such a situation, one may anticipate anisotropic properties and in particular, an anisotropic electron-phonon coupling. The resistivity anisotropy, of the order of 10-50, is much smaller than in layered oxides, indicating a substantial inter-layer hopping Laszlo . Transport properties show relatively small anomalies at $`T_{CDW}`$, suggesting that only a small portion of the FS becomes gapped in the CDW state. In addition, the 2H dichalcogenides become better conducting in the CDW state, indicating a higher degree of coherence.

We have previously studied the quasi-particle (QP) self-energy for the first member of the series, TaSe2, and found that it was strongly influenced by the CDW transition Valla-TaSe . Sharp structure in both the dispersion (”kink”) and the scattering rate at $``$70 meV in the CDW state weakened and shifted to $``$35 meV upon transition into the normal state. Similar ”kinks”, have been previously identified in conventional metals and attributed to electron-phonon coupling Mo-Be . As the energy scale of the ”kink” in the CDW state of TaSe2 was too large for phonons, we suggested that the excitation was electronic in origin: a fluctuation of the CDW order parameter. In this letter, we show that the self-energy of NbSe2 ($`T_{CDW}35`$ K and $`T_C=7.2`$ K) is less sensitive to the CDW transition and that it is dominated by electron-phonon coupling. Measured directly for the first time, the coupling is very anisotropic, with the largest value on the inner K-H centered sheet.

The experiments reported here were carried out on a high-resolution photoemission facility based on undulator beam line U13UB at the National Synchrotron Light Source and a Scienta SES-200 electron spectrometer, which simultaneously collects a large angular window $`14^{}`$ of the emitted photoelectrons. The combined instrumental energy resolution was set to $``$4 meV for low temperature studies ($`T<15`$K) and to $``$6 meV elsewhere. The angular resolution was better than $`\pm 0.1^{}`$ translating into a momentum resolution of $`\pm 0.0025\AA ^1`$ at the 15.2 eV photon energy used in the present study. Samples, grown by the standard iodine vapor transport method, were mounted on a liquid He cryostat and cleaved in-situ in the UHV chamber with base pressure $`3\times 10^9`$ Pa. The temperature was measured using a calibrated silicon sensor mounted near the sample.

Figure 1 shows the photoemission intensity, recorded at $`T=`$10 K, as a function of binding energy and momentum along three different momentum lines in the Brillouin zone. Nine Fermi crossings are included: three pairs on the double-layer split Fermi sheets centered around $`\mathrm{\Gamma }`$ and three crossings on the split sheets centered at the K point. A characteristic change in the QP velocity (”kink”) can be easily identified in all crossings. The kinks are also accompanied by a sharp change in the QP widths at the ”kink” energy. These observations are indicative of (bosonic) excitations interacting with the QPs. The ”kink” occurs roughly at the energy of the excitation involved in the coupling. Compared to TaSe2, the excitation spectrum is limited to significantly lower energies. It also appears that the kink is not unique; its strength and energy depend on k, being different for different crossings. To quantify this, we have used momentum distribution curves (MDC) to extract the band dispersions for all the states shown. The dispersion curves provide direct information on the real part of the self- energy, Re$`\mathrm{\Sigma }`$. As shown in the figure, the high-energy part of the extracted dispersion can be fitted with a parabola that crosses through $`k_F`$, whereas the low energy part ($`\omega <1520`$ meV) is fitted with a straight line. If we assume that the parabola represents the ”non-interacting” dispersion (i.e. the dispersion in the absence of interactions that cause the kink), then the slopes of these two lines at $`\omega =0`$ may be used to directly extract the coupling constant, $`\lambda =v_F^0/v_F1`$, where $`v_F^0`$ is the ”non-interacting” (bare) Fermi velocity and $`v_F`$ is the renormalized one.

The ”non-interacting” parabolas are further subtracted from the measured dispersions to extract Re$`\mathrm{\Sigma }(\omega )`$. The results are shown in Fig. 2(a) for several crossings from Fig. 1. Re$`\mathrm{\Sigma }(\omega )`$ gives the same coupling constant ($`\lambda =[(\text{Re}\mathrm{\Sigma })/\omega ]_0`$), but also provides additional information about the spectrum of excitations interacting with the QPs. It is obvious from Fig. 2(a) that not only is the magnitude of Re$`\mathrm{\Sigma }`$ different for different states, but also the peaks are at different energies, ranging from $``$13 meV to $``$35 meV. Various experimental and theoretical studies have shown that the phonon spectrum is fully consistent with these energies, with acoustical phonon branches laying below $`\omega `$12 meV, and optical branches spanning the region $`15<\omega <40`$ meV phonons . Shifts of the Re$`\mathrm{\Sigma }`$ maxima would further suggest that some electronic states are coupled predominantly to acoustic modes while others couple more strongly to the optical modes, even though the states are sometimes very close in momentum (compare points 4 and 5, for example). A strong $`k`$-dependence of $`\mathrm{\Sigma }`$ would complicate the MDC line-shape in the energy region where the $`k`$-dependence exists. Although we have detected some deviation from Lorentzian lineshapes, we were not able to precisely determine the MDC’s line-shape due to the overlap of bi-layer split states. It is interesting that in spite of these differences in $`\mathrm{\Sigma }`$, the resulting coupling constant does not vary much, $`\lambda 0.85\pm 0.15`$, within the experimental uncertainty. The only exception is the inner K-centered sheet (point 6), where $`\lambda 1.9\pm 0.2`$. We have completed several measurements on different samples, always with similar results for $`\lambda `$ at, or close to that Fermi point. This seemingly too large coupling constant is actually in a good agreement with the large measured value of linear specific heat coefficient, $`\gamma 18.5`$ mJmol-1K-2 Garoche ; specific-heat , which is proportional to the renormalized density of states (DOS) at the Fermi level, $`N(0)(1+\lambda )`$, through $`\gamma =(1/3)\pi ^2k_B^2N(0)(1+\lambda )`$. Band structure calculations give the ”bare” DOS $`N(0)2.8`$ states eV-1 unit cell-1 Corcoran , suggesting $`\lambda 1.8`$. However even this might be an underestimate for our state as $`\gamma `$ measures an average over the FS, weighted by each state’s DOS. A similar value for $`\lambda `$ is obtained from $`c`$-axis optical conductivity Sasa  suggesting that the $`c`$-axis transport is probably dominated by the K-H centered cylinders with largest warping.

It is instructive that in TaSe2 the CDW gap opens up in the same region of the FS ARPES-Kpoint , while the $`\mathrm{\Gamma }`$-A centered Fermi cylinders remain ungapped, and gain coherence in the CDW state Valla-TaSe . Therefore, it seems plausible that both the superconductivity and the CDW state originate from the inner K sheet and are driven by strong electron-phonon coupling. This seems to be in line with the original suggestion of Wilson Wilson1  that the self-nesting of the inner K sheet drives the CDW in 2H dichalcogenides. A lack of CDW gap on the $`\mathrm{\Gamma }`$ centered sheets in all of 2H dichalcogenides studied in ARPES suggests that these sheets support neither the self-nesting nor the nesting which would mix them with the K-centered sheets. In particular, the f-wave symmetry of the CDW gap Castro  may be ruled out. The relative strength of the CDW and superconducting ordering is determined by the nesting properties of the inner K cylinder, while the upper limit for $`T_C`$ (when the CDW is destroyed by applying pressure, for example) is given by $`\lambda `$. Nesting weakens with increasing 3D character (increased warping with $`k_Z`$) under pressure and on moving from TaSe2 to NbS2. $`\lambda `$ increases from TaSe2 to NbSe2 NbS2  and is essentially pressure independent. In agreement with previous ARPES studies Straub ; Yokoya ; Tonjes , we do not see a CDW gap in NbSe2, suggesting that the nested portion of the FS is very small and was not sampled in any study. As there is a non-trivial $`k_Z`$-dispersion (warping), it is possible that the in-plane $`k_F`$ might be tuned into the nesting and that the gap opens only near certain $`k_Z`$. Note that the energy splitting between the double walled sheets is larger for K-centered sheets. A similar $`k`$-dependence is also expected for the interlayer hopping, $`t_{}`$, that produces the warping. Additionally, as Fermi velocities are larger for $`\mathrm{\Gamma }`$\- centered sheets, it is reasonable to expect that the in-plane $`k_F`$ varies with $`k_Z`$ much less on the $`\mathrm{\Gamma }`$-cylinders than on the K-cylinders (the change in the in-plane Fermi momentum is approximately given by $`\mathrm{\Delta }k_Ft_{}/v_F`$). The measured FSs centered at $`\mathrm{\Gamma }`$ are too large at the sampled $`k_Z`$, and we do not expect them to ever reach the self-nesting condition $`2k_F=q_{CDW}`$. On the other hand, the inner K-centered sheet seems to be very close to producing the required nesting. A more detailed mapping is needed to explore the nesting properties and eventual CDW gap opening in this region. According to STM studies STM , the CDW gap is large ($`\mathrm{\Delta }_{CDW}35`$ meV) and should be easily measurable in ARPES. The overall electronic properties in NbSe2 are much less sensitive to the CDW transition than in TaSe2. Even the CDW induced structure in the self-energy that existed in TaSe2 is absent in NbSe2. Both the ”kink” and the scattering rate are remarkably insensitive to the CDW (See Fig. 2b).

An interesting question is whether the anisotropic electron-phonon coupling constant $`\lambda `$ would be projected into the magnitude of the superconducting gap. A recent photoemission study Yokoya  has shown that the superconducting gap is indeed anisotropic, being quite uniform, $`\mathrm{\Delta }1`$ meV, on the Nb-derived Fermi cylinders, but reduced beyond the experimental sensitivity ($`\mathrm{\Delta }0`$) on the 3D pancake-like FS. No data for the inner K-H cylinder has been reported. In Fig. 3 we show the spectral intensity at two Fermi points, 6 and 7 for several temperatures. The coupling constant differs by a factor of $`2`$ at these two points and yet, the superconducting gap, that clearly opens up below $`T_C`$ at both Fermi points, is the same within the experimental error. This is in contrast to MgB2, where the ”hot” regions are gapped by correspondingly larger gaps Tu . The equalizing of the superconducting gap for states with different coupling strengths but similar symmetries represents a $`k`$-space analogous of the proximity effect proximity .

NbSe2 has several properties in common with other layered materials. One of the most obvious peculiarities is the fact that the transport becomes more metallic in the CDW state even though a portion of the FS is ”destroyed” (gapped). We suggest that this is a manifestation of anisotropic Fermi surfaces where strongly coupled portions, or ”hot spots”, generally play a negative role in normal state conductivities, acting as scattering ”sinks” for the remaining, less renormalized regions. These ”hot spots” drive the system into an ordered state, but it is the ”cold spots” that usually dominate the conductivities. This duality has recently been detected in MgB2 where the normal state transport shows an extremely weak coupling, while the superconducting properties are dictated by strongly coupled portions on the FS Tu . A Similar situation also exists in the cuprates where the in-plane normal state properties are dominated by near-nodal region Valla-SciencePRL . The antinodal regions reduce the in-plane conductivities, while the opening of the pseudogap improves them.

Another similarity with the cuprates exists in that the pseudogap ($`T^{}`$) and $`T_C`$ lines in the $`Tx`$ phase diagram of cuprates show similar behavior to the $`T_{CDW}`$ and $`T_C`$ lines in the $`Tpressure`$ phase diagram of dichalcogenides. As $`T^{}(T_{CDW})`$ goes down, $`T_C`$ increases with doping (pressure). This analogy suggests that the cuprate phase diagram might be shaped by similar competing orders. However, the analogy is no longer valid in the overdoped regime where $`T_C`$ turns back down, even when the pseudogap no longer exists, suggesting that some other, more exotic coupling mechanism might be acting in the cuprates. Here, we would like to point out that the high energy of the ”kinks” observed in TaSe2 Valla-TaSe  and recently in NbSe3 NbSe3  rules out the phonon scenario and represents an unambiguous evidence that the ”kinks” may indeed be caused by other mechanisms. Consequently, the ”kinks” observed in the cuprates PDJ-Lanzara  do not necessarily reflect the electron-phonon coupling .

In conclusion, we have detected a strong anisotropy of the self-energy in a layered dichalcogenide, 2H-NbSe2, with the electron-phonon coupling constant $`\lambda `$ ranging from 0.8 to 1.9 on Nb-derived sheets. The strongest coupling has been found on the inner K-H cylinder, which plays a central role in both CDW and SC transitions in 2H-dicalcogenides. The anisotropy in coupling strength does not induce the anisotropy in the SC gap.

###### Acknowledgements.

The authors would like to acknowledge useful discussions with Tim Kidd, Saša Dordević, Laszlo Forró and Antonio Castro Neto. The BNL program was supported by the US DOE under contract number DE-AC02-98CH10886. The Boston University program was supported by the US DOE under DE-FG02-98ER45680. The sample preparation in Lausanne was supported by the NCCR research pool MaNEP of the Swiss NSF.


format: text

text:
# Database of optical properties of cosmic dust analogues (DOP)

## 1 Introduction

Dust has been detected in most astronomical objects — from the Solar system to active galactic nuclei, and everywhere the dust grains were found to play an important role in different physical and chemical processes. Modelling of these objects should include accurate consideration of scattering, absorption and emission of radiation by the dust grains.

So far in such a modelling one mainly utilized the model of homogeneous spherical particles, though there are many evidences that the interstellar grains are mostly non-spherical and inhomogeneous. Today enough possibilities exist to avoid the use of the spherical model. To help scientists to make this step, we have developed a database including useful information, data, codes, references and links and made it accessible via the Internet. As far as we know there are no analogous databases in the World Wide Web.

In this paper our database having the address http://www.astro.spbu.ru/DOP is described first briefly, and then partly in more detail. The database contains various introductory notes written for students and beginners and more specific information on light scattering theory and its applications. The DOP includes many pieces of the review “Optics of cosmic dust” by Voshchinnikov (Vo02 (2002)) and materials being used in lectures on physics of interstellar matter in the St. Petersburg University.

## 2 General overview of the database

The information, resources, and links presented in the DOP are collected in three main blocks:

1) those related with the theoretical aspects — scatterer models, methods of their realization, bibliography;

2) those connected with calculations (and hence applications) — optical constants, light scattering codes, graphic and tabular libraries of the optical properties, special tools;

3) those coupled with other topics (e.g., radiative transfer modelling).

The database has the following sections and subsections:

* Theoretical aspects

+ Definitions of the optical characteristics of non-spherical scatterers
    + Scatterer models and their parameters
    + Exact and approximate methods of light scattering theory
    + Light scattering experiments
    + Bibliography of light scattering works
* Application aspects

+ Optical constants of materials
    + Light scattering codes
    + Benchmark results
    + A graphic library of the optical properties
    + Self-training algorithm for calculation of the optical properties of fractal-like aggregates
* Related resources

+ Radiative transfer tools
    + Miscellaneous

Let us now turn to a brief description of the subsections.

### 2.1 DOP subsections on theoretical aspects (part A)

Here we first define the cross-sections, albedo, scattering matrix and other basic optical characteristics of non-spherical scatterers and compare them with those for spherical particles. The idea is to supplement and compress information given in other sources of such definitions (see Mishchenko et al., Mi00 (2000), mishetal02 (2002))

The parameters of non-spherical scatterers (size, shape, etc.) are briefly discussed in the next subsection and the light scattering models popular in astronomy (spheres, cylinders, spheroids, aggregates) are listed in there.

To realize any scattering model, i.e. to obtain its optical properties, one should select a proper exact or approximate light scattering method. We make an overview of various light scattering methods and cite several other useful reviews of this kind. A special very detailed review of the methods was prepared by N.G. Khlebtsov (see Sect. 3.1 for more details)

The optical properties of cosmic dust analogues can be also measured in laboratory. References to two available WWW databases of such measurements are presented.

The last subsection is devoted to bibliography on light scattering theory and its various applications which is considered in detail in the Sect. 3.2.

### 2.2 DOP subsections on application aspects (part B)

Here we start with definitions of optical constants (refractive indices or dielectric functions) of materials and link the Jena-St. Petersburg Database of Optical Constants for astronomy (JPDOC) developed earlier. This database described by Henning et al. (1999) and Jäger et al. (2003) consists of original laboratory data, references and links to resources related to measurements and calculations of the optical constants for materials important for astronomical applications. More information on this part is presented in Sect. 3.3.

The next subsection contains a set of light scattering codes developed by us and used in astronomical applications. Links to other collections of such codes are also given.

To test a code not used before one often needs some well checked results. The optical characteristics of spheroidal particles of different size, shape, and refractive index calculated by several methods are given in our benchmark subsection.

An important part of the database is the graphic library of the optical properties of different scatterer models — homogeneous and core-mantle spheres, homogeneous and core-mantle cylinders, homogeneous spheroids. In another part of this subsection the properties of dust grains which can be estimated from the observations of the interstellar extinction and polarization are discussed (see Sect. 3.4 for more details).

The DOP also includes description and computer realization of a self-training algorithm to predict the optical properties of fractal-like aggregates. As this tool was not yet properly outlined in the literature, Sect. 3.5 is devoted to its brief description.

### 2.3 DOP subsections with related resources (part C)

The optical properties are very often used in radiative transfer (RT) calculations. Therefore, we concern this point as well as list the RT methods and the majority of the RT codes used in astronomy.

The last subsection contains links to different related resources (like other databases of light scattering codes, electronic newsletters, etc.) which might be of some use in the DOP context.

## 3 A more detailed description of some parts of the database

### 3.1 Review of light scattering methods for non-spherical particles (A3)

The review consists of two parts — an bibliographic overview of history and development of several basic methods and a comparison of all available approaches to solution of the light scattering problem for non-spherical particles.

The overview made by N.G. Khlebtsov concerns both the exact (separation of variables, point matching, integral equations, coupled dipoles, T-matrix) and approximate \[Rayleigh, Rayleigh–Gans–Debye (RGD) and its generalizations, anomalous diffraction, eikonal, geometrical optics (GO), perturbation, etc.\] methods. In this review the history of development of the methods is traced in very careful manner. The works on the applicability range and various applications of the methods are referred as well, and all together there are over 500 cited papers.

The comparison of methods is made in form of a table for 11 main approaches:
— separation of variables method,
— discretized Mie formalism,
— finite element method,
— finite difference method,
— point matching method,
— generalized multipole technique,
— extended boundary conditions method,
— method of moments,
— coupled dipoles method,
— Fredholm integral equation method,
— ray tracing/Monte Carlo method.

These approaches are classified and their ability to treat particles of different shape and structure as well as anisotropic and chiral scatterers and ensembles of particles is considered. For each method, references to the ansatz and some recent reviews are given. Comments to the table show the relationships and differences of the approaches, briefly concern their classifications and features, sources of available codes, and refer to works on comparison of the methods and other reviews.

We have also started the numerical investigation of applicability ranges of most popular methods using the spheroidal model of scatterers. An example of the results is given in Fig. 1 (see Il’in et al., Il02 (2002) for more details).

### 3.2 Bibliography on light scattering and its applications (A5)

Besides the review, the DOP also includes a special searchable bibliographic database on light scattering theory and its various applications. This bibliography being collected for over 15 years by V.A. Babenko consists of the following parts:
spherical homogeneous particles — methodological questions, functions, calculations, internal fields, resonances, near field, etc.
inhomogeneous spheres — layered, with radial and chaotic internal inhomogeneities, with inclusions, resonances, etc.
approximations for homogeneous spheres — GO, Rayleigh, RGD, eikonal, van de Hulst (anomalous diffraction), and so on approximations;
complicated problems — beam scattering, scattering by charged spheres, by rough, anisotropic and chiral particles, etc.;
cylindrical particles (see below);
non-spherical objects — works on different approaches and methods;
specific scattering media — interstellar particles, marine particles, metallic zoles, powders, etc.;
carbon;
ice;
water;
agglomerates, fractals;
non-linear and mechanical effects in aerosols;
vaporization and explosion of particles;
inverse problems;
laboratory experiments;
laser remote sensing;
radiative transfer;
— etc.

To demonstrate the wideness of the questions considered we shall exemplify a full list of subsections for the section on cylindrical particles:
— circular infinite cylinders;
— oblique incidence and orientation averaging;
— multi-layered cylinders;
— GO approximation (including hexagonal cylinders);
— cylinders with smooth radial inhomogeneity;
— gyrotropic, gyroelectric and anisotropic cylinders;
— internal and near fields;
— experiments on cylinders;
— morphological resonances;
— van de Hulst approximation for cylinders;
— Rayleigh and RGD approximations;
— Hart-Montroll and S-approximations;
— non-circular infinite cylinders;
— circular cylinder of finite length and discs.

The total number of included references exceeds 8000 and each subsection contains several dozens of references.

### 3.3 Optical constants for cosmic dust analogues (B1)

This part of the DOP includes some notes introducing the optical constants required by astronomy, and links the optical constants database (JPDOC). The JPDOC consists of two main parts. One is the bibliography on optical constants of astronomically interesting materials:
— amorphous/glassy/crystalline silicates of different kinds;
— silicon, SiO, crystalline/fused SiO2;
— metals: Fe, Mg, and others;
— oxides: FeO, Fe2O3, Fe3O4, MgO, Al2O3, MgAl2O4;
— sulfides: FeS, MgS, SiS2;
— carbides: SiC, FeC, TiC;
— carbonaceous species: diamonds, graphite, coals, kerogens, HAC, glassy/amorphous carbon, PAHs and so on;
— organics: tholin, “organic refractory”, etc.
— ices: H2O, CO, CO2, NH3, HCN, etc. and their mixtures;
— FeSi, CaCO3 and some other materials.
We tried to include all papers published since the beginning of the previous century (all together over 700 references).

Another part is a set of data specially measured in the Jena laboratory for astronomical purposes. Note that although various terrestrial analogues of cosmic solids have been studied in chemical and physical laboratories, many of these experiments neither took into account the specifics of cosmic dust materials (composition, lattice structure, etc.) and conditions (low temperature, processing, etc.), nor covered the wavelength intervals of the current astrophysical interest.

A summary of the data obtained in Jena and made available is given in Table 1; data examples are presented in Fig. 2 (see Jäger et al., 2003 for more examples).

### 3.4 Graphic library of optical properties of non-spherical scatterers (B4)

This DOP section consists of two parts — a basic graphic library and special notes on optical properties of cosmic dust analogues and their connection with observed manifestations of interstellar dust.

In the basic library we show how different optical characteristics (extinction, scattering and absorption efficiencies, albedo, asymmetry factor, etc.) vary with scatterer parameters: the real ($`n`$) and imaginary ($`k`$) parts of the refractive index, size (the refractive index is kept independent of radiation wavelength), shape, structure. So far we have considered homogeneous and core-mantle spheres and infinite cylinders, and homogeneous spheroids, but are planning to involve more scatterer models. The figures in the DOP are given for weakly absorbing ($`n=1.31.7`$, $`k=00.3`$) and highly absorbing ($`n=1.53.5`$, $`k=0.52.5`$) particles (see Fig. 3 for examples). These refractive indices are typical of astronomically interesting materials in visual.

The special notes are the chapters of the review of Voshchinnikov (Vo02 (2002)) which cover the following questions of observed phenomena modelling — extinction efficiencies: general behaviour and deviations, wavelength dependence, the 2175 A feature, absolute extinction and abundances, etc.; polarization efficiencies: size/shape/orientation effects, linear polarization: wavelength dependence, circular polarization: change of sign and so on. Some more chapters will be included in the near future.

### 3.5 Self-training algorithm for fractal-like aggregates (B5)

This part of the DOP was developed in Kharkov University by V.P. Tishkovets, S.A. Beletsky and P.V. Litvinov. Its idea is as follows — as calculations of the optical properties of subparticle aggregates need a huge amount of processor time, it is worthwhile to construct a self-training algorithm that could learn to predict the properties analyzing a limited set of these data.

To reduce the number of free parameters, fractal-like aggregates of equal size subparticles were selected. The structure of the aggregates can be characterized by two parameters: the fractal dimension $`D`$ and the prefactor constant $`\rho `$ (see, e.g., Feder, 1988). Other parameters are the number of subparticles $`N`$, the size parameter $`x`$ and the real and imaginary parts of the refractive index of subparticles.

The algorithm was based on an artificial neural network named ’multi-layered perceptron’ (see, e.g., Miller, 1990 for more details). A model of perceptron with 6 input and 30 output neurons and two layers with 30 neuron in each between the input and output ones was used. The algorithm was realized in the form of Pascal and Fortran codes. The region of its training was $`n=1.41.7`$; $`k=0.0010.1`$; $`x=1.5`$; $`D=3`$; $`\rho =8`$; $`N50`$. Accurate calculations were done by the Mackowski and Mishchenko (1996) code. The results were the expansion coefficients of the scattering matrix elements in series of the generalized spherical functions for particles of different $`N`$, refractive index, etc. Figure 4 shows the accuracy of the perceptron results.

The codes and data archive as well as documentation are available via the DOP. Our work with the perceptron on modelling of the optical properties of interplanetary dust grains has demonstrated its advantages due to extremely small size of the codes and very fast data access.

## 4 Concluding remarks

We have introduced a database of the optical properties of cosmic dust analogues. It is designed for student and scientists to help in applications of non-spherical scatterer models in astronomy but can be useful in other scientific fields as well.

A further development of the database is planed in near future. It will include a graphics interface for the JPDOC, a special tool for on-line calculations of different optical properties of non-spherical layered particles, more chapters from Voshchinnikov’s (Vo02 (2002)) review and its continuation, and so on.

###### Acknowledgements.

The authors thank V.G. Farafonov, S.I. Grachev and D.I. Nagirner for useful discussions and comments. The work was partly supported by the INTAS (Open Call 1999 grant N 652), the Russian Ministry of High Education, the Russian federal program Astronomy and that for prominent scientific schools.


text:
# A corner transfer matrix renormalisation group investigation of the vertex-interacting self-avoiding walk model

## 1 Introduction

Self-avoiding walk models are of interest as models of polymers in solution as well as realisations of $`O(n)`$ spin models in the limit $`n0`$.

The self-avoiding walk problem is a difficult problem numerically because of the intrinsic long ranged nature of the problem; self-avoiding walks are extended objects. The numerical methods most used to date for the study of such models are Monte Carlo simulations, series expansions and transfer matrix calculations.

Monte Carlo simulations are hampered by problems of trapping in the low temperature, or equivalently large interaction, phases. These may be relieved using the fluctuating bond method, which also enables the investigation of the dense walk phases. This method has the drawback of modifying the model in a way which is not always convenient, particularly if frustration effects are of interest. Current Monte Carlo methods, such as the PERM method and parallel tempering, relieve these difficulties to some extent, but at the cost of computational complexity.

Series expansions are limited by the lengths of walks that may be enumerated. In practice these are limited to a few tens of lattice steps. Very technical methods may be used to achieve longer series expansions of the partition function (of the order of 60 lattice steps), and in some cases the estimate of the critical point is very impressive. If correlation length exponents are required, then the range of walk lengths open to enumeration are of the order of 30 steps. This becomes increasingly limiting as the complexity of the model increases.

Transfer matrices have the advantage over Monte Carlo simulations of giving numerically exact values of the free energy, correlation lengths, polymer densities,etc., for strips of infinite length, but of finite width. Finite size analysis methods then allow the extrapolation of the finite size estimates. The widths of lattice accessible are limited: for the self-avoiding walk model widths up to 12 or 14 are possible (see for example ), but this again drops as the complexity of the model increases (see for example ). It often occurs that the number of data points one may use in finite size analysis becomes dangerously small for finite size scaling techniques to be of use, particularly if odd/even parity effects are present.

In this article we apply a recent extension of the corner transfer matrix renormalisation group (CTMRG) method to a class of self-avoiding walks introduced by Blöte and Nienhuis, which we shall refer to as the vertex-interacting self-avoiding walk model, to reflect the fact that the monomer-monomer interactions are added for doubly visited sites (for a full description of the model see Section 2.1). We will show that the CTMRG method extends the possibilities of transfer matrices to larger system sizes, at the cost of introducing a phase space truncation scheme. The accuracy of the results is however controllable through a parameter $`m`$ corresponding to the number of states kept from one iteration to the next.

The other important result in this paper relates to the vertex-interacting self-avoiding walk model. Whilst it has the standard collapse transition as a special case (where the walk is obliged to turn after each step), it is believed to have a generic collapse transition in a different universality class. We show that this collapse transition has a first order character, despite the fact that there exists critical exponent $`\nu =12/23`$ which has been exactly determined by Warnaar et al.

## 2 The corner transfer matrix renormalisation method

Transfer matrices, as a numerical method, have a major advantage over Monte Carlo methods in that they sum over all possible configurations, and so are not hampered by trapping effects at low temperatures (or high densities in polymer problems) which plague Monte Carlo simulations of models with competing interactions (frustrated models).

The major inconvenience is the large amount of computer memory required to store the matrices, limiting the maximum lattice width in practical calculations. This is alleviated in part by finite-size analysis methods which work well, at least for the simpler models, coupled with efficient extrapolation methods.

In two-dimensional classical spin models, or one-dimensional quantum spin models, there exist iterative approximation schemes which, by successive pruning of phase space, permit calculations for much larger lattice widths for given computer resources. These methods are the Density Matrix Renormalisation Group method (DMRG) and the Corner Transfer Matrix Renormalisation Group method (CTMRG). These methods have little recourse to the renormalisation group, and are better thought of as iterative matrix approximation schemes. The basic idea is that the (approximate) transfer matrix for a lattice of size $`2N`$ is calculated from the transfer matrix for a lattice of size $`2N2`$ by inserting extra lattice sites. A change of basis is then performed, such that the transfer matrix may be projected onto a smaller basis, with the smallest possible loss of information.

In a recent paper we showed how the CTMRG method, first introduced for classical spin systems by Nishino, may be applied to self-avoiding walk type models. This requires a non-trivial modification of the method since the iteration scheme requires all interactions to be local. The self-avoiding walk problem, by its very nature, is non-local. In the remainder of this section the details of the method, absent in reference, are provided. In particular we show how this non-locality may be overcome, and describe the implementation of CTMRG for this class of problem.

### 2.1 The self-avoiding walk, the $`O(n)`$ model and vertex models with complex weights

It is known that spin models with $`O(n)`$ invariant spins have high temperature expansions consisting of closed loops. The self-avoiding loop (SAL) model is related to the high temperature expansion of the $`O(n=0)`$ model. In a similar way, the high temperature expansion of the susceptibility leads to a model with one open path, and a gas of loops. The $`n0`$ limit eliminates the loops to leave a self-avoiding walk. The partition function for a self-avoiding walk is thus equivalent to the susceptibility of an $`O(n=0)`$ spin model.

The connection between the $`O(n)`$ class of models and better known models is as follows: $`n=1`$ corresponds to the Ising model, $`n=2`$ to the XY model and $`n=3`$ to the classical Heisenberg model. The partition functions of these models are given by:

$$𝒵=\underset{\{\stackrel{}{s}_i\}}{}\mathrm{exp}\left(\frac{1}{2}\beta J\underset{i,j}{}\stackrel{}{s}_i\stackrel{}{s}_j\right),$$
(1)
where $`i,j`$ refers to a sum over nearest neighbour spins. Whilst the high-temperature expansion of the Ising model is relatively easy, the expansion for other values of $`n`$ is a little cumbersome. In what follows we shall consider a slightly different $`O(n)`$ invariant model, due to Nienhuis, with the partition function given by

$$𝒵_{O(n)}=\underset{\{s_i^\alpha \}}{}\underset{i,j}{}\left(1+K\underset{\alpha =1}{\overset{n}{}}s_i^\alpha s_j^\alpha \right),$$
(2)
where the spins, $`s_i^\alpha =\pm 1`$, are placed on the lattice bonds. A diagrammatic expansion of Equation 2 follows if we identify the $`1`$ as the weight of an empty bond between the sites $`i`$ and $`j`$ and the $`K`$ as the weight of an occupied bond. Since $`s_i^\alpha =0`$ and $`\left(s_i^\alpha \right)^2=1`$, the only terms in the expansion of Equation 2 which survive the sum over $`\{s_i^\alpha \}`$ are those terms consisting of closed loops where the ‘colour’ $`\alpha `$ is conserved around each loop. Different loops may have different ‘colours’. The sum over $`\alpha `$ gives a factor $`n`$ for each loop. An equivalent expression for Equation 2 may now be given in terms of loop graphs $`𝒢`$:

$$𝒵_{O(n)}=\underset{𝒢}{}n^{l(𝒢)}K^{b(𝒢)},$$
(3)
where $`l`$ is the number of loops and $`b`$ is the number of occupied bonds. The partition for self-avoiding loops is then just

$$𝒵_{SAL}=\underset{n0}{lim}\frac{1}{n}𝒵_{O(n)}.$$
(4)
The parameter $`n`$ is now a fugacity controlling the number of loops, and in this representation need no longer be taken as an integer. This fugacity corresponds to a long-range interaction, since the loops may be of any size. This non-locality is undesirable for our purposes. We would like to express $`n`$ as a product over local weights despite the variable size of the loops. The only real numbers with this property are $`0`$ and 1. This leads us to introduce complex local weights, $`w_i`$. If $`|w_i|=1`$, the product of the $`w_i`$ around a loop will also be of modulus 1. The final weight must, however, be real. This may be achieved if for every complex weight, we also have its complex conjugate. Each loop may be followed clockwise or anticlockwise. The loops in Equation 3 are not oriented, but may be oriented by associating $`2^{l(𝒢)}`$ oriented graphs with each non-oriented graph. The loop fugacity is not required to be the same for the two orientations of the loop, and so a fugacity $`n_+`$ ($`n_{}`$) is associated with (anti)clockwise oriented loops. The partition function may then be written

$`𝒵_{O(n)}`$ $`=`$ $`{\displaystyle \underset{𝒢}{}}(n_++n_{})^lK^b,`$ (5)
$`=`$ $`{\displaystyle \underset{𝒢^{}}{}}n_+^{l_+}n_{}^l_{}K^b`$ (6)
where $`𝒢^{}`$ is the set of oriented loop graphs and $`l_+`$ ($`l_{}`$) is the number of (anti)clockwise oriented loops. On the square lattice there must be four more corners with one orientation, compared to the other orientation, in order to close a loop. Using this fact we may set $`w_i=\mathrm{exp}(i\theta /4)`$ for a clockwise corner, $`w_i=\mathrm{exp}(i\theta /4)`$ for an anticlockwise corner and $`w_i=1`$ otherwise. This leads to $`n_+=\mathrm{exp}(i\theta )`$ and $`n_{}=\mathrm{exp}(i\theta )`$, and hence $`n=2\mathrm{cos}(\theta )`$. The oriented loop factor has now been broken up into local weights and the partition function may be rewritten in terms of a vertex model:

$$𝒵=\underset{𝒢^{}}{}\underset{i}{}v_i$$
(7)
where $`v_i`$ is the weight of the vertex found at site $`i`$. The possible vertices are shown in Figure 1. The derivation given here is only for the simplest case, but we may freely change the weights of the vertex configurations in order to generate different interactions in the original model. Notably, to have exactly the standard self-avoiding walk model it is necessary to set the weights of vertices 14 to 19 to zero, and set $`n=0`$, or equivalently set $`\theta =\pi /2`$.

A vertex model may be thought of as a spin model where the spins live on the bonds of the lattice. Each spin has three states corresponding to the empty lattice bond, a bond oriented in the positive $`x`$ ($`y`$) direction and a bond oriented in the negative $`x`$ ($`y`$) direction. The vertex weights may be written as a vertex function $`W(\sigma _i,\sigma _j,\sigma _k,\sigma _l)`$, where the tensor elements correspond to the vertex weights of Figure 1. The vertex function was defined using the following convention:

$$\sigma =\{\begin{array}{cc}1\hfill & \text{if the bond has an arrow to the right or upwards,}\hfill \\ 0\hfill & \text{if the bond has no arrow,}\hfill \\ 1\hfill & \text{if the bond has an arrow to the left or down.}\hfill \end{array}$$
(8)
### 2.2 The vertex-interacting self-avoiding walk model

In this section we apply the CTMRG method explicitly to the vertex-interacting self-avoiding walk model introduced by Blöte and Nienhuis. This model consists of a random walk which is allowed to collide at a site, but is not allowed to cross or occupy a bond more than once. An interaction energy $`\epsilon <0`$ is introduced for each collision, which replaces the attractive nearest-neighbour interaction present in the standard $`\mathrm{\Theta }`$-point model. A step fugacity, $`K`$, is introduced in order to control the average length of the walk. A stiffness is also introduced by weighting sites sitting on a bend in the walk differently from sites sitting on straight portions of the walk. Following the convention of Blöte and Nienhuis, we choose to add a weight $`p`$ to the straight segments of the walk. The choice $`p=0`$ then corresponds to a walk in which straight segments are eliminated (the walk is obliged to turn through a right angle after each step). The partition function for the model is then given by

$$𝒵=\underset{\text{walks}}{}(Kp)^{N_s}p^{N_c}\tau ^{N_I},$$
(9)
where $`\tau =\mathrm{exp}(\beta \epsilon )`$, $`N_s`$ is the length of the walk (number of steps), $`N_I`$ is the number of doubly occupied sites and $`N_c`$ is the number of corners.

### 2.3 The corner transfer matrix renormalisation group method

Following Baxter, the partition function of a two-dimensional lattice model may be written in terms of the product of four matrices representing the four quarters of the lattice, see Figure 2. The inputs and outputs of the matrices are the configurations at the seams of the four quarters. These matrices are known as corner transfer matrices. In general the four matrices are different, but may often be related by lattice symmetries. For our model the four matrices are the same up to a complex conjugation operation.

In general it is difficult to explicitly calculate these matrices for systems with a large number of sites. This is where the CTMRG method comes in; the matrices for larger lattices are calculated iteratively from smaller lattices. This is done as follows. A prototype system consisting of a small number of spins, four in the case of a vertex model, is set up exactly, see Figure 3. At each iteration the system is enlarged by adding spins, and this enlarged system is then projected back onto the prototype system in some optimal way, so as to minimise the loss of information. The spins of the prototype system have $`m`$ states, where $`m`$ determines the amount of information which may be carried forwards at each iteration.

In addition to the vertex function $`W(\sigma _i,\sigma _j,\sigma _k,\sigma _l)`$, we also introduce the vertex functions $`P_N(\sigma _i,\xi _j,\xi _k)`$ and $`C_N(\xi _i,\xi _j)`$ where the $`\xi _i`$ are the $`m`$-state spins of the prototype system. The matrix $`C_N`$ corresponds to the estimate of the corner transfer matrix after $`N`$ iterations. The partition function for a $`2N\times 2N`$ lattice is given by

$$Z_N=\mathrm{Tr}\left(C_NC_N^{}\right)^2,$$
(10)
where the $``$ denotes complex conjugation. The lattice is enlarged by introducing additional vertices, as shown in Figure 4. A corner matrix for the enlarged lattice is given by

$$C_{N+1}^{}=WP_NP_N^{}C_N,$$
(11)
where $``$ indicates that the appropriate contractions over indices are performed. In a similar way we introduce an enlarged vertex function

$$P_{N+1}^{}=WP_N.$$
(12)
By considering carefully the symmetry of the vertex function $`W`$ it may be shown that the partition function for the enlarged lattice $`2(N+1)\times 2(N+1)`$ is given by

$$Z_{N+1}^{}=\mathrm{Tr}(C_{N+1}^{}C_{N+1}^{})^2.$$
(13)
We now wish to project the enlarged model, defined through $`C_{N+1}^{}`$ and $`P_{N+1}^{}`$, back onto the smaller prototype system, giving $`P_{N+1}`$ and $`C_{N+1}`$.

Since the configuration space is smaller after projection, some information is lost. The trick is to minimise this information loss. This is done by performing a change of basis in some optimal manner. The partition function $`Z_{N+1}^{}`$ is invariant by any change of basis, but the partial trace, giving $`Z_{N+1}`$ after the change of basis, is not. For a model with real weights it is clear that the optimal choice will minimise $`\epsilon =|Z_{N+1}^{}Z_{N+1}|`$. This is done by diagonalising the matrix $`\stackrel{~}{z}=(C_{N+1}^{}C_{N+1}^{})^2`$ and keeping the basis vectors corresponding to the largest eigenvalues. In a model with real weights $`\stackrel{~}{z}`$ is just the density matrix, up to a constant. In the case considered here the eigenvalues of $`\stackrel{~}{z}`$ are complex, coming in conjugate pairs. Minimising $`\epsilon `$ corresponds to taking the basis vectors for the eigenvalues with the largest real parts, since the imaginary parts cancel. This may well involve dropping terms contributing to $`Z_{N+1}^{}`$ with larger moduli than other terms which are kept. It is not clear how this effect propagates with the number of iterations. Another choice is to keep the terms contributing to $`Z_{N+1}^{}`$ with the largest absolute values. This is found to give better results in practice. To complete the iteration scheme, we must also calculate $`P_{N+1}`$. This is done by projecting $`P_{N+1}^{}`$ onto $`P_{N+1}`$ using the same basis vectors. As this iterative process is repeated, the partition function and other thermodynamic quantities are calculated for larger and larger lattices. The value of $`m`$ defines the size of the configurational space kept from one iteration to the next. The larger the value of $`m`$, the better the approximation.

The iteration scheme not only gives an approximation of the partition function for a square system, but also gives an approximation for the transfer matrix from which an infinite strip may be constructed. The transfer matrix is given directly in terms of the vertex functions $`P_N`$ and $`W`$ as

$$𝒯=P_NWWP_N^{}$$
(14)
for even lattice widths (width $`L=2N+2`$) and

$$𝒯=P_NWP_N^{}$$
(15)
for odd lattice widths (width $`L=2N+1`$). Free energies and correlation lengths for strips may then also be calculated, and the arsenal of finite-size scaling methods available for strips may be used. CTMRG is most efficient, however, when it is used to calculate directly one-point functions, such as the density, and this is how we use it in this article.

The density is most simply calculated by selecting one of the central lattice bonds. The density $`\rho `$ is then given by $`\rho =|\sigma |`$, where $`\sigma `$ is the value of the spin on the selected bond. In terms of the corner transfer matrices,

$$\rho _L=\frac{\mathrm{Tr}(\mathrm{\Lambda }C_N^{}C_N^{}C_N^{}C_N^{})}{\mathrm{Tr}(C_N^{}C_N^{}C_N^{}C_N^{})}$$
(16)
where $`\mathrm{\Lambda }`$ is a matrix whose elements take the value $`|\sigma |`$ and $`L=2N`$ is the linear dimension of the lattice.

In spin models the one-point functions are often calculated by the infinite lattice method, where the CTMRG method is iterated for a fixed value of $`m`$ until the quantity of interest (magnetisation or energy per spin) no longer changes. This occurs because the number of eigenstates kept at each iteration is no longer sufficient to propagate the two point correlations throughout the system. The limiting value of the one point correlation function may then be studied as a function of $`m`$. In such a manner the magnetisation or energy per site may be estimated for the infinite system directly. This method relies on the rather nice way in which CTMRG breaks down for spin models. Unfortunately the complex nature of the vertex models at the heart of our calculations gives rise to more complicated behaviour as $`m`$ is increased. Figure 5 shows the density of the self-avoiding walk model ($`p=1`$, $`\tau =0`$ and $`n=0`$, or equivalently $`\theta =\pi /2`$) as a function of $`L`$ for different values of $`m`$. It is clearly seen that there is no limiting value for the density as $`L`$ is increased. Once the value of $`m`$ is insufficient, $`\rho `$ has a chaotic oscillatory behaviour with $`L`$. Instead of fixing $`m`$ and varying $`L`$, we choose to do the opposite. Taking the density as an example, we define $`\rho _L(m)`$ as the density for a system of linear dimension $`L`$ for a given value of $`m`$. The limiting value of $`\rho _L(m)`$ is calculated by increasing the value of $`m`$ for fixed $`L`$ until the required accuracy for $`\rho _L`$ is reached. This is then studied as a function of $`L`$ using finite-size scaling methods, described below.

## 3 Results

### 3.1 The vertex-interacting model with $`p=0`$

When $`p=0`$ the model is expected to have a critical behaviour equivalent to the standard $`\mathrm{\Theta }`$-point model; for small enough $`\tau `$ the average length of the walk diverges as a critical value of the step fugacity, $`K_c`$, is reached. At $`K_c`$ the density is still zero and the gyration radius, $`R_G`$, (or any other characteristic measure of the size of the walk) diverges as $`KK_c^{}`$ following a power law defining the correlation length exponent $`\nu `$:

$$\xi (K_cK)^\nu .$$
(17)
For $`\tau `$ smaller than some tricritical value, $`\tau _{tc}`$, the walk is in the same universality class as the usual self-avoiding walk, and $`\nu =3/4`$. This behaviour terminates at the tricritical point $`(K_{tc},\tau _{tc})`$, where $`\nu =4/7`$ as for the standard $`\mathrm{\Theta }`$-point. For values of $`\tau >\tau _{tc}`$ the length of the walk diverges discontinuously at a first order transition. The phase diagram of the $`p=0`$ vertex-interacting model is qualitatively the same as for the $`\mathrm{\Theta }`$ point, and is shown schematically in Figure 6. Unlike for the standard $`\mathrm{\Theta }`$-point model, the location of the tricritical point is known exactly to be $`(K_{tc}=1/2,\tau _{tc}=2)`$.

The idea is now to study this model using the CTMRG method outlined. As mentioned in section 2.3, CTMRG gives as a byproduct an approximation of the transfer matrix for a strip of width $`L`$ through (14) or (15), depending on the parity of L, from which the phase diagram and critical exponent $`\nu `$ could be calculated using the phenomenological renormalisation group due to Nightingale. This turns out, at least for the case at hand, not to be the most efficient way of approaching the problem, due to the additional computational effort required to calculate the relevant eigenvalues of $`𝒯`$. The CTMRG method is far more efficient for the calculation of one-point functions, such as the bond density, which may be calculated in the middle of the system. With this in mind we use a phenomenological renormalisation scheme based directly on the scaling behaviour of density.

The scaling ansatz states that the singular part of the free energy is a homogeneous function of its arguments. This leads to the following finite size scaling expression:

$$f_s(K,L)=L^d\stackrel{~}{f}(|KK_c|L^{1/\nu }),$$
(18)
where $`d`$ is the spatial dimension, here $`d=2`$. Taking the first derivative,

$$\rho _s(K,L)=L^{1/\nu 2}\stackrel{~}{\rho }(|KK_c|L^{1/\nu }),$$
(19)
In general

$$\rho (K,L)=\rho _{\mathrm{}}(K)+\rho _s(K,L),$$
(20)
and so it is necessary to know the value of $`\rho _{\mathrm{}}(K)`$ in order to exploit fully the scaling behaviour. The density of a finite walk on an infinite lattice is zero, and so $`\rho _{\mathrm{}}(K<K_c)=0`$. This enables the setting up of a phenomenological renormalisation group method using the function

$$\phi _{L,L^{}}(K)=\frac{\mathrm{log}\left(\rho _s(K,L)/\rho _s(K,L^{})\right)}{\mathrm{log}\left(L/L^{}\right)}.$$
(21)
Using (19) we find that $`\phi _{L,L^{}}=1/\nu 2`$ independently of $`L`$ and $`L^{}`$. Naturally there are additional finite size corrections which should be taken into account, but the conclusion is that if the function $`\phi _{L,L^{}}(K)`$ is plotted for various values of $`L`$ and $`L^{}`$ then it will converge to a fixed point given by $`\phi (K_c)=1/\nu 2`$. In what follows we have set $`L^{}=L2`$ and looked for solutions of the equation

$$\phi _{L,L2}(K_c^L)=\phi _{L2,L4}(K_c^L).$$
(22)
If such solutions exist then $`K_c=lim_L\mathrm{}K_c^L`$ and $`\nu =lim_L\mathrm{}1/(2+\phi _{L,L2}(K_c^L))`$.

Preliminary results using CTMRG were recently reported in reference ; notably the density profiles and values for $`K_c`$ and $`\nu `$ are given for various values of $`\tau `$ with $`p=0`$. When $`p=0`$, the vertex-interacting $`\mathrm{\Theta }`$ point is known to occur exactly at $`K_{tc}=1/2`$ and $`\tau _{tc}=2`$. In Figure 7 we present results using this scheme where we have fixed $`p=0`$ and $`\tau =\tau _{tc}=2`$. The results converge rather quickly to the exact values, $`K_c=1/2`$ and $`\nu _\mathrm{\Theta }=4/7`$. Whilst here the tricritical values of $`K`$, $`\tau `$ and $`\nu `$ are known exactly, this is in general not the case. It is interesting to know how well the method works with no external input. As we pass through $`\tau _{tc}`$ the value of $`\nu `$ changes from $`\nu _{\mathrm{SAW}}=3/4`$ for $`\tau <\tau _{tc}`$ to $`\nu _\mathrm{\Theta }=4/7`$ for $`\tau =\tau _{tc}`$ and on to $`\nu =1/2`$ for $`\tau >\tau _{tc}`$. This last value is simply a reflection of the fact that the density is nonzero, and so $`R_GN^{1/2}`$. Clearly the finite size estimates, $`\nu _L(\tau )`$, will be continuous functions of $`\tau `$, tending to the step function only in the infinite lattice limit. If the different $`\nu _L(\tau )`$ cross, then these intersections necessarily tend, in the limit of an infinite system, to the correct values: $`\nu _{tc}=\nu _\theta `$ and $`\tau _{tc}`$. This is a standard method used to determine numerically the location, for example, of the $`\mathrm{\Theta }`$ point in the standard model (see for example ). The estimates derived from this method are shown in Figure 8 plotted as a function of $`1/L`$. The estimates converge convincingly to the exact values $`\tau _{tc}=2`$, $`K_{tc}=1/2`$ and $`\nu _{tc}=4/7`$.

### 3.2 The vertex-interacting model with $`p=1`$

What is interesting about the vertex-interacting model is that the phase diagram is very different when $`p0`$. Naïvely one would expect the $`p=1`$ model to also be equivalent to the $`\mathrm{\Theta }`$-point model; the model consists of a walk which is essentially self-avoiding with short-ranged attractive interactions. This is in fact not the case; an additional transition arises in the dense walk phase and the nature of the transition at the (now) multicritical point changes. There are no exact results for the location of the multicritical point $`(K_{mc},\tau _{mc})`$ as was the case for $`p=0`$. There exist, however, exact predictions for the universality classes of the different phase transitions. The low-density phase transition is expected to be in the self-avoiding walk universality class (with $`\nu =3/4`$). The high-density critical line is expected to be in the Ising universality class ($`\nu =1`$). The value of $`\nu _{mc}`$ at the multicritical point has been found from an exact result for the 19 vertex model with complex weights used in our calculations. It is predicted that $`\nu _{mc}=12/23`$. In the remainder of this section we calculate the phase diagram and compare critical exponent estimates calculated from the CTMRG method with the known results.

In figure 9 the density is plotted as a function of $`K`$ for $`p=1`$ and $`\tau =1`$. The presence of two phase transitions may clearly be seen, the first between the zero-density phase and a dense walk phase, and the second between two dense walk phases.

For the first transition, as $`K`$ is increased, it is possible to use the scaling arguments given in section 3.1 since the density in the thermodynamic limit is zero. For the second transition, however, the density is non-zero, and a priori unknown. To apply the scaling arguments it is necessary to estimate the density in the infinite lattice limit, $`\rho _{\mathrm{}}(K)`$. Alternatively it is possible to elimitate $`\rho _{\mathrm{}}`$ by using two lattice sizes. Along with (19) and (20) this gives a new finite-size scaling function

$$\stackrel{~}{\phi }_L(K)=\frac{\mathrm{log}\left(\frac{\rho (K,L)\rho (K,L2)}{\rho (K,L2)\rho (K,L4)}\right)}{\mathrm{log}\left(\frac{L}{L2}\right)},$$
(23)
which, as before, may be used to find estimates of the critical lines by looking for solutions of the equation

$$\stackrel{~}{\phi }_L(K_c^L)=\stackrel{~}{\phi }_{L2}(K_c^L).$$
(24)
In the limit $`L\mathrm{}`$ the solutions of this equation tend to the correct fixed points, such that $`\stackrel{~}{\phi }(K_c^L)1/\nu 3`$. This again enables the calculation of the critical exponent $`\nu `$.

The high-density transition line is expected to be in the Ising universality class, with an exponent $`\nu =1`$. It is important to note that the exponent is used to represent both the critical exponent related to the correlation length and the inverse of the Haussdorf dimension, $`d_H`$. The Haussdorf dimension is the fractal dimension defined through the scaling of the density with the distance from the centre of mass of the walk. This leads to the scaling law:

$$R_GN^{1/d_H},$$
(25)
where $`R_G`$ is the radius of gyration, or any other characteristic measure of the linear size of the walk. Since in the low density regime the gyration radius is proportional to the correlation length, then as the critical line is approached, and $`L\mathrm{}`$, we may identify $`\nu =1/d_H`$. It is common practice to define $`\nu `$ from (25), whence $`\nu =1/d=1/2`$ on the first-order transition line, where the fractal and space dimension become the same. The $`\nu `$ associated with the high-density transition will necessarily not be defined in this geometrical fashion, but corresponds to the true thermodynamic critical exponent related to the correlation length. The estimate of $`K_c`$ and the associated estimates for $`\nu `$ are given in Figure 10 for $`p=1`$ and $`\tau =1`$. The evolution of $`\nu `$ is not monotonic, emphasising the need to go to large enough lattice sizes before having reliable estimates of the exponent. It is clear that $`\nu `$ converges to the expected value $`\nu =1`$, and we estimate the critical fugacity in this case to be $`K_c=1.326\pm 0.002`$.

Figure 11(a) shows the phase diagram where both transition lines are calculated using the scaling function given in (23). In practice there are usually two solutions to (23) corresponding to the same transition; both converge towards the same value of $`K_c`$ as the size of the lattice is increased. However, one of these solutions evolves more slowly than the other as system size is increased, and this is taken as being the relevant solution for the calculation of $`\nu `$ and $`K_c`$. As $`\tau `$ is increased these two solutions come closer (for fixed lattice size), and join at a certain value of $`\tau `$. Beyond this value of $`\tau `$ there is nolonger a solution corresponding to the transition in question (there may be solutions corresponding to other transitions). This vanishing point moves to larger $`\tau `$ as the lattice size is increased. The phase diagram shown in Figure 11(a) is calculated for system sizes up to $`L=30`$.

Figure 11(b) shows the estimates for the zero-density phase transition lines calculated using (21). As $`\tau `$ approaches its multicritical value, the finite sized system feels the presence of the high-density phase transition, leading to a peak in the estimated transition line. This peak reduces in size rapidly as the system size is increased.

Estimates of the location of the multi-critical point, ($`K_{mc},\tau _{mc}`$), may again be found by looking for crossings in the finite-size estimates of $`\nu (\tau )`$ calculated along the self-avoiding walk transition line using (21) and (22). These estimates are given in figure 12. The location of the multicritical point is estimated as $`K_{mc}=0.3408\pm 0.0002`$, $`\tau _{mc}=4.69\pm 0.01`$. These estimates agree very well with the estimates given by Guo, Nienhuis and Blöte, who give $`z_{mc}=K_{mc}^2\tau _{mc}0.54`$, however the finite-size estimates of $`\nu _L(\tau )`$ cross at a value too low to be consistent with the exact value at the multicritical point ($`\nu _L\mathrm{}0.51`$ compared with $`\nu _{\mathrm{exact}}=12/23=0.5217\mathrm{}`$). This leads to the conclusion that either the scaling method is flawed or $`\rho _{\mathrm{}}0`$ at the multicritical point. In Figure 13 the density is plotted as a function of $`\tau `$ along the critical (self-avoiding walk) line. At the multicritical point the density has a clear finite limit as the size is increased, indicating that the transition is first order as approached in this direction. As further evidence of the first order character of the transition when approached following the self-avoiding walk line we plot in Figure 14 the finite-size estimates of $`K_c`$ and $`\rho `$ for different values of $`\tau `$ near $`\tau _{mc}`$. From Figure 11 it is seen that for $`\tau `$ close to $`\tau _{mc}`$ the system is influenced by the high-density transition. As the system size is increased, the estimates transit to the zone of influence of the zero-density transition, where they tend to their limiting values. This transit is clearly seen in Figure 14 for $`\tau =4.3<\tau _{mc}`$. When $`\tau =4.6`$ it is seen that the transit begins, but the sizes we could consider in a practical time scale were not large enough to see the end of the transition. It is however reasonable to suppose that for $`\tau =4.6`$ the limiting density will be zero and occur as $`L\mathrm{}`$. No such transit is seen for $`\tau =4.7>\tau _{mc}`$, which displays a non-zero limiting density. The precise value of this limiting density is of course not reliable, since the location of the transition point, found assuming $`\rho _{\mathrm{}}=0`$, will not be precise. The determination of the multicritical point is however expected to be correct since $`\rho _{\mathrm{}}=0`$ along the entire self-avoiding walk line, as must be the case since the fractal dimension $`d_H=1/\nu <2`$.

There exist five branches of solutions for the Blöte-Nienuis $`O(n)`$ model on the square lattice, through its mapping onto the 19 vertex model, for which the (multi)critical points are exactly known, as well as many of the critical exponents. The multicritical point here is thought to correspond to the same universality class as the $`n=0`$ solution belonging to the third branch. The location of this point is given exactly as:

$$\begin{array}{ccc}\hfill z=K^2\tau & =& \left\{2\left[12\mathrm{sin}(\theta /2)\right]\left[1+2\mathrm{sin}(\theta /2)\right]^2\right\}^1\hfill \\ \hfill K& =& 4z\mathrm{sin}(\theta /2)\mathrm{cos}(\pi /4\theta /4)\hfill \\ \hfill pK& =& z\left[1+2\mathrm{sin}(\theta /2)\right]\hfill \\ \hfill \theta & =& \pi /4\hfill \end{array}\}.$$
(26)
This gives the location of the exactly known multicritical point as $`K_{mc}^{}=0.446933\mathrm{}`$, $`p_{mc}^{}=0.275899\mathrm{}`$ and $`\tau _{mc}^{}=2.630986\mathrm{}`$. We checked the determination of this point using CTMRG fixing $`p=p_{mc}^{}`$ and $`\tau =\tau _{mc}^{}`$ and looking for the limiting value of $`K_L`$. This is shown in Figure 15 and is clearly consistent with the exact known value. The density was then calculated firstly at the finite-size estimate of the critical point and then fixing simultaneously the values of $`K`$, $`p`$ and $`\tau `$ to the values given in (26) (see Figure 16). In both cases a conservative estimate gives $`\rho _{\mathrm{}}=0.43\pm 0.01`$. These results indicate that the “multicritical” point is fairly exotic in that it corresponds generically to a first-order transition but with a critical behaviour if approached in a specific direction. We have tried to recover the critical exponent $`12/23`$ by using the scaling behaviour of the density at $`K_{mc}^{},\tau _{mc}^{}`$ and $`p_{mc}^{}`$. Whilst the results are not inconsistent with $`\nu =12/23`$, they are not of very good quality and are also not inconsistent with a value $`\nu 0.54`$ as originally found by Blöte and Nienhuis using transfer matrices. Whilst in principle there is no limit to either the precision or lattice sizes obtainable using CTMRG, there is a practical limitation in terms of computing time and resources.

## 4 Conclusion

In this paper we have presented a version of the corner transfer matrix renormalisation group method suitable for use with self-avoiding walk type models and related $`O(n)`$ symmetric models. We have applied it to the vertex-interacting self-avoiding walk model. The method was shown to reproduce the known results, and has given new insight into the critical nature of the model when $`p0`$.

Self-avoiding walk models are used as models for polymers in dilute solution, and it is known that the lattice $`\theta `$-point model provides a good description of what happens for real homopolymers in solution. The standard textbook reason for not taking into account the rigidity of the polymer is that the rigidity simply changes the persistence length, or in terms of the lattice model the number of real monomers per step of the walk. The critical behaviour is unaffected in the infinite walk limit. This should be true as long as the fractal dimension $`d_H`$ is smaller than the spatial dimension (i.e. $`d_H<2`$ here). However, when $`d_H=d`$ it may be expected that the polymer feels the presence of the lattice. This will then give rise to competition between the rigidity and the underlying lattice structure, reminiscent of frustration in spin models (indeed the critical behaviour of the Nienhuis $`O(n)`$ model on the square lattice has been related to the fully-frustrated XY model). If such a claim is true, then the critical behaviour observed will be universal (in the sense of not depending on the lattice) and independent of rigidity if $`d_H<d`$, and may be non-universal (in the same restricted sense) and may depend on the rigidity if $`d_H=d`$. One could therefore argue that the simple observation that the critical behaviour of the collapse transition changes as $`p`$ is varied implies that the density is non-zero at the collapse transition when $`p0`$. In interpreting this transition as due to frustration, it should be noted that it does not exist for the $`O(n)`$ model on the hexagonal lattice.

There are other models which have qualitatively similar phase diagrams. Most well known is the self-avoiding trail model. A trail is defined to be self-avoiding on the bonds, but is allowed to visit sites more than once, and in particular is allowed to intersect itself. The attractive interaction is again associated with doubly visited sites (intersection or collision). It has been conjectured that the collapse transition for trails occurs exactly at $`\tau =3`$ with an exponent $`\nu =1/2`$; it is also conjectured that limiting density is zero. This conjecture resulted from a Monte-Carlo enumeration of extremely long walks of the kinetic trail model, believed to be equivalent to the trail model at the collapse transition . As discussed above, in this case what is being calculated is in fact the fractal dimension $`d_H`$. The relation $`\nu =1/d_H`$ is only valid when the critical length may be measured using the gyration radius. As argued above, the validity of this identification may be questioned when $`d_H=d`$. It would be of interest, then, to check by other means the possibility that the criticality is in fact not of the same type as observed here for $`p0`$. This question is currently under investigation.

Another model with a qualitatively similar phase diagram is the Hydrogen-bonding model. This model includes interactions between nearest-neighbour sites, which may be taken into account in the CTMRG method by mapping the walk model onto a 4-state vertex model with complex weights. There are, however, important differences, notably the high-$`\tau `$ dense-walk phase has $`\rho 1`$ throughout, whereas here $`\rho 1`$ except in the limit $`\tau \mathrm{}`$. To what extent this may affect the nature of the phase transitions is not clear; early results indicate that perhaps the dense-walk phase transition is not in the Ising universality class.

## References


text:
# Chapter 1 Shower Studies at around 10¹⁸ eV with the Surface Detector of the Pierre Auger Observatory

## Chapter 1 Shower Studies at around $`10^{18}`$ eV with the Surface Detector of the Pierre Auger Observatory

### Abstract

Three months of data taken by the Engineering Array Surface Detector of the Pierre Auger Observatory have been analyzed with the aim of setting a procedure for the study of showers at $`E_010^{18}`$ eV. We present the results, concerning: (i) the selection procedure and the check of its efficiency; (ii) the event rate and its stability; (iii) the characteristics of such events (i.e., angular direction, energy estimator $`S(1000)`$, rise time vs core distance, correlation with pressure) showing their consistency with the expectations. The performances of the whole Surface Detector of the Auger Observatory are discussed.

### 1. Introduction

The study of cosmic rays around $`10^{18}`$ eV with the surface detector of the Auger Observatory will be of main significance for: a) the verification and extension of the CR anisotropy measurement reported by AGASA , SUGAR  and Fly’s Eye  from the directions of the Galactic Center and Disk; b) the measurement of the CR flux in an energy range interconnected with other experiments; c) the possibility of improving the cross calibrations with the Auger fluorescence detector, which naturally operates at lower threshold. The Auger southern site surface detector, under construction in the province of Mendoza, Argentina, will consist of 1600 water Cherenkov detectors (WCD), on a triangular grid with a 1.5 km spacing, covering an area of about 3000 km2; a subset of it, the “Engineering Array” (EA), including 32 completely instrumented stations, has been in operation since 2002. We present here the analysis of three months of data taken between May and July 2002, to set a procedure for the selection of showers around $`10^{18}`$ eV.

### 2. Analysis and results

To define the most appropriate trigger configuration to select events with $`E_010^{18}`$ eV, we conducted a preliminary study, using the data from a WCD prototype operated in correlation with AGASA ($`E_0>10^{18}`$ eV and $`\theta <45^{}`$) . The analysis  consisted of the simulation of the Auger trigger and its application to the WCD data after proper conversion to the Auger scale. This allowed the derivation of the single tank trigger efficiency as a function of its rate, $`f_t`$, core distance, $`r`$, and energy, $`E_0`$. The array trigger efficiency, $`\eta `$, obtained from such data, is shown in Tab.1 for different $`f_t`$ and $`E_0`$.

Due to the low efficiency of 4-fold mode, 3-fold events have to be used. The background rate due to 3-fold random coincidences for the whole Auger detector is given in Table 2 together with the rate of expected CR events, $`N_{exp}`$, with $`10^{18}<E_0<210^{18}`$ eV, evaluated using the CR spectrum , the Auger acceptance ($`A\mathrm{\Omega }4700`$ sr km2, $`\theta <45^{}`$) and $`\eta `$.

The extraction algorithm is firstly based on the selection of 3 nearby triggered tanks in a equilateral pattern. For the Auger-EA ($`f_t20`$ Hz and coincidence window $`\mathrm{\Delta }t25\mu `$s) the measured rate, $`N_{sel}30`$/day, is consistent with that expected from chance coincidences.

The measured rate is then compared with the expected one due to CRs. Taking into account the contributions from different energies ($`2.510^{17}<E_0<410^{18}`$ eV), we evaluate it using spectrum , including the Auger-EA acceptance ($`A\mathrm{\Omega }51`$ sr km2, $`\theta <45^{}`$) and the array trigger efficiency $`\eta `$. $`\eta `$ is obtained using the WCD prototype data (converted to the Auger-EA scale), in correlation with the AGASA EAS, for the different energy intervals, with an extrapolation to $`E_0<10^{18}`$ eV. The total predicted rate amounts to $`N_{exp}=5\pm 1`$/day, consistent with the measured one, within the systematic uncertainties (we remind that the results are based on trigger efficiencies evaluated at the AGASA energy scale).

The extracted events have been reconstructed, following the prescriptions given in . The distributions of the arrival angles $`\theta `$ and $`\varphi `$ are given in fig. 2 (top left and right, respectively), showing good consistency with the expectations. The risetimes from 10% to 50% levels of the signal ($`t_{50}t_{10}`$) vs core distance ($`r`$) in four intervals of $`\theta `$ are shown in the same figure (bottom right), being in good agreement with simulated data  and real ones from Haverah Park . Concerning primary energy, an estimator, $`S(1000)`$ (i.e., the signal in Vertical Equivalent Muons at $`r=1000`$m) corrected for $`\theta `$, is shown in fig. 2 (bottom left). The conversion from $`S(1000)`$ to primary energy is subject of intense study by the Auger Collaboration (see e.g. ): currently, the peak of the $`S(1000)`$ distribution corresponds to primary energies ranging from 0.4 to 0.9$`10^{18}`$ eV (consistent with the present knowledge of the primary intensity at such energies, when compared with the measured rates). Lastly, the behavior of the event rate vs atmospheric pressure, $`P`$, has been studied. The pressure coefficient obtained is $`k=1/N(dN/dP)=0.011\pm 0.005`$ mbar-1 (significant at 2$`\sigma `$), in agreement with Haverah Park .

### 3. Conclusions

Three months of data taken by the Auger-SD EA in a stable configuration have been analyzed, with the following outcome:
\- an algorithm for the extraction of events around $`10^{18}`$ eV, based on geometry and single tank trigger condition corresponding to $`f_t`$2 Hz (to reduce the random coincidences to the percent level of the physical event rate) has been tuned and tested using a WCD prototype operated in connection with AGASA;
\- the event rate is consistent with the expectations from CR showers due to primaries of energy around $`10^{18}`$ eV;
\- the shower parameters (distributions of zenith and azimuthal angles, $`S(1000)`$, rising time vs core distance and vs zenith angle, correlation with pressure) are consistent with the expected ones;
\- assuming for the full Auger-SD $`\eta =0.23`$, and using the spectrum , we get $`N_{exp}(10^{18}\text{eV}<E_0<210^{18}\text{eV})2.810^4`$ ev/yr (in the AGASA energy scale). The detection level of the AGASA anisotropy ($`A4`$%), assuming poissonian fluctuations, would thus be $``$ 4.5 s.d./yr (in harmonic analysis). Concerning stability, representing a further requirement for the measurement, within the present statistics, non poissonian daily instability effects can be verified to be $`<`$ 6-7%.

### 4. References

1. Hayashida N. et al. 1999, Astropart. Phys. 10, 303

2. Bellido J.A. et al. 2001, Astropart. Phys. 15, 167

3. Bird D.J. et al. 1999, Ap. J. 511, 739

4. Sakaki N. et al. 1997, Proc. 25th ICRC 5, 217

5. Ghia P.L. and Navarra G. 2003, Auger-GAP note 2003-007

6. Nagano M. and Watson A.A. 2000, Rev. Mod. Phys. 72, 689

7. Billoir P. 2002, Auger-GAP note 2002-044, 2002-075

8. Sciutto S. 2002, private communication

9. Watson A.A. and Wilson J.G. 1974, J.Phys.A 7, 1199, Walker R. and Watson A.A. 1981, J.Phys.G 7, 1297

10. $`k=0.007\pm 0.001`$ mbar-1, Edge D.M. et al. 1978, J.Phys.G 4, 133


text:
# Dual kinetic balance approach to basis set expansions for the Dirac equation

## Abstract

A new approach to finite basis sets for the Dirac equation is developed. It solves the problem of spurious states and, as a result, improves the convergence properties of basis set calculations. The efficiency of the method is demonstrated for finite basis sets constructed from B splines by calculating the one-loop self-energy correction for a hydrogenlike ion.

At present, a great variety of calculations in atomic physics and quantum chemistry are based on finite basis sets. First successful attempts to utilize finite basis sets in relativistic quantum mechanics were made many years ago . Application of B splines for constructing basis sets  provided new impact to this field. Nowdays, B splines are widely employed in computational atomic and molecular physics (see  and references therein).

In contrast to the nonrelativistic case, the use of finite basis sets in the relativistic theory is generally accompanied by the occurrence of spurious states (see, e.g., ). For the attractive Coulomb potential, spurious states appear for $`\kappa >0`$ as the lowest bound states with non-physical energies ($`\kappa =(1)^{j+l+1/2}(j+1/2)`$ is the quantum number determined by the angular momentum and the parity of the state). For the point nucleus, their energies coincide with the physical energies of the lowest bound states with $`\kappa <0`$. The wave functions of these states oscillate rapidly and, therefore, in many cases they may be disregarded in practical atomic calculations. However, since the presence of the spurious states disturbs the spectrum, it worsens the convergence properties of the basis set calculations in some cases. For this reason, despite of the early promising results , the finite basis set method has not been extensively employed in calculations of radiative corrections. To date, most of these calculations are performed by means of analytical or numerical representations for the Coulomb-Green function (see, e.g., ) or by the space discretization method , in which the spurious states are eliminated from the very beginning (see  for details). The problem of spurious states is especially demanding in calculations of atoms in strong external fields and for molecules, where, generally speaking, they can not be selected as the lowest bound-state energies and, therefore, can not be eliminated without employing special methods. Furthermore, even in those cases where spurious states do not play any significant role, their presence generates difficulties for providing an adequate estimate of the accuracy of the calculations.

A number of proposals for solving the problem of spurious states were presented previously . Our analysis of these methods indicates, however, that none of them can be considered as completely satisfactory. We find that treating the problem by a particular choice of the boundary conditions in the B-spline method  does not provide any improvement in calculations of radiative corrections. The other methods either require considerable modifications of the standard numerical procedure or they are strongly limited to a specific choice of finite basis sets. Their applicability to calculations of radiative corrections thus remains questionable. For this reason, we suggest a new method for solving the problem of spurious states. The efficiency of the method is demonstrated for finite basis sets constructed from B splines by calculating the one-loop self-energy correction for a hydrogenlike ion.

For the case of a central field $`V(r)`$, the Dirac wave function is conveniently represented by

$`\psi (𝐫)={\displaystyle \frac{1}{r}}\left(\begin{array}{c}G(r)\mathrm{\Omega }_{\kappa m}(𝐧)\\ iF(r)\mathrm{\Omega }_{\kappa m}(𝐧)\end{array}\right),`$ (3)
where $`𝐧=𝐫/r`$. With this representation, the radial Dirac equations can be written as

$`H_\kappa \varphi =E\varphi ,`$ (4)
where (in units: $`\mathrm{}=1`$)

$`H_\kappa =\left(\begin{array}{cc}mc^2+V& c\left[{\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right]\\ c\left[{\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right]& mc^2+V\end{array}\right)`$ (7)
and

$`\varphi (r)=\left(\begin{array}{c}G(r)\\ F(r)\end{array}\right)`$ (10)
is the two-component radial wave function. The scalar product of the two-component functions is defined by

$`\varphi _a|\varphi _b={\displaystyle _0^{\mathrm{}}}𝑑r[G_a(r)G_b(r)+F_a(r)F_b(r)].`$ (11)
The radial Dirac equations can be derived from an action principle $`\delta S=0`$ with

$`S=\varphi |H_\kappa |\varphi E\varphi |\varphi ,`$ (12)
if proper boundary conditions for $`G(r)`$ and $`F(r)`$ are implemented. The functions $`\varphi (r)`$ can be approximated by

$`\varphi (r)={\displaystyle \underset{i=1}{\overset{2n}{}}}c_iu_i(r),`$ (13)
where the two-component functions $`u_i(r)`$ are assumed to be square integrable, linearly independent, and satisfying proper boundary condition at $`r=0`$. The variational principle reduces to the following algebraic equations

$`dS/dc_i=0,i=1,2,\mathrm{},2n.`$ (14)
This leads to a generalized eigenvalue problem

$`K_{ik}c_k=EB_{ik}c_k,`$ (15)
where $`K_{ik}=(u_i|H_\kappa |u_k+u_k|H_\kappa |u_i)/2`$, $`B_{ik}=u_i|u_k`$, and the summation over repeated indices is implied.

Let us first demonstrate that the widely applied choice

$`u_i(r)`$ $`=`$ $`\left(\begin{array}{c}\pi _i(r)\\ 0\end{array}\right),i=1,\mathrm{},n,`$ (18)
$`u_i(r)`$ $`=`$ $`\left(\begin{array}{c}0\\ \pi _{in}(r)\end{array}\right),i=n+1,\mathrm{},2n,`$ (21)
where $`\{\pi _i(r)\}_{i=1}^n`$ are square integrable functions satisfying the boundary condition $`\pi _i(0)=0`$, results in the occurrence of spurious states. In this case, equation (15) takes the form

$`(mc^2+VE)_{ik}p_k+cD_{ik}q_k`$ $`=`$ $`0,`$ (22)
$`c(D^{})_{ik}p_k+(mc^2+VE)_{ik}q_k`$ $`=`$ $`0,`$ (23)
where $`(\pm mc^2+VE)`$, $`D`$, and $`D^{}`$ are $`n\times n`$ matrices with elements

$`(\pm mc^2+VE)_{ik}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑r\pi _i(r)(\pm mc^2+VE)\pi _k(r),`$ (24)
$`D_{ik}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑r\pi _i(r)\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\pi _k(r),`$ (25)
$`(D^{})_{ik}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑r\pi _i(r)\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\pi _k(r),`$ (26)
and $`p_i=c_i`$, $`q_i=c_{i+n}`$ for $`i=1,2,\mathrm{},n`$. Let us consider the nonrelativistic limit ($`c\mathrm{}`$) and introduce vectors $`P`$ and $`Q`$ with components $`\{p_i\}_{i=1}^n`$ and $`\{q_i\}_{i=1}^n`$, respectively. Then Eq. (23) yields $`Q=(1/2mc)D^{}P`$. Substituting this expression into Eq. (22), we obtain

$`DD^{}P+2m(mc^2+VE)P=0.`$ (27)
For the pure Coulomb field, $`V(r)=\lambda /r`$ ($`\lambda >0`$), introducing the matrix $`C_{ik}=D_{ik}(m\lambda /\kappa )\delta _{ik}`$, Eq. (27) reduces to

$`C_\kappa C_\kappa ^{}P=ϵP,`$ (28)
where $`ϵ=[2m(Emc^2)+m^2\lambda ^2/\kappa ^2]`$ and the dependence of the $`C`$ matrix on $`\kappa `$ is explicitly indicated. Taking into account that $`C_\kappa =C_\kappa ^{}`$, we find that the corresponding equation for $`\kappa ^{}=\kappa `$ can be written as

$`C_\kappa ^{}C_\kappa P^{}=ϵ^{}P^{}.`$ (29)
Multiplying Eq. (28) with $`C_\kappa ^{}`$, we obtain that each nonzero eigenvalue of $`C_\kappa C_\kappa ^{}`$ is an eigenvalue of $`C_\kappa ^{}C_\kappa `$. Evidently, the inverse statement can be proven in a similar manner and the dimension of a nonzero eigenvalue subspace is the same for $`C_\kappa C_\kappa ^{}`$ and $`C_\kappa ^{}C_\kappa `$. Accordingly, the spectra of $`C_\kappa C_\kappa ^{}`$ and $`C_\kappa ^{}C_\kappa `$ may differ only by the dimension of the zero eigenvalue subspace. For finite matrices, the dimension of the subspace with $`ϵ=0`$ is the same for $`C_\kappa C_\kappa ^{}`$ and $`C_\kappa ^{}C_\kappa `$, since the total number of eigenvectors as well as the dimension of the nonzero eigenvalues subspace is the same for $`C_\kappa C_\kappa ^{}`$ and $`C_\kappa ^{}C_\kappa `$. Therefore, the finite matrices $`C_\kappa C_\kappa ^{}`$ and $`C_\kappa ^{}C_\kappa =C_\kappa C_\kappa ^{}`$ have an identical spectrum. Conversely, we know that the exact analytical solution of the Dirac equation for the Coulomb potential yields different lowest bound-state energies for $`\kappa <0`$ and $`\kappa >0`$. This is due to the fact that within the exact (infinite dimension) treatment the subspace with $`ϵ=0`$ may have different dimensions for $`\kappa <0`$ and $`\kappa >0`$ cases. This can easily be checked by solving the equation

$`(d/dr+\kappa /rm\lambda /\kappa )G(r)=0,`$ (30)
which in case of finite dimensions is equivalent to the equation $`C_\kappa ^{}P=0`$. Solving equation (30) yields $`G(r)=A_0r^\kappa \mathrm{exp}(m\lambda /\kappa r)`$. For $`\kappa <0`$, this solution has the proper behaviour at $`r0`$ and at $`r\mathrm{}`$. However, this does not hold for $`\kappa >0`$. Concluding, in the approximation of finite dimensions, our proof clearly indicates the presence of spurious states with energy $`Emc^2=m\lambda ^2/2\kappa ^2`$ corresponding to $`ϵ=0`$. It can be shown that this result remains valid for the full relativistic theory as well. It is obvious that spurious states must occur for any other potential one is dealing with in atomic calculations.

To eliminate the spurious states, in Ref.  ”kinetically balanced” Slater type functions were employed. Within this method, for $`\kappa >0`$ the lower components in equation (21) are replaced by functions $`\rho _i(r)`$ which, in the nonrelativistic limit, are related to the upper components $`\pi _i(r)`$ in equation (18) via

$`\rho _i(r)(1/2mc)(d/dr+\kappa /r)\pi _i(r).`$ (31)
This method provides a high accuracy in calculations of bound-state energies in atoms for both sign of $`\kappa `$, if an extended nuclear charge distribution is introduced and proper boundary conditions are implemented . However, since the basis functions are ”kinetically balanced” for positive energy states only, the application of this method to calculations of the QED corrections might be problematic . The equivalent treatment of positive and negative energy states would provide, in particular, the well-known symmetry properties of the Dirac spectrum under the transformations $`VV`$, $`\kappa \kappa `$, and $`GF`$. It is evident that the ”kinetically balanced” functions do not meet with this requirement.

In the original version of the B-spline method , to achieve that the first positive-energy states $`\kappa >0`$ correspond to physical bound states, an additional term had to be introduced in the Hamiltonian, which formally implements the so-called MIT boundary condition : $`G(R)=F(R)`$, where $`R`$ is the cavity radius, together with the condition $`G(0)=0`$. However, since the presence of the additional term does not imply any practical advantages, it is usually omitted in calculations. Instead, the boundary conditions are generally implemented by eliminating the first and the last basis function which are the only ones that do not vanish at $`r=0`$ and $`r=R`$, respectively. This method was successfully employed for calculations of the two-photon exchange diagrams within the rigorous QED approach  and for relativistic calculations of the recoil effect . However, its application to calculations of pure radiative corrections  was less successful, compared to the other methods . We conjecture that this would not be the case if the spurious states were eliminated in a more natural manner than it was done in .

It is known (see, e.g., ) that the case of the pure Coulomb potential requires generally a special care in implementing finite basis set methods. This is due to the singularity of the Coulomb potential at $`r0`$. However, in practical calculations it is standard to modify the potential to account for the finite nuclear size, which eliminates this problem. For this reason and for simplicity, we restrict our consideration to the case of a finite nuclear charge distribution, bearing in mind that the limit of a point nucleus can be treated by extrapolating a series of calculations for extended nuclei to vanishing nuclear size. For extended nuclei we propose to employ the following basis set

$`u_i(r)`$ $`=`$ $`\left(\begin{array}{c}\pi _i(r)\\ \text{}{\displaystyle \frac{1}{2mc}}\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\pi _i(r)\end{array}\right),in,`$ (34)
$`\text{}u_i(r)`$ $`=`$ $`\left(\begin{array}{c}{\displaystyle \frac{1}{2mc}}\left({\displaystyle \frac{d}{dr}}{\displaystyle \frac{\kappa }{r}}\right)\pi _{in}(r)\\ \text{}\pi _{in}(r)\end{array}\right),in+1,`$ (37)
where the linearly independent functions $`\{\pi _i(r)\}_{i=1}^n`$ are assumed to be square integrable and satisfying the proper boundary condition at $`r=0`$. We state that this basis set satisfies the following requirements:

1. It is symmetric with respect to the replacement $`\kappa \kappa `$ and the interchange of the upper and lower components.
2. The functions $`u_1`$, …, $`u_n`$ provide the correct relation between upper and lower components for $`|Emc^2|,|V(r)|2mc^2`$, while the functions $`u_{n+1}`$, …, $`u_{2n}`$ do the same for $`|E+mc^2|,|V(r)|2mc^2`$.
3. Calculations utilizing the standard finite basis set determined by Eqs. (18) and (21) can be easily adopted when employing the basis (34)-(37).
4. No spurious states occur for attractive as well as for repulsive potentials.

The properties 1 - 3 follow immediately from definitions (34) and (37). The absence of spurious states can be explained as follows. Performing similar steps as for the derivation of Eq. (27), for $`|Emc^2|2mc^2`$ we obtain

$`{\displaystyle \frac{1}{2m}}LP+(V+mc^2E)P=0,`$ (38)
where

$`L_{ik}={\displaystyle _0^{\mathrm{}}}𝑑r\pi _i(r)\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\pi _k(r).`$ (39)
Eq. (38) takes the form of the ordinary Schrödinger equation with $`l=|\kappa |+(\kappa /|\kappa |1)/2`$ in the finite basis representation. As is known, it generates no spurious states. The region $`|E+mc^2|2mc^2`$, where spurious states may exist for repulsive potentials and for $`\kappa <0`$, can be considered similarly. In this case, we obtain the equation

$`{\displaystyle \frac{1}{2m}}MQ+(V+mc^2+E)Q=0,`$ (40)
where

$`M_{ik}={\displaystyle _0^{\mathrm{}}}𝑑r\pi _i(r)\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)\pi _k(r).`$ (41)
Eq. (40) has also the form of the ordinary Schrödinger equation but with $`l^{}=|\kappa |(\kappa /|\kappa |+1)/2`$. It transforms into equation (38) under the replacements $`\kappa \kappa `$, $`VV`$, $`EE`$, $`QP`$ and does not generate any spurious states. This is a consequence of the equivalent treatment of the positive and negative energy states. For this reason, the new basis may be termed conventionally as dual kinetic-balance (DKB) basis.

The validity of statement 4 has also been proven by numerical calculations with $`\pi _i(r)=B_i(r)`$, where $`B_i(r)`$ are the B splines defined on the interval $`(0,R)`$ as in Ref. . The first and the last spline function have been omitted. Standard test calculations (see, e.g., ) show that this basis satisfies suitable completeness criteria as $`n\mathrm{}`$.

Finally, let us consider the calculation of the one-loop self-energy (SE) correction to the ground-state energy of a hydrogenlike ion employing the new basis set. Generally, the SE correction is expanded into the zero-, one-, and many-potential terms. The ultraviolet divergences in the zero- and one-potential terms and in the counterterm cancel each other and their evaluation can be performed according to the formulas presented in Refs. . As to the many-potential term, although it does not contain any ultraviolet divergences, its calculation is most difficult since it involves the summation over the whole Dirac-Coulomb spectrum. In Table 1, we compare our results obtained for the many-potential term for $`Z=20`$ employing the DKB basis set (34), (37) with $`\pi _i(r)=B_i(r)`$, the old basis (18), (21) with the same $`\pi _i(r)`$, and the results of a calculation using the analytical representation for the Coulomb-Green function. The shell model for the nuclear charge distribution has been used with the radius $`R=3.478`$ fm. In the basis set calculations, the contributions with $`|\kappa |10`$ were obtained by an extrapolation. Adding the zero- and one-potential terms to the many-potential term yields 0.06409 a.u. for the old basis, 0.06426 a.u. for the DKB basis, and 0.06425(1) a.u. for the Coulomb-Green function calculation. This comparison clearly demonstrates a significant improvement in accuracy, if the DKB basis is employed instead of the old one. More extensive calculations employing the DKB basis will be presented in forthcoming papers.

This work was supported in part by RFBR (Grant No. 01-02-17248), by the program ”Russian Universities” (Grant No. UR.01.01.072), and by the Russian Ministry of Education (Grant No. E02-3.1-49). The work of V.M.S. was supported by the Alexander von Humboldt Stiftung. V.A.Y. acknowledges the support of the foundation ”Dynasty” and the International Center for Fundamental Physics. G.P. and G.S. acknowledge financial support by the BMBF, DFG, and GSI.