text:
# The small-scale structure of the Magellanic Stream

## 1 Introduction

The Magellanic Stream (MS) is a thin, $``$ 10 wide (Putman & Gibson, 1999), tail of neutral hydrogen (HI), emanating from the Magellanic Clouds (MCs: the Large Magellanic Cloud, LMC, and the Small Magellanic Cloud, SMC) and trailing away for almost 100 across the sky ($`60`$$`<`$ Dec $`<+15`$). This huge HI structure is the most fascinating signature of the wild past interaction of our Galaxy with the MCs, and the MCs with each other. At the same time, streaming so far away from the MCs into the Galactic Halo, the MS is an excellent probe of the properties of the outer Galactic Halo, for which otherwise we have very few observational clues.

For many years the MS was viewed as a complex of six discrete concentrations (Mathewson et al., 1974). The new HI Parkes All-Sky Survey (Putman & Gibson 1999) reveals the more complex nature of the MS, with a fascinating network of filaments and clumps. No stars have been found so far in the MS, making it the closest medium to a primordial environment that exists in the Local Group today.

Hα emission appears to be routinely detected throughout the MS (Weiner et al., 2002), suggesting that the MS is being ionized by photons escaping from the Galaxy (Bland-Hawthorn & Maloney, 1999). While this mechanism is in accord with the H$`\alpha `$ observations of many high velocity clouds (HVCs), the predicted H$`\alpha `$ flux for the MS is significantly lower than what is found observationally (Weiner et al., 2002; Bland-Hawthorn & Maloney, 2002), suggesting that an additional source of ionization must be invoked. Several possible mechanisms have been suggested, including an interaction between the MS gas and the hot Galactic Halo gas, turbulent mixing, the existence of young massive stars embedded in the MS (Bland-Hawthorn & Putman, 2001), shocks, self-interaction of the MS gas, and magnetic fields (Konz et al., 2001).

The MS appears to be the result of interaction between the Galaxy and the MCs, but there is no consensus as to the exact form of this interaction (Putman, 2000b). Theories have swung back and forth on the relative importance of tidal stripping (Murai & Fujimoto, 1980; Gardiner & Noguchi, 1996) and various kinds of gas dynamical interactions (Mathewson et al., 1987). The recent discovery by Putman et al. (1998) of a counter HI stream, leading the MCs, lends weight to the tidal origin theory.

The work presented in this paper was motivated by the following two questions.

(1) Is the MS dissipating into the Halo?
The most fundamental issue about the origin and structure of the MS is to what extent interaction with the Galactic Halo determines or influences the MS gas. This problem becomes particularly important at the extreme northern end of the MS, because without pressure from an external medium the MS clouds should dissipate on a short time scale. However, if the external pressure is sufficient to confine the MS clouds, then the gaseous Galactic Halo has a much higher scale height than most theories predict. This paper addresses several issues in relation to this particular question.

(2) Is a hierarchy of structures present in this almost primordial environment?
In the Galaxy, a hierarchy of structure is present in the diffuse interstellar medium (ISM), down to very small scales. This is mainly ascribed to interstellar (IS) turbulence. This process is seen both in the ionized medium, traced by pulsar scintillation and dispersion variations, and in the atomic gas, traced by the spatial power spectrum of HI emission and absorption fluctuations (Fiedler et al., 1994; Green, 1993; Lazarian & Pogosyan, 2000; Dickey et al., 2001). Recently it has become possible to extend these studies to other galaxies, at least to our neighbors (Westpfahl et al., 1999; Stanimirovic et al., 1999; Elmegreen et al., 2001). The results suggest that the processes we see at work in the solar neighborhood are quite general. As the IS turbulence is driven, at least partially, by stellar processes such as winds and supernova remnants, it would be quite a surprise if the same turbulence spectrum was seen in an environment without any stars at all, like the MS. We have addressed this question briefly in Stanimirovic et al. (2002) and will concentrate on it further in a subsequent paper.

This paper is organized as follows: In Section 2 we summarize previous HI observations of the MS, and especially its northern part. As parameters such as the MS age and its distance come up often throughout this paper, a brief description of various theories for the MS formation is also outlined here. Section 3 describes new HI observations of two regions at the tip of the MS conducted with the Arecibo telescope. The HI data and observational results obtained are presented in Section 4. Several clumps and their properties are shown in Section 5. In Section 6 we discuss clump confinement issues and the constraints they place on the density of the Galactic Halo. Expected Halo densities, as well as several theoretical approaches for considering interactions between an MS cloud and the ambient medium, are discussed. An alternative explanation for the velocity field of one of the regions is also presented. We summarize our results in Section 7.

## 2 Previous HI observations and theoretical models

### 2.1 Lower resolution observations

First HI detections of the MS were obtained by Dieter (1965) and Wannier & Wrixon (1972), however Mathewson et al. (1974) were the first to discover a coherent large-scale structure which they named the Magellanic Stream. Observations with better resolution and sensitivity followed by Cohen & Davies (1975), Mirabel et al. (1979), Erkes et al. (1980), Mirabel (1981) and Mathewson & Ford (1984). These observations revealed the large-scale morphology of the MS, that is traditionally viewed as a long filament, containing a bead-like sequence of six discrete clouds: MS I, near the MCs at Dec $`60`$, to MS VI at the very northern tip at Dec $`+15`$ (Mathewson et al., 1977), see the top panel in Fig. 1. Mirabel (1981) in particular focused on the tip of the MS, searching for evidence of its disintegration. Wayte (1989) used the Parkes 64-m dish to uncover complex and chaotic motions right at the MS tip.

The most recent HI observations of the MS, with spatial resolution of 15 arcmin and velocity resolution of 26 km s-1, were undertaken with the Parkes Multibeam system (Putman, 2000b; Putman et al., 2002). These, as well as the recent observations obtained with the same instrument but with higher velocity resolution of 1 km s-1 by Brüns et al. (2000), revolutionized the classical picture of the MS. Instead of six concentrations of gas, two striking, long and distinct filaments are seen that run together for most of the length of the MS but merge in several places, giving the appearance of a double-helix structure. North of Dec $`0`$, the two filaments break into a network of clumps and smaller filaments (see the bottom panel in Fig. 1). A global velocity gradient is observed along the MS, from $`+350`$ km s-1 (LSR) near the MC’s to $`450`$ km s-1 (LSR) at the very northern tip.

### 2.2 Theoretical models

The two leading hypotheses for explaining the origin of the MS are based on tidal interaction versus ram pressure stripping, and are applied in several models. These two families of models predict different ages and distances for the MS. As the predicted parameters are important for the calculations in Sections 5 and 6 we carry out calculations for both sets of model predictions.

The tidal hypothesis (Murai & Fujimoto, 1980; Lin & Lynden-Bell, 1982; Gardiner et al., 1994; Gardiner & Noguchi, 1996; Yoshizawa & Noguchi, 1999) — invokes the idea of gravitational stripping of the SMC gas by the Galaxy. According to N-body simulations, the MS was formed 1.5 Gyr ago in a tidal encounter between the LMC and the SMC at their perigalactic passage. This encounter drew gas out of the SMC that later evolved into a leading bridge and a trailing tail. The MS follows the MCs on an almost polar orbit with a transverse velocity of 220 km s-1 (LSR). Simulations show that the MS gas occupies a range of distances, from 45 kpc in the vicinity of MS III, to 60 – 70 kpc at MS VI.

One of the main drawbacks of the tidal model used to be the observation of only one tail of gas behind the MCs. Recent Parkes observations (Putman et al., 1998) have discovered the ‘Leading Arm Feature’ (LAF), a counter-stream that leads the direction of motion of the MCs, as expected in the tidal models. The absence of stars in the MS is another troubling issue for tidal models, as both stars and gas should be affected equally by tides. However, the most recent N-body simulations by (Yoshizawa & Noguchi, 1999) show that a very compact initial configuration of the SMC stellar disk could result in only gas being disrupted, while the stars are left unaffected.

The ram-pressure hypothesis (Moore & Davis, 1994; Heller & Rohlfs, 1994; Sofue, 1994) — invokes the idea of ram-pressure stripping of gas from the MCs by an extended halo of diffuse gas around the Galaxy. Here, the MCs entered the extended ionized halo of the Galaxy some 500 Myr ago at a galactocentric distance of 65 kpc. Material was then stripped from the SMC and the intercloud region and began to fall toward the Galaxy. In this model, the gas with the lowest column density, that lost the most orbital angular momentum, has fallen the furthest toward the Galaxy. This gas is currently located at the tip of the MS, at a galactocentric radius of 25 kpc. The model requires that the diffuse Galactic Halo has densities of $`5\times 10^5`$ cm-3 at distance of 65 kpc in order to strip gas, but not decelerate it too much.

The tidal $`+`$ drag model (Gardiner 1999; Moore & Davis, in preparation) The recent observational confirmation of the LAF provides a new set of parameters for theoretical models. In order to better constrain the observational characteristics of the LAF, Gardiner (1999) introduced a weak non-gravitational term (a drag term) into the classical tidal model. The new simulations show interesting results: while the LAF is reproduced, simulated particles form a complete ring from the LAF, connecting the LAF to two trailing streams. These streams have different galactocentric velocities starting from MS I to MS V, and then meet together again around MS VI (see Fig. 4 from Gardiner 1999). This simulation predicts that large-scale velocity splitting should be seen throughout most of the MS. The success of this model, based on both gravitational and non-gravitational forces, encourages further developments in the field that are eagerly anticipated (Moore, in preparation; Maddison, Kawata, & Gibson 2002).

## 3 Observations and data reduction

The observations were obtained with the 305-m Arecibo telescope111The Arecibo Observatory is part of the National Astronomy and Ionosphere Center, which is operated by Cornell University under a cooperative agreement with the National Science Foundation. in June and July of 2000 and June of 2001. The Gregorian feed was used with the narrow 21-cm (L-band) receiver. The illuminated part of the 305-m dish covers an area of about 210$`\times `$240 m, giving a beam FWHM of approximately 3.1 $`\times `$ 3.4 arcmin (Heiles et al., 2001). Two observing bandwidths of 12.5 and 6.25 MHz were used simultaneously on two correlator boards, each with 1024 channels, for each circular polarization. The resultant velocity resolution is 2.6 and 1.3 km s-1, respectively.

Data were taken in the on-the-fly mapping mode, by driving in RA, at a rate of 4 arcsec per sec, and stepping in Dec by 2 arcmin between strips. The correlator was set to make one scan per strip by recording data continuously every 30 sec, which is effectively every 2 arcmin on the sky. As a result the data were slightly undersampled. At the beginning of each scan, a noise diode of known temperature was fired for 3 sec for calibration.

Two regions of interest were selected from Mirabel et al. (1979): one centered at the northern tip of the MS, $`(l,b)=(86.7,42.8)`$, and the other in the MS V region, $`(l,b)=(94.4,50.6)`$ (Fig. 1, these regions are outlined in the lower panel). Below we refer to these areas as MS VI and MS V, with the understanding that our maps do not cover the full extents of those clouds. Both regions were mosaiced with many small, overlapping maps, each of 10$`\times `$10 arcmin, which take $``$ 1 h to observe. The mosaic pieces were combined during the gridding process.

We removed a linear baseline by fitting ranges of emission-free channels on both sides of the line. The gain correction, as a function of both azimuth and zenith angle, was then applied. All spectra were subsequently convolved onto a sky grid using a Gaussian convolution function with FWHM of 2.6 arcmin, broadening the angular resolution to $``$ 4 arcmin. A gridding correction factor of 1.45 (Sramek & Schwab, 1989) was applied to account for the new beam size (S. Stanimirovic 2002, in preparation).

The relationship between flux density and brightness temperature is given by $`S(\mathrm{Jy}\mathrm{beam}^1)=0.11T_\mathrm{b}(\mathrm{K})`$. The noise level is 0.04 K per 1.3 km s-1 wide channel, which is equivalent to a column density of $``$ 1017 cm-2. The peak brightness temperature is 1.5 K in the MS VI region, and 1.1 K in MS V. For comparison, Mirabel et al. (1979) obtained a maximum column density of $``$ 4$`\times `$1019 cm-2 in the MS VI region using the Jodrell Bank telescope. In the same region we derive a maximum column density of 4.3$`\times `$1019 cm-2.

## 4 Observational Results

### 4.1 Morphology and velocity field of the MS VI region

The MS VI data is shown in Fig. 2 where several velocity channels through the spectral-line data cube are presented. Fig. 3 shows position-velocity cuts. Fig. 4 shows maps of the HI column density and the intensity-weighted velocity field (panels a and b). This region has several elongated clumps, inter-connected and forming a long curved filament. Often, sub-clumps are found nested inside larger clumps. We identify individual clumps and describe their properties in Section 5.1. Several interesting features appear in the data.

A large loop-like feature, centered at RA 23h 07m 26s, Dec 1247 30′′, with radius of 50 arcmin, is visible between $`v_{\mathrm{LSR}}=360`$ and $`v_{\mathrm{LSR}}=385`$ km s-1. The loop is present through $``$15 velocity channels but its size appears to stay the same from channel to channel. Although at first glance this loop is reminiscent of HI holes associated with expanding shells of gas seen in many galaxies, because no systematic change of the radius with velocity is found this cannot be interpreted as an expanding shell. Velocity profiles along the loop are broad but always single peaked.

It is interesting to note that the most prominent clumps have a dominant orientation, from south-west to north-east, with position angle of $`40`$ (measured counter-clock wise from north). This corresponds to a position angle of $`90`$ in the Magellanic coordinate system, defined by Wakker (2001) and is parallel to the main axis of the MS (Putman, 2000a). Similar orientation of clumps in MS VI was noted previously by Mirabel et al. (1979).

The total HI mass in the mapped region is 3.7$`\times `$ 105 M, if the distance of 60 kpc is assumed, and it comprises about 28% of the HI mass in the whole MS VI region estimated by Putman (2000a).

The velocity profiles in this region are single-peaked but quite broad with the dispersion being $``$ 15 km s-1. Several RA-velocity diagrams are shown in Fig. 3. At about Dec 1305 an interesting velocity gradient of almost 40 km s-1 (LSR) is obvious. This gradient stays almost the same in the galactocentric velocity system, meaning that a spatial variation of the projection of the LSR velocity relative to the Galactic Center cannot explain the velocity field in this relatively small, 40$`\times `$40 arcmin, field.

The velocity field shown in Fig. 4 is the first moment map of the HI cube. The two prominent features in this map are north-south ridge-lines at velocity $`365`$ to $`375`$ km s-1 in the east, and velocity $`395`$ to $`400`$ km s-1 in the west. The points labeled A and B in Fig. 4 indicate the extent of the east ridge which may continue beyond the edge of our map. Points labeled C, D and E indicate the extent of the west ridge. The steep velocity gradient evident in the middle panels of Fig. 3 is in the region bounded by points A, B, C and D in Fig. 4. This region alone resembles the velocity pattern of a rotating disk, however the velocity field of the larger area does not fit this interpretation. We discuss this further in Section 6.6.

The velocity field in Fig. 4 is most likely caused by a combination of blending of smaller clumps along the line of sight and an interaction between the MS clouds and the ambient medium. With the Arecibo beam we are able to resolve several clumps blended together and forming a velocity gradient, see Fig. 3, especially around Dec 13 08. Some of these clumps show tails in velocity and sharp edges. This may be due to interaction with the surrounding medium. We discuss clump morphology in Section 6.7. Mirabel et al. (1979) saw a similar phenomenon in another MS VI region, while Wayte (1989) attributed broad velocity profiles to the breaking of the MS clouds into sub-clouds with different velocities.

### 4.2 Morphology and velocity field of the MS V region

Velocity channel maps presented in Fig. 5 show the general morphology of the MS V region. Several dominant large-scale features are visible. In the $`320`$ km s-1 plane there is a round cloud centered around RA 23h 39m, Dec 0735 00′′. This region shows single-peaked velocity profiles. A long, straight filament stretches north-south for almost one degree, at RA 23h 38m 30s from Dec 8 to 9, at velocity of $`327`$ km s-1. At higher velocities of $`340`$ to $`355`$ km s-1, the most prominent feature is a curved filament, centered at RA 23h 39m 30s, Dec 0740. This consists of a core superposed on an arc that extends for almost 1.3.

North of Dec = 0720 all line profiles show double peaks. This is shown in the position-velocity diagrams, Fig. 6, where a single velocity component is present only in the upper left panel. All other panels show double profiles. The two velocity peaks are roughly centered at $`320`$ and $`360`$ km s-1 and run almost in parallel.

Double-peaked velocity structure has been seen previously in MS V by Mirabel et al. (1979) at very similar central velocities. Wayte (1989) also noted this velocity bifurcation in several regions in the MS. The most recent HIPASS data, with very coarse velocity resolution of about 26 km s-1, show that spatially different filaments are common throughout the whole MS and most likely have different origins (Putman, 2000a). Similar line-splitting is found in several compact HVCs by Braun & Burton (2000), who suggest that this structure is caused by organized outflows or inflows.

What causes the dual velocity structure seen in MS V? The velocity separation between the two peaks of $`40`$ km s-1 is reminiscent of velocity features seen in the SMC. This was formerly interpreted as a signature of the SMC being ripped into two clouds which began during its last close encounter with the LMC some 200 Myr ago (Mathewson & Ford 1984). Recently, high spatial and velocity resolution observations of the SMC (Staveley-Smith et al., 1997; Stanimirovic et al., 1999) discovered more complex velocity profiles and numerous expanding shells of gas. This line-splitting in the SMC is now attributed primarily to the effect of expanding shells. Multiple structures in velocity are also seen in the Magellanic Bridge region (McGee & Newton, 1986). No evidence for shells has been found in the MS.

We have fitted overlapping Gaussians to the blended velocity components, and the results show that their velocity separation is quite constant over the area of MS V, as seen in Fig. 6.

#### 4.2.1 Gaussian decomposition

To model the velocity profiles with one or a superposition of two Gaussian functions, the data cube was binned spatially in squares of 5 by 5 pixels. As an example, Fig. 7 shows one velocity profile, its decomposition into two Gaussian functions, and the residuals after subtracting the fit from the data. The mean rms value of the residuals in the whole region is close to the noise level, with its maximum below the 2-$`\sigma `$ level. The two separate components are centered around $`350`$ km s-1, the ‘higher’ velocity component (‘H’), and $`320`$ km s-1, the ‘lower’ velocity component (‘L’). The integrated intensity for each velocity component is shown in Fig. 8. The higher velocity component is the core-arc feature described above, plus a diffuse feature to the north-west. The lower velocity component is the round cloud and the north-south filament. The HI mass of the higher velocity component is 2.3$`\times 10^5`$ M. The HI mass of the lower velocity component is 3.2$`\times 10^5`$ M (if the distance of 60 kpc is assumed). These together make up only 15% of the HI mass of the entire MS V region, 3.9$`\times 10^6`$ M (Putman, 2000a).

A slight gradient in velocity of about 10 km s-1 is seen from the north to the south and is shown by the slight tilt of the ridge lines in Fig. 6. These trace roughly $`360`$ to $`350`$ km s-1 for the higher velocity component, and $`325`$ to $`315`$ km s-1 for the lower velocity component, from right to left on the lower middle panels of Fig. 6. The dispersion of the lower velocity component has everywhere a value of about 10 km s-1. The dispersion of the higher velocity component shows a significant variation with position. In the north, the velocity profiles are very narrow, with typical velocity dispersion of about 6 km s-1. The dispersion increases towards the south reaching about 13 km s-1 along the core-arc feature.

We have experimented with fitting more than two velocity components to the data, however the resultant fits did not improve the residuals significantly. Neither MS V nor MS VI regions show narrow line components indicative of cold HI (with temperature of 30 – 300 K). This is in agreement with theoretical predictions by Wolfire et al. (1995), that no cold cores should be found in clouds in the Galactic Halo at heights above 20 kpc.

## 5 Clump analyses

### 5.1 Clump identification

In order to study the properties of the MS clumps and the indirect clues they can provide about the properties of the Galactic Halo, we have attempted to identify individual clumps. We used two different approaches.

(i) To identify smaller clumps and sub-clumps within larger clumps, the three-dimensional clump-finding algorithm, ‘clfind’ coded in IDL (Williams et al., 1994) was applied to both data sets. A spatial binning by 4 pixels, and a velocity binning by 2 pixels was applied to reduce the computing time. ‘Clfind’ views data as a set of contour levels. It locates local maxima and defines clumps by tracing them by connecting pixels, at each contour level, that are within one resolution element of each other. The algorithm does not assume any particular geometry for the clumps. The most important input parameter for ‘clfind’ is the contour increment. This was set to twice the noise level, as recommended in Williams et al. (1994). All selected clumps were further vigorously inspected by eye by overlaying the ‘clfind’ output on the original data using the karma visualization package (Gooch, 1995). Only a few (seven) clumps that are isolated and fully inside the edges of our map were considered for further analysis. Those clumps, three in MS V and four in MS VI, are shown in Figs. 9 and 10. For each clump the central position, central velocity, clump size ($`R_\mathrm{c}`$) measured from the integrated intensity map, FWHM velocity line width ($`\mathrm{\Delta }v`$), and peak brightness temperature were estimated and are listed in Table 1.

(ii) The clumps found by ‘clfind’ may or may not be distinct physical structures. On larger scales, all of MS VI and the ‘L’ and ‘H’ components of MS V could be considered as clumps themselves. The corresponding quantities for these three largest clumps are also listed in Table 1.

### 5.2 Observed and derived properties

Fig. 11 shows the angular size and line width for all the clumps. A correlation between line width and cloud size has been found for interstellar clouds in a number of surveys (references in Vazquez-Semadeni et al. 1997). This is one of the famous Larson scaling laws: $`\mathrm{\Delta }vR_\mathrm{c}^\beta `$. For the case of self-gravitating molecular clouds $`\beta =0.5`$ (Myers & Goodman, 1988), a slightly flatter slope of 0.4 was found by Falgarone et al. (1992) in molecular clouds that are very little contaminated by star formation. We overlay the Larson relation in Fig. 11 for the assumed distance of 60 kpc. While all sub-clumps and one of the largest clumps (the whole of MS VI) appear to follow the correlation, two other large clumps have $`\mathrm{\Delta }v`$ lower by about 10 km s-1 than the predicted values.

We now derive several important clump properties. These properties depend on the assumed distance of the MS. Figs. 12 through 16 have double axes corresponding to the extreme distances of 60 and 20 kpc (see Section 2.2). For all these figures, the left-hand and bottom axes correspond to the far distance, and the right-hand and top axes correspond to the near distance.

The total HI mass is derived by integrating the brightness temperature over the entire line width. The HI mass is plotted as a function of cloud size in Fig. 12. For a given size our noise level sets a minimum detectable cloud mass. Using the narrowest measured clump line widths of $``$ 10 km s-1 we get the threshold indicated by the dashed line on Fig. 12. Fig. 12 suggests that there may be many more structures like our clumps in MS V and MS VI which could be found by a more sensitive survey.

Assuming that clouds are spherically symmetric, using their HI mass and radius we can estimate an average HI volume density, $`n_\mathrm{c}`$. We find the values of $`n_\mathrm{c}`$ shown in Fig. 13, where the 3-$`\sigma `$ detection threshold is again shown by the dashed line. A typical value for the small clumps is 0.05 cm-3 and for the large clumps 0.02 cm-3, for the 60 kpc distance assumption. For the 20 kpc distance the volume density increases by a factor of three, as shown by the right-hand scale on Fig. 13.

## 6 Discussion

### 6.1 Cloud confinement

We now turn to cloud stability issues and the survival of the MS so far away from the MCs. There are several questions to ask about the HI structures which we find.

(1) Are the MS clumps gravitationally confined?
In order to be gravitationally confined, a clump’s escape velocity must be larger than $`\mathrm{\Delta }v`$/2. This requires that the clump’s total mass is $`M_{\mathrm{grav}}R\mathrm{\Delta }v^2/8G`$. Fig. 14 compares the measured HI masses of our clumps with this hypothetical gravitational mass. For all our points $`30<M_{\mathrm{grav}}/M_{\mathrm{HI}}<250`$, if the distance of 60 kpc is assumed, and $`120<M_{\mathrm{grav}}/M_{\mathrm{HI}}<750`$, for the 20 kpc distance assumption. Thus gravitational confinement would require unreasonable amounts of dark matter.

(2) Are the MS clumps unbound and freely expanding?
In the tidal model the MS age is 1.5 Gyr, while in the ram pressure model the age is 500 Myr. If the clouds are freely expanding their expansion time is given by: $`t_{\mathrm{exp}}=2R/\mathrm{\Delta }v`$. These numbers are much less than the MS age. Fig. 15 plots the expansion age as a function of clump radius. The expansion age as a fraction of the MS age is indicated by the horizontal lines which apply to either model. This huge discrepancy between the expansion time and the age is the fundamental problem of the gas dynamics of the MS, as pointed out first by Mirabel et al. (1979). Barring the unreasonable dark matter halos required for gravitational confinement the simplest explanation is that the MS clouds are confined by a hot gaseous Galactic Halo.

### 6.2 Density of the Milky Way Halo

If the MS is confined by the Galactic Halo, then the external pressure, $`P_\mathrm{h}=kn_\mathrm{h}T_\mathrm{h}`$, must equal the clouds’ internal pressure $`P_\mathrm{c}=kn_\mathrm{c}T_\mathrm{c}`$: $`P_\mathrm{c}P_\mathrm{h}`$. The relevant value of $`T_\mathrm{c}`$ is not the microscopic kinetic temperature but the total random motion of the atoms which we describe as a temperature defined by the line-width:

$$T_\mathrm{c}=\frac{\mathrm{\Delta }v^2}{8\mathrm{ln}2}\frac{k}{m_\mathrm{H}}=21.9\mathrm{K}\times \left(\frac{\mathrm{\Delta }v}{\mathrm{km}\mathrm{s}^1}\right)^2.$$
(1)
Assuming a halo temperature of 106 K and using $`\mathrm{\Delta }v`$ for the clumps in Table 2 and the densities $`n_\mathrm{c}`$ illustrated on Fig. 13, we get Halo densities $`n_\mathrm{h}`$ shown in Fig. 16 as a function of distance from the Galactic Plane, $`z`$. Typical values for the near distance ($`z15`$ kpc) are $`n_\mathrm{h}=10^3`$ cm-3 and for the far distance ($`z45`$ kpc) are $`n_\mathrm{h}=3\times 10^4`$ cm-3. Note that these estimates are upper limits on the Halo density as magnetic, turbulent or ram pressure contributions were not included.

### 6.3 Expected Halo density

Current theories of galaxy formation predict the existence of a hot Galactic Halo containing the leftover gas which was unable to cool since the formation of the Galaxy. However, the density structure of the Halo is still very uncertain, especially at large distances from the Galactic Plane. Here we compare our results with several previous estimates, both observational and theoretical.

#### 6.3.1 Observational constraints

The only direct observations of the hot Halo gas come from measurements of the soft X-ray background (Snowden et al., 1997). These X-ray maps reveal hot gas with a temperature $`10^{6.6}`$ K, a scale-height of 1.9 kpc and an electron density of $``$ 3.5$`\times 10^3`$ cm-3. Indirect observational constraints on the Halo density are more common. These are based on dispersion measure measurements of pulsars in globular clusters, or on the observational properties of HVCs, including the MS clouds.

Mirabel et al. (1979) estimate that $`n_\mathrm{h}=7\times 10^4`$$`7\times 10^5`$ cm-3 for the distance range 10 – 100 kpc from HI observations of the MS and under similar confinement considerations as discussed in this paper. Similar values, $`n_\mathrm{h}10^4`$ cm-3, were suggested from Hα observations along the MS by Weiner & Williams (1996). Quilis & Moore (2001) applied the terminal velocity model by Benjamin & Danly (1997) for HVCs moving at $`300`$ km s-1, at distances 10 – 100 kpc and with HI column density of $`10^{20}`$ cm-2. This results in the requirement that $`n_\mathrm{h}(10\mathrm{kpc})<3\times 10^3`$ cm-3 and $`n_\mathrm{h}(100\mathrm{kpc})<3\times 10^4`$ cm-3, under the assumption that these are pure gas clouds. Our results agree well with the estimates from these three studies.

Another constraint for the Halo density at the distance of 130 kpc of $`10^5`$ cm-3 comes from HI observations of HVC 125$`+`$41$``$207 by Brüns et al. (2001). A significantly lower value for the upper limit on the Halo density, $`n_\mathrm{h}(50\mathrm{kpc})<10^5`$ cm-3, was imposed by Murali (2000) by allowing the whole MS VI cloud to survive at the tip of the MS for about 500 Myr. We look more closely into this approach in Section 6.4.

#### 6.3.2 Theoretical predictions

Benjamin & Danly (1997) considers a Galactic Halo consisting of three components: the warm ionized layer of HII parameterized by Reynolds (1993); the mean HI density distribution given by Dickey & Lockman (1990); and the hot Halo component with $`T_\mathrm{h}=10^6`$ K prescribed by Wolfire et al. (1995). While the first two components were parameterized observationally, the third component was constructed theoretically assuming an idealized isothermal halo in hydrostatic equilibrium, with the mid-plane density matching X-ray emission data by Garmire et al. (1992). The combination of all three components results in the analytic representation for $`n_\mathrm{h}`$ as a function of $`z`$ given by Benjamin & Danly (1997):

$$n_\mathrm{h}(z)=1.1\times 10^3\left[1+\left(z_{\mathrm{kpc}}^2/19.6\right)\right]^{1.35}\mathrm{cm}^3.$$
(2)
At $`z=15`$ kpc the expected density is $`n_\mathrm{h}=4\times 10^5`$ cm-3, and at $`z=45`$ kpc $`n_\mathrm{h}=2\times 10^6`$ cm-3. These values are 25 – 150 times lower than those found in Section 6.2.

The ram-pressure model for the MS origin by Moore & Davis (1994) requires slightly higher values for the Halo density in order to match the kinematics of MS clouds:

$$n_\mathrm{h}=0.002/\left[1+\left(R/12\mathrm{kpc}\right)\right]^2\mathrm{cm}^3.$$
(3)
Hence, at a distance of 20 kpc this yields a gaseous halo with $`n_\mathrm{h}=3\times 10^4`$ cm-3, while at a distance of 60 kpc $`n_\mathrm{h}=6\times 10^5`$ cm-3. Our estimates are much closer to the values predicted by this model but still several times higher.

### 6.4 Nature of Galactic infall

There are (at least) three different approaches present in the literature concerning the dynamics of clouds moving through an ambient medium. As understanding of the nature of this interaction can in turn provide clues about the properties of the Halo, here we discuss these approaches briefly.

(1) The mass loading approach
As the MS moves through the hot ambient Halo gas its clouds are subject to a drag force. Murali (2000) showed that this motion is actually dominated not by drag but by strong heating from accretion. If the Halo density is too high, heating due to accretion can cause cloud evaporation. This process is usually referred to as mass-loaded flow and it is common in comet interactions with the solar wind (Flammer, 1991). The Halo gas ionizes the cloud’s leading edges and picks up the ‘new material’ in a very asymmetric way that is constrained by the magnetic field: magnetic field lines slip past the cloud therefore allowing ablation of material only from the poles. Murali (2000) accounted for these effects and estimated that for the survival of the whole MS IV region for 500 Myr, at a distance of 50 kpc, $`n_\mathrm{h}<10^5`$ cm-3 is required. The key factor in this approach is the cooling resulting from a cloud’s mass-loss.

Following Murali’s calculations (their equation 1) for the MS V and MS VI clouds studied here, drag and accretion set a very low upper limit for the Halo gas density: $`n_\mathrm{h}10^7`$ cm-3 at $`z=45`$ kpc, and $`n_\mathrm{h}10^6`$ cm-3 at $`z=15`$ kpc, if mass-loss is not included. These values are significantly lower than theoretical predictions. Unless a significant cooling due to mass-loss is introduced, it is very hard for the hot Halo gas to exist in this model.

If the mass-loss rate, Ṁ$`=2.5\times 10^2`$ M yr-1 (Murali 2000) is included, the requirements are comparable to our estimates: $`n_\mathrm{h}3\times 10^4`$ cm-3 at $`z=45`$ kpc and $`n_\mathrm{h}3\times 10^3`$ cm-3 at $`z=15`$ kpc. A slightly lower rate of Ṁ$`=10^3`$ M yr-1 results in ten times lower Halo densities. However, the mass-loss rates of Ṁ$`=10^3`$ or few$`\times 10^2`$ M yr-1 throughout the whole cloud life-time are too high in comparison with the cloud HI masses discussed in Section 5.2. Such mass-loss rates would imply that cloud HI masses were 10 – 1000 times higher at the time of the MS formation. Smaller rates that are in agreement with the current HI masses are not sufficient to allow the existence of the Halo gas. Detailed simulations of the mass-loaded mechanism may help to constrain parameters such as mass-loss rate.

(2) Ballistic approach
Gregori et al. (1999) showed results of a 3-D study of a ballistic interaction of a moderately supersonic dense cloud with a warm magnetized medium. Since clouds are supersonic, their motion leads to the formation of a forward bow shock and a reverse crushing shock propagating through the cloud. As the cloud moves, its surface is subject to several instability mechanisms, the most disruptive of these being the Rayleigh-Taylor instability (R-T) that develops at the interface between two fluids when the lighter fluid accelerates the heavier one. These instabilities slowly disrupt the entire cloud and can dramatically change its morphology. In contrast to the mass-loaded flows, magnetic field lines here stay trapped, causing the development of strong magnetic pressure at the leading edge of the cloud.

The time-scale over which this dramatic change occurs is several times the so called ‘crushing time’, $`\tau _{\mathrm{CR}}=2R_\mathrm{c}\chi ^{1/2}/Mc_\mathrm{s}`$, with $`M`$ being the Mach number, $`\chi `$ being the ratio of cloud to ambient density, and $`c_\mathrm{s}`$ being the sound speed. During a time shorter than $`\tau _{\mathrm{CR}}`$, the magnetic pressure is comparable to the ram pressure, while afterwards ram pressure dominates. For the MS clouds, $`M2`$, and assuming $`\chi =100`$ from Section 6.2, $`\tau _{\mathrm{CR}}10`$ – 30 Myr, depending on whether the near or far distance to the clouds is assumed. As $`\mathrm{a}\mathrm{few}\times \tau _{\mathrm{CR}}`$ is a significant fraction of the MS age this suggests that during a large fraction of a cloud’s life-time magnetic pressure may play a significant role. Inclusion of magnetic pressure in the calculations in Section 6.2 may somewhat reduce the required values for the Halo density.

(3) Hydrodynamical approach
Full hydrodynamical, 3-D simulations of HVCs moving through a diffuse hot gaseous component were performed recently by Quilis & Moore (2001). They investigated both pure gas clouds, being in pressure equilibrium with the external hot medium, and dark matter dominated clouds with an additional potential field that maintains their dynamical equilibrium. Particular cometary morphology, seen in the case of many HVCs as well as in many MS clouds, was reproduced for both cloud types but under the condition that the diffuse medium has density $`n_\mathrm{h}>10^4`$ cm-3.

### 6.5 Terminal velocity and neutral gas fraction

At the tip of the MS, clouds have had plenty of time to decelerate by a drag force (Murali 2000). The present velocity of the MS clouds is moderately supersonic (220 km s-1 vs. sound speed cs of $`100`$ km s-1) suggesting that clouds are most likely now at the terminal velocity (Benjamin and Danly 1997). This means that clouds have stopped decelerating and that their boundary layers are subject to Kelvin-Helmholtz instabilities, while Rayleigh-Taylor instabilities are not important at this stage. Benjamin & Danly (1997) show that for a population of clouds with known distances, terminal velocities and column densities, it is possible to constrain the mean density of the gaseous halo and cloud neutral fraction $`f_\mathrm{c}=N_{\mathrm{HI}}/(N_{\mathrm{HI}}+N_{\mathrm{HII}})`$:

$$C_\mathrm{D}f_\mathrm{c}n_\mathrm{h}=\frac{2N_{\mathrm{HI}}g(z)}{v_\mathrm{T}^2},$$
(4)
here $`C_\mathrm{D}`$ is the drag coefficient, $`g(z)`$ is the gravitational acceleration, and $`v_\mathrm{T}`$ is the cloud’s terminal velocity. The mean column density for the MS V and MS VI clouds is $`\mathrm{log}N_{\mathrm{HI}}19`$ cm-3. Assuming that $`v_\mathrm{T}=220`$ km s-1, the right-hand side of the equation does not depend on distance and has a mean value of $`C_\mathrm{D}f_\mathrm{c}\rho _\mathrm{h}6\times 10^4`$ cm-3. If we assume $`C_\mathrm{D}=1`$ and apply Halo densities from Section 6.2, an estimate of the cloud’s neutral fraction can be obtained: at $`z=15`$ kpc $`f_\mathrm{c}=0.6`$, while at $`z=45`$ kpc $`f_\mathrm{c}1`$.

### 6.6 An alternative, gravitational confinement

An alternative explanation for cloud confinement is that they are gravitationally bound. This paradigm has been proposed for several compact HVCs (CHVCs) which are in some respects similar to MS V and MS VI clouds. An additional hint for this is provided by the peculiar velocity field seen in the MS VI region (Figs. 3 and 4), and discussed in Section 4.1, that may be interpreted as being due to differential rotation. In the case of three CHVCs, Braun & Burton (2000) found evidence for such orbital motion. In order to compare the MS VI cloud with these CHVCs we have attempted to model the velocity field north of Dec 1230 using the tilted ring algorithm. The derived rotation curve has a quick and almost linear rise to $``$ 15 km s-1, and then flattens out with large uncertainties.

This curve was used to obtain an estimate of the dynamical mass: $`M_{\mathrm{dyn}}=1.5\times 10^7`$ M or $`M_{\mathrm{dyn}}=5\times 10^6`$ M, depending whether the far or near distance is assumed. This further implies dark-to-visible mass ratios of 50 to 180, respectively, after including a contribution by helium of 40% by mass to the observed HI mass. Therefore, if interpreted as a rotating HI cloud, this would require a very dark matter halo.

No model for the formation of the MS has ever predicted the existence of dark-matter dominated clouds in this gas filament. Even the parent galaxies, the MCs, are not rich in dark matter. It is therefore hard to believe that rotational support is realistic. This was the conclusion of Section 6.1. However the derived rotation curve is very similar to the rotation curves of CHVCs derived by Braun & Burton (2000). Dynamical masses of these CHVCs are about 10 times higher than the dynamical mass of the MS VI region, while the dark-to-visible mass ratio of the MS VI region is a few times higher than that for CHVCs. As the velocity field of MS VI is most likely the result of clump blending along the line-of-sight, it is interesting to note how well this effect can mimic differential rotation.

### 6.7 On the clump morphology

Several clumps in the MS data cube show gradients in the brightness and velocity distribution, exhibiting a cometary or head-tail morphology. The best example is the south-eastern side of the MS V cloud at $`340`$ to $`355`$ km s-1, see Fig. 17. This bow-shock-like feature may be a strong compression front, centered at RA 23h 40m, Dec 0734 with two tails swept back on either side. This is strongly suggestive of a supersonic interaction between the HI cloud and an external, low density, ionized medium.

Similar cometary morphology is seen in many HVCs in surveys of Meyerdierks (1991) and Brüns et al. (2000). Putman (2000a) found many such clouds in the MS, which she speculated may be in the process of evaporation. The 3-D hydrodynamical simulations of an HVC moving through a diffuse hot medium by Quilis & Moore (2001) show a bow shock in front of the cloud, gas compression at the front edge of the cloud, and visible or invisible tails behind the cloud, depending on the density of the ambient medium. The double-tail morphology seen here is more similar to the simulations of the interaction between a disk-like structure with the surrounding medium by Quilis et al. (2000). These simulations were the first successful attempt to incorporate complex turbulent and viscous stripping at the interface of the cold and hot gaseous components. They predict a compression front and two tails folding on both sides. The similarity of the bow-shock-like feature with the simulations by Quilis et al. (2000) demonstrates the importance of turbulent and viscous mixing at the boundary layers when considering cloud – ambient medium interactions and should be fully explored in the future.

## 7 Conclusions

In this paper, we have presented new HI observations of two regions, MS V and MS VI, at the northern tip of the MS using the Arecibo telescope. The new data sets show the complex morphology of the MS with numerous inter-connected clumps. We have found double velocity structures in the MS V region. This region also contains an interesting bow-shock-like feature that is strongly suggestive of an interaction between the MS and an external medium. The MS VI data set shows a large velocity gradient, caused most likely by blending of small clumps along the line of sight. This region resembles the velocity field of a rotating disk, but we find that rotation is not a likely explanation.

Several MS clumps have been isolated to investigate their confinement mechanism. We show that unreasonably large amounts of dark matter are required in order for clumps to be gravitationally confined. Clumps do not seem to be in free expansion either. The easiest way to explain the clump properties is with external pressure confinement by the hot Galactic Halo. In this scenario we place an upper limit on the Halo density: $`n_\mathrm{h}=10^3`$ cm-3 at $`z15`$ kpc, and $`n_\mathrm{h}=3\times 10^4`$ cm-3 at $`z45`$ kpc. These values agree well with several previous indirect observations of the Halo, but are significantly higher than the values predicted theoretically for an isothermal hot Halo by Wolfire et al. (1995). Our results are closer to the expected Halo densities of the Moore & Davis model.

Although cloud evaporation is most likely an important process at the tip of the MS, the mass-loaded approach for cloud interaction with the ambient medium requires very high mass-loss rates to enable the existence of the hot Halo gas. The ballistic consideration of clouds in the Halo suggests that during a large fraction of the cloud lifetime magnetic pressure may have a significant role in cloud evolution. Inclusion of magnetic pressure would reduce somewhat the required values for the Halo density.

Hydrodynamical modeling of cloud interaction with the ambient medium also requires high Halo densities to reproduce cometary features commonly seen in HI observations of HVCs. If we assume that MS clouds are at their terminal velocity, the cloud’s neutral fraction can be estimated at $`z=15`$ kpc $`f_\mathrm{c}=0.6`$, while at $`z=45`$ kpc $`f_\mathrm{c}1`$.

We wish to thank several colleagues for their valuable input and discussions. We are in debt to Phil Perillat and Chris Salter for helpful suggestions during observations and data reduction. We are grateful to Abby Hedden for eagerly helping with observations and to Amanda Kirchner for assisting during the June 2000 observing run. We thank Mary Putman for stimulating discussions and for providing the HIPASS MS data at the northern tip prior to publication. Discussions with Tom Jones, Bon-Chul Koo, Robert Braun and Butler Burton are greatly appreciated. We also thank Chris Salter and Lister Staveley-Smith for a careful reading of the paper which resulted in many improvements.


format: text

text:
# Tight-binding parameters from the full-potential linear muffin-tin orbital method: A feasibility study on NiAl

## I Introduction

There has been a growing interest in recent years to construct methods of electronic-structure calculation in which the required computational time scales linearly with the size of the system, the so-called $`O(N)`$ method. The proposed approaches abandon the k-space method used in convential band-structure calculations in favor of working directly in real space.  Instead of calculating the energy bands, one instead focuses on local properties, e.g., the local density of states. Extensive properties, like the total energy, are obtained by integrating over space instead of over the Brillouin zone (hence the linear scaling with the size of the system). This real-space approach is particularly relevant to cases in which lack of perfect crystalline symmetry is an essential part of the problem (e.g., in surfaces, impurities, defects, and amorphous materials) since in those cases Brillouin zone and energy bands are simply non-existent. The tight-binding (TB) method is among the most popular implementations of this real-space strategy.

The tight-binding approach to electronic structure has a long history which dates back to the classic work of Slater and Koster half a century ago.  It has been found especially useful in the modeling of the important technological material silicon. The conventional means to obtain the TB parameters is to follow the original procedure of Slater and Koster namely to fit the TB energy bands to the ones obtained from accurate self-consistent calculations. The parameters for elemental solids obtained in this way have been compiled by Papaconstantopoulos. This database has also been actively enlarged to include parameters for various alloys. Recent interest in real-space electronic structure methods has supplied a renewed vigor to this field of TB parametrization.

Despite its seemingly simple principle, the fitting process that one uses in obtaining the TB parameters is not always straightforward in practice. Starting from a guess set of TB parameters one calculates the TB energy bands which are then compared with the accurate bands from a self-consistent calculation. One can proceed to minimize the merit function by using standard nonlinear optimization procedures. This, in general, is quite a straightforward procedure for simple systems with one atom per unit cell. However, for more complicated multiatom unit cells with more intricate “spaghetti” of bands, the number of independent parameters on which the merit function depends grows rapidly. In the optimization process for this case, the merit function can get trapped easily in a local minimum and the output TB bands have little resemblance to the bands that they are intended to fit. Getting around this problem by choosing a good initial set of parameters, unfortunately, still comprises more art than science. Another drawback of the TB parameters that are obtained from fitting is that one normally fits only the bands below the Fermi energy and a few bands above it. In general the parameters obtained depend on the number of bands that are used in the fitting. This feature makes it difficult to attach a physical meaning to the parameters.

Andersen and coworkers have proposed a method to obtain TB parameters directly from his linear muffin-tin orbital (LMTO) method. This method has also been extended by other groups. The method outputs the orthogonal hamiltonian (hopping) parameters between highly-localized (screened) muffin-tin orbitals in the two-center approximation. Since each hopping parameter effectively vanishes beyond second-nearest neighbors, TB calculations using these parameters can be highly efficient. Moreover, the obtained parameters are also independent of structure (transferable), a property which is highly desirable in TB applications. Most of these nice properties, however, depend on the atomic sphere approximation (ASA) which is used to derive these results. This approximation has been known to work rather well for close-packed structures but is not very accurate for open structures, although this can be remedied to some extent by using the so-called empty spheres. The TB-LMTO method assigns a minimal base (at most one s, three p, and five d orbitals per atom) to each atom in the crystal. In the LMTO-ASA method, the muffin-tin orbitals are set to have zero kinetic energy ($`\kappa ^2=0`$) in the interstitial region.

Full-potential (FP) implementations of linear electronic-structure methods, e.g., the FPLMTO method or the FP Linear Augmented Plane Wave (FPLAPW) method, are among the most accurate electronic-structure methods. It is therefore desirable to have a direct method of extracting practical TB parameters from these methods. The basis functions in the FPLAPW method are solutions to Schrödinger equation inside the MT sphere and plane waves in the interstitial region. Since the basis functions are not localized, its matrix elements therefore cannot be mapped into the TB form. (It should be noted, however, that it is possible to project a FPLAPW calculation into localized basis.) The basis functions in the FPLMTO method, on the other hand, can be chosen (by setting the value of the kinetic energy $`\kappa ^2`$) to have a decaying tail outside its MT sphere. Direct mapping to TB parameters is therefore possible for the FPLMTO method. To our knowledge this direct approach has been seriously studied only by McMahan and Klepeis who used it to study Si and Si/B phases.

In this paper we use this direct approach and extract TB (hamiltonian and overlap) parameters for the intermetallic alloy NiAl from a FPLMTO method. These parameters are then used as input to a real-space calculation of the local density of states using the recursion method. To our knowledge, this is the first test of this direct method to a transition-metal alloy. Our choice of NiAl is determined by our previous experience with it and by the fact that NiAl has a simple B2 structure in which each atom in the unit cell has the full cubic point group. The latter is of practical relevance since the symmetry allows us to reduce the number of independent TB parameters that need to be calculated and saved in the database. Since the TB parameters that we extract from the FPLMTO method are to be used in a recursion calculation, we do not need to extract the (two-center) Slater-Koster parameters as was done by McMahan and Klepeis.  Our motivation in performing this calculation is twofold. First, one hopes to obtain a general idea of the level of numerical accuracy that can be expected to result from a combination method like this. Second, eventually one would like to use the TB parameters extracted from this method in applications where it is too expensive computationally to use direct ab-initio methods, e.g., in studying the effects of a small concentration of impurity atoms to static and dynamic properties of alloys. Since the TB method is a non-self-consistent method, the plausibility of its result can only be judged by how well it reproduces certain “benchmark” results. In this paper the “benchmarks” are the local density of states that can also be obtained from the FPLMTO method.

Section II of this paper discusses the specific implementation of the FPLMTO method that we use for this work. The choice of basis is of the utmost importance in a TB calculation. Since the FPLMTO basis functions are well-known for being complicated we have chosen to discuss it in some details in Sec. II A. Practical points regarding the transformation to real space are listed in Sec. II B. The recursion calculation and the results that we obtained from it are discussed in Sec. III; and major points of the paper are discussed in a summary section at the end of the paper.

## II FPLMTO Method

In this paper we have used an implementation of the FPLMTO method that is developed by Wills and Price.. Several other implementations are available in the literature and, although the main formalism is largely the same, each implementation is different in its details. The FPLMTO uses a fixed basis set in calculating its hamiltonian and overlap matrix elements; it is therefore similar to the linear combination of atomic orbitals (LCAO) method. However, while in the LCAO method the basis set is determined, and fixed, at the beginning, the basis set in the FPLMTO method is determined self-consistently. The FPLMTO method is very similar to the Augmented Spherical Wave (ASW) of Williams, Kübler, and Gelatt.

One of the strengths of the FPLMTO method is that it allows unlimited number of basis functions to be assigned to each atom in the unit cell. This provides unlimited variational flexibility in finding the ground state configuration. The practical requirement of TB, however, requires us to use a minimal basis set. Furthermore, although FPLMTO method allows the use of localized ($`\kappa ^2<0`$) or unlocalized ($`\kappa ^2>0`$) basis, only the localized basis should be used for TB purposes. Williams et al. have used $`\kappa ^20.2`$ Ry exclusively in their ASW paper. They have also warned that the use of a single fixed $`\kappa `$ is rather restrictive and cited the work of Gunnarsson et al.  that shows the variation of eigenenergies with $`\kappa `$ in a molecular calculation.

### A Basis Functions

In the FPLMTO method the crystal is divided into non-overlapping muffin-tin (MT) spheres surrounding each atomic sites. The radius of the MT sphere that is centered at $`𝝉_\alpha `$, the position of the $`\alpha `$-th atom in the unit cell, is denoted by $`s_\alpha `$. On this muffin-tin geometry we construct a muffin-tin potential. This potential is spherically symmetric within each muffin tin and equal to a constant $`\mathrm{V}_0`$, the muffin-tin zero, at the interstitial region. The potential inside the MT, $`v_{\mathrm{MT}}(r)`$, does not have to join continuously to the interstitial potential $`\mathrm{V}_0`$ at the MT radius.

We emphasize that in the FPLMTO method the MT potential is used solely to construct the basis functions. The final output charge density and potential, which are obtained from a self-consistent calculation, do not necessarily have the form of a MT potential. During each iteration in the self-consistency loop, the LMTO basis functions are used to solve the Schrödinger (or Dirac) equation in a variational manner similar to the LCAO method. The resulting eigenfunctions are then used to calculate the new charge density and potential. The important difference from the LCAO method is that here, in addition to the charge density and potential, one also updates the basis functions on every step. This update is performed by constructing a new MT potential from the output potential. The constant potential $`\mathrm{V}_0`$ is given the value of the average potential at the interstitial, while the spherically-symmetric potential $`v_{\mathrm{MT}}(r)`$ is obtained from the angular average of the potential inside the muffin tin. Note that we average the potential rather than the charge density as used in an alternative scheme.

The final output MT orbitals therefore provide the “best” set of basis functions (within the variational freedom imposed by the choice of the LMTO parameters) for performing an LCAO-like calculation on the system under consideration. This view of the MT orbitals as basis functions for an LCAO-type calculation has appeared from time to time in the literature, under the name of the linear combination of muffin-tin orbitals (LCMTO) method. It has been particularly useful for calculations of the electronic structure of systems that do not possess lattice periodicity, e.g., surfaces, impurities, and atomic clusters. Unfortunately, the difficulty in performing the necessary three-dimensional integrations in this method has made it less popular relative to other methods that utilize orbitals with simpler analytical properties like the Slater- or Gaussian-type orbitals.

Each basis function in the FPLMTO method is centered around an atomic site that we will call its host site. With respect to this site, all other atomic sites will be referred to as remote sites. The basis function is constructed as a continuous patchwork of three different parts: (1) inside the host MT sphere; (2) at the interstitial region; (3) inside all remote MT spheres. Thus the basis function centered at atomic site $`𝝉_\alpha `$ inside the unit cell with lattice vector $`𝐑_i`$ has the form:

$$\varphi _{i\alpha }(𝐫)=\varphi _{i\alpha }(i\alpha |𝐫)+\varphi _{i\alpha }(I|𝐫)+\underset{i^{}\alpha ^{}}{}^{}\varphi _{i\alpha }(i^{}\alpha ^{}|𝐫).$$
(1)
Here the primed summation means that it should be carried out over all $`(i^{}\alpha ^{})(i\alpha ).`$ The MTO also carries other indices: $`L(lm)`$ which are the angular-momentum quantum numbers of the part of the MTO inside the host sphere $`\varphi _{i\alpha }(i\alpha |𝐫)`$; the tail parameter $`\kappa `$ which controls the behavior of the orbital at the interstitial region; and possibly also the principal quantum number $`n`$ which controls the number of radial nodes of $`\varphi _{i\alpha }(i\alpha |𝐫)`$. The last is usually not necessary since one normally assigns only one principal quantum number for each $`(lm)`$. In any case, for clarity we will suppress these indices unless their presence is really necessary.

An MTO is essentially an augmented spherical wave (ASW). This means that it is constructed by defining it to be a spherical wave (spherical Bessel, Neumann, or Hankel function, i.e. solution of Helmholtz equation in spherical coordinates) in the interstitial region:

$$\varphi _{i\alpha }(I|𝐫)=K_l(\kappa ,r_{i\alpha })i^lY_{lm}(\widehat{𝐫}_{i\alpha })\mathrm{\Theta }(I|𝐫),$$
(2)
with a radial part:

$$K_l(\kappa ,r)=\kappa ^{l+1}\times \{\begin{array}{cc}n_l(\kappa r)ij_l(\kappa r),\hfill & \kappa ^2<0\text{,}\hfill \\ n_l(\kappa r),\hfill & \kappa ^2>0\text{.}\hfill \end{array}$$
(3)
Here $`𝐫_{i\alpha }(𝐫𝐑_{i\alpha })`$, $`𝐑_{i\alpha }(𝐑_i+𝝉_\alpha )`$, $`r_{i\alpha }|𝐫_{i\alpha }|`$, and $`\widehat{𝐫}_{i\alpha }𝐫_{i\alpha }/r_{i\alpha }`$. The masking function $`\mathrm{\Theta }(I|𝐫)`$ is equal to 1 if $`𝐫`$ is in the interstitial region and 0 otherwise. For $`\kappa ^2<0`$, $`\kappa =i|\kappa |`$, this spherical wave is usually expressed in terms of the spherical Hankel function of the first kind:

$$K_l(\kappa ,r)=i\kappa ^{l+1}h_l^{(+)}(\kappa r)=\frac{e^z}{r}\times \{\begin{array}{cc}1,\hfill & l=0\text{,}\hfill \\ |\kappa |(1+z^1),\hfill & l=1\text{,}\hfill \\ |\kappa |^2(1+3z^1+3z^2),\hfill & l=2\text{,}\hfill \end{array}$$
(4)
where $`z|\kappa |r`$. Since this spherical wave is the envelope of the full MTO we see that, for negative $`\kappa ^2`$, the MTO has a tail that decays exponentially with distance from its host site. The decay rate is controlled directly by the magnitude of $`\kappa `$. This feature makes it suitable for tight-binding calculations.

Inside the MT spheres, the envelope spherical wave is augmented (i.e. replaced) with a solution of the Schrödinger equation for the MT potential. The augmentation process is similar to the augmentation of the plane wave in the Augmented Plane Wave (APW) method. The difference, of course, is that here we are augmenting a spherical wave instead of a plane wave. Note that the basis functions in the APW method are not suitable for tight binding since their plane-wave envelopes are not localized in space.

As we approach the host MT sphere from the outside, the basis function is determined by Eq.(2). Note that, relative to the origin at the host site, this function is separable, i.e. it is a product of a radial and an angular part. Inside the host MT sphere, this function is augmented with a solution of the Schrödinger equation for the MT potential. Since this potential is spherically symmetric inside the MT, the basis function inside the host MT is also separable:

$$\varphi _{i\alpha }(i\alpha |𝐫)=\left[A\phi _{i\alpha l}(e_{nl},r_{i\alpha })+B\dot{\phi }_{i\alpha l}(e_{nl},r_{i\alpha })\right]i^lY_{lm}(\widehat{𝐫}_{i\alpha })\mathrm{\Theta }(i\alpha |𝐫),$$
(5)
where the MT masking function $`\mathrm{\Theta }(i\alpha |𝐫)`$ is equal to 1 if $`𝐫`$ is inside the MT sphere $`(i\alpha )`$ and 0 otherwise. The energy dependence of the radial solution has been approximated with a linear combination of the radial solution $`\phi _{i\alpha l}`$ and its energy derivative $`\dot{\phi }_{i\alpha l}`$ ($`\phi _{i\alpha l}/\epsilon `$) calculated at a fixed chosen energy $`e_{nl}`$. The coefficients are obtained by requiring the function to match continuously and smoothly with the interstitial function at the MT radius:

$$\left[\begin{array}{cc}\phi _{i\alpha l}(s_{i\alpha })& \dot{\phi }_{i\alpha l}(s_{i\alpha })\\ \phi _{i\alpha l}^{}(s_{i\alpha })& \dot{\phi }_{i\alpha l}^{}(s_{i\alpha })\end{array}\right]\left[\begin{array}{c}A\\ B\end{array}\right]=\left[\begin{array}{c}K_l(\kappa ,s_{i\alpha })\\ K_l^{}(\kappa ,s_{i\alpha })\end{array}\right].$$
(6)
Here a prime denotes differentiation with respect to the radial coordinate: $`\phi ^{}\phi /r`$.

The augmentation process at the remote spheres proceeds rather differently from the one for the host sphere. As we approach a remote MT sphere from the interstitial region, the basis function is again determined by Eq.(2). The important difference from the host-sphere case is that relative to the origin at the remote site the interstitial function is not separable. We can, however, express the spherical Hankel function at the interstitial region as a sum of separable functions in the neighbourhood of the remote site:

$$K_l(\kappa ,r_{i\alpha })i^lY_L(\widehat{𝐫}_{i\alpha })=\underset{L^{}}{}J_l^{}(\kappa ,r_{i^{}\alpha ^{}})i^l^{}Y_L^{}(\widehat{𝐫}_{i^{}\alpha ^{}})S_{L^{}L}(\kappa ,𝐑_{i^{}\alpha ^{}}𝐑_{i\alpha }),r_{i^{}\alpha ^{}}<|𝐑_{i^{}\alpha ^{}}𝐑_{i\alpha }|,$$
(7)
with the spherical Bessel function as the radial part in each term:

$$J_l(\kappa ,r)=\kappa ^lj_l(\kappa r)=\frac{1}{z}\times \{\begin{array}{cc}\mathrm{sinh}z,\hfill & l=0\text{,}\hfill \\ |\kappa |^1(\mathrm{cosh}zz^1\mathrm{sinh}z),\hfill & l=1\text{,}\hfill \\ |\kappa |^2[(3z^2+1)\mathrm{sinh}z3z^1\mathrm{cosh}z],\hfill & l=2\text{.}\hfill \end{array}$$
(8)
The addition theorem for the Bessel functions, Eq.(7), is a consequence of the fact that both sides of the identity are solution of the translationally-invariant Helmholtz equation. The spherical Hankel function, $`K_l(\kappa ,r)`$, is regular everywhere except at the origin, while the spherical Bessel function, $`J_l(\kappa ,r)`$, is regular everywhere except at infinity. The domain of validity of Eq.(7) is just the domain where both sides of the identity are regular.

The radial part of each term in the one-center expansion, Eq.(7), is centered at the remote site $`𝐑_{i^{}\alpha ^{}}`$. Thus the $`(\phi ,\dot{\phi })`$ augmentation of the interstitial function, Eq.(2), at this remote site can be performed by augmenting each $`J_l(\kappa ,r_{i^{}\alpha ^{}})`$ with the radial part of the solution of the Schrödinger equation for the spherically-symmetric MT potential centered at $`𝐑_{i^{}\alpha ^{}}`$. Thus we replace:

$$J_l(\kappa ,r_{i^{}\alpha ^{}})\left[C_l\phi _{i^{}\alpha ^{}l}(e_{nl},r_{i^{}\alpha ^{}})+D_l\dot{\phi }_{i^{}\alpha ^{}l}(e_{nl},r_{i^{}\alpha ^{}})\right].$$
(9)
The coefficients are again obtained by requiring continuity in the value of the function and its first derivative with those of the interstitial envelope function at the MT radius:

$$\left[\begin{array}{cc}\phi _{i^{}\alpha ^{}l}(s_{i^{}\alpha ^{}})& \dot{\phi }_{i^{}\alpha ^{}l}(s_{i^{}\alpha ^{}})\\ \phi _{i^{}\alpha ^{}l}^{}(s_{i^{}\alpha ^{}})& \dot{\phi }_{i^{}\alpha ^{}l}^{}(s_{i^{}\alpha ^{}})\end{array}\right]\left[\begin{array}{c}C_l\\ D_l\end{array}\right]=\left[\begin{array}{c}J_l(\kappa ,s_{i^{}\alpha ^{}})\\ J_l^{}(\kappa ,s_{i^{}\alpha ^{}})\end{array}\right].$$
(10)
The part of basis function inside a remote MT sphere is thus given by:

$$\varphi _{i\alpha }(i^{}\alpha ^{}|𝐫)=\mathrm{\Theta }(i^{}\alpha ^{}|𝐫)\underset{L^{}}{\overset{ll_m}{}}\left[C_l^{}\phi _{i^{}\alpha ^{}l^{}}(r_{i^{}\alpha ^{}})+D_l^{}\dot{\phi }_{i^{}\alpha ^{}l^{}}(r_{i^{}\alpha ^{}})\right]i^l^{}Y_L^{}(\widehat{𝐫}_{i^{}\alpha ^{}})S_{L^{}L}(\kappa ,𝐑_{i^{}\alpha ^{}}𝐑_{i\alpha }).$$
(11)
The summation on the r.h.s. of the equation ideally runs to infinity. Practically, a converged total energy is achieved using $`ll_m68`$ in most cases. Note that, inside the remote sphere, an angular-momentum expansion is necessary since the resulting function needs to match the envelope interstitial function at its MT sphere. Relative to this remote site, the envelope function contains multiple harmonics since it is centered at a different site.

The MT orbital, constructed from Eqs.(2), (5), and (11), is continuous and smooth everywhere. Note, however, that it is not normalized and it does not have a separable form:

$$\varphi _{i\alpha }(𝐫)u(r_{i\alpha })Y_{lm}(\widehat{𝐫}_{i\alpha }),$$
(12)
in other words, the MTO does not have good angular-momentum quantum numbers. Although the parts of the MTO inside the host MT and at the interstitials are separable, the presence of the remote sites breaks the perfect spherical symmetry around the host site, and therefore also removes the separability. The MTO, however, does obey exactly the point symmetry of its host site. Moreover, since the non-separable contributions come from the remote sites, where the amplitude of the exponentially-decaying envelope function is small, to a good approximation the MTO does transform as if it has the quantum numbers $`(lm)`$ of the part of the orbital inside the host sphere. The fact that the MTO is not fully separable should be noted especially if one intends to map the FPLMTO result into a TB form using the Slater-Koster parametrization method. In the SK-LCAO parametrization method the basis functions are assumed to be atomic orbitals; these orbitals are separable since they are derived from a single-atom spherically-symmetric potential.

### B Fourier Transform

In applications of the FPLMTO method to ideal crystalline systems, instead of using the basis functions in real space, one works with the Bloch functions:

$$\psi _\alpha (𝐤,𝐫)=\frac{1}{\sqrt{N}}\underset{i}{}\mathrm{exp}(i𝐤𝐑_i)\varphi _{i\alpha }(𝐫).$$
(13)
The energy bands are obtained by solving the eigenvalue equation:

$$det|H_{\alpha \beta }(𝐤)ES_{\alpha \beta }(𝐤)|=0.$$
(14)
The rank of the matrices is equal to the total number of the orbitals used for all the atoms within the unit cell. The hamiltonian matrix is given by integration over an arbitrarily chosen unit cell $`𝐑_n`$:

$$H_{\alpha \beta }(𝐤)=\psi _\alpha (𝐤,𝐫)|H|\psi _\beta (𝐤,𝐫)=\frac{1}{N}\underset{ij}{}e^{i𝐤(𝐑_j𝐑_i)}\underset{n}{}\left[\varphi _{i\alpha }(I_n)|H|\varphi _{j\beta }(I_n)+\underset{\gamma }{}\varphi _{i\alpha }(n\gamma )|H|\varphi _{j\beta }(n\gamma )\right],$$
(15)
and similar expression is obtained for the overlap matrix by replacing the hamiltonian operator $`H`$ with identity. The vectors have been defined by:

$$𝐫|\varphi _{i\alpha }(n\gamma )=\varphi _{i\alpha }(n\gamma |𝐫),$$
(16)
which is the part of the basis function, centered at the host site $`i\alpha `$, inside the MT sphere $`n\gamma `$ or the interstitial region $`I_n`$, which is defined to be the interstitial region within the unit cell $`𝐑_n`$. Cross terms are absent from the square bracket in Eq.(15) since the overlap of their masking functions is zero.

It should be mentioned that in the actual FPLMTO formalism, the matrix elements are not computed using a multicenter expansion (over all lattice sites) as shown in Eq.(13). Rather, using the addition theorem for the tail of the MT orbital, Eq.(7), the multicenter expansion is first transformed into a one-center expansion (over angular momenta at a single site), and this sum is the one that is actually computed.  This, however, is merely a computational contrivance, albeit a pivotal one in the FPLMTO method, and should not mask the true structure of the basis function in k-space as a simple Bloch sum of the MT orbitals.

In the FPLMTO method of Wills et al., the fact that each term of the matrix elements is to be calculated only over a certain region, either inside the MT sphere or the interstitial region as seen in Eq.(15), is used to facilitate its computation. Instead of using the true MTO, which is complicated, one uses a pseudo basis orbital. This orbital is equal to the true MTO inside the relevant integration region but outside the region it is replaced with a smooth function which is chosen to facilitate the computation. The integration can then be performed efficiently by working in the momentum space.

Using the translational properties of the Bloch functions, the matrix element can be written as:

$$H_{\alpha \beta }(𝐤)=\underset{j}{}\mathrm{exp}(i𝐤𝐑_j)H_{\alpha \beta }(𝐑_j),$$
(17)
with the matrix elements in real space:

$$H_{\alpha \beta }(𝐑_j)=\underset{n}{}\left[\varphi _{0\alpha }(I_n)|H|\varphi _{j\beta }(I_n)+\underset{\gamma }{}\varphi _{0\alpha }(n\gamma )|H|\varphi _{j\beta }(n\gamma )\right].$$
(18)
To calculate the inverse transform of Eq.(17) we use the fact that the real-space matrix element $`H_{\alpha \beta }(𝐑)`$ decays exponentially with the distance $`|𝐑|`$. It is thus reasonable to ignore matrix elements between orbitals separated by a distance greater than a certain cutoff radius $`R_c`$. In practice this can be performed by working with a finite crystal which contains $`N_x`$ unit cells along the $`x`$-direction such that $`N_xa_x>2R_c`$, with $`a_x`$ being the lattice constant along the $`x`$-direction, and similar ranges for the $`y`$ and $`z`$ directions. The real-space matrix elements can then be obtained from a discrete Fourier transform:

$$H_{\alpha \beta }(𝐑_j)=\frac{1}{N}\underset{𝐤}{}\mathrm{exp}(i𝐤𝐑_j)H_{\alpha \beta }(𝐤).$$
(19)
Here $`N=N_xN_yN_z`$ and $`k_i=\pi n_i/(N_ia_i)`$ with $`N_i<n_iN_i`$ for $`i\{x,y,z\}`$.

Some practical remarks should be mentioned regarding this calculation of the inverse Fourier transform. First, the output matrix elements from the FPLMTO are calculated using the spherical harmonic $`i^lY_{lm}(\widehat{𝐫})`$ as the angular part of the basis function inside its host sphere. Usually it is easier to work with real basis functions and this can be achieved by making linear combinations of the spherical harmonics. Second, the basis functions are not normalized. The normalized matrix elements are obtained by using the overlap matrix elements:

$$H_{\alpha \beta }^{}(𝐑)=S_{\alpha \alpha }(\mathrm{𝟎})^{\frac{1}{2}}H_{\alpha \beta }(𝐑)S_{\beta \beta }(\mathrm{𝟎})^{\frac{1}{2}}.$$
(20)
Third, the actual displacement vector connecting the centers of the two orbitals in $`H_{\alpha \beta }(𝐑)`$ is not $`𝐑`$ but is instead $`(𝐑+𝝉_\beta 𝝉_\alpha )`$. Instead of the original hamiltonian matrix in Eq.(17), it is preferable to work with a modified matrix that is related to the original by a unitary transformation:

$$H_{\alpha \beta }^{\prime \prime }(𝐤)=e^{i𝐤(𝝉_\beta 𝝉_\alpha )}H_{\alpha \beta }(𝐤)=\underset{j}{}e^{i𝐤(𝐑_j+𝝉_\beta 𝝉_\alpha )}H_{\alpha \beta }(𝐑_j).$$
(21)
This positions the origin at the center of orbital $`(0\alpha )`$ and allows the matrix elements to obey the point symmetry of the corresponding site. It may therefore be used to reduce the number of matrix elements that need to be stored. Note that the eigenvalues are unaffected by this unitary transformation.

### C Löwdin Orthogonalization

The real-space hamiltonian matrix elements in Eq.(19) (together with the overlap matrix elements $`S_{\alpha \beta }`$ which can be obtained similarly) can be used directly in a nonorthogonal recursion calculation. Unfortunately, this requires us to obtain the inverse of the overlap matrix in real space, which is a non-trivial computational task. In this paper we instead use an orthogonal hamiltonian which is obtained by Löwdin orthogonalization of the original hamiltonian and overlap matrices.  An alternative scheme is to use the chemical pseudopotential approach and work with a non-hermitian matrix $`S^1H`$. The advantage of Löwdin orthogonalization is that it is a symmetry transformation: angular symmetry is preserved, i.e., orbitals with angular momentum $`l`$ are not mixed by the transformation with orbitals with $`l^{}l`$.  This means that the real-space matrix elements can still be parametrized by Slater-Koster parametrization, e.g., using the procedure proposed by McMahan and Klepeis.

Since we are working with a perfect crystal, the Löwdin orthogonalization process can be applied to the momentum-space matrix elements. One first diagonalize the overlap matrix:

$$D=U^{}SU,$$
(22)
where $`U`$ is the matrix containing the eigenvectors of $`S`$ as its columns, and $`D`$ is a diagonal matrix containing the eigenvalues of $`S`$. Using these matrices, we define a new matrix:

$$A=UD^{\frac{1}{2}}U^{}.$$
(23)
Since the overlap matrix is a positive-definite hermitian matrix, its eigenvalues are all positive real numbers and the matrix $`D^{\frac{1}{2}}`$ is well defined. In practice, however, small negative eigenvalues of $`S`$ may sometimes appear due to the finite machine precision used in the computation. This problem is especially relevant if the diagonalization of the overlap matrix is performed in real space and when one uses an overlap matrix which is necessarily approximate due to the finite cutoff imposed on the range of the matrix elements or due to the Slater-Koster two-center approximation used in parametrizing them. One commonly-used fix to this problem is to add a diagonal matrix with small elements to make the overlap matrix positive definite. The orthogonal Löwdin hamiltonian matrix is obtained by sandwiching the original hamiltonian matrix with $`A`$ matrices:

$$H^{(L)}=AHA.$$
(24)
The nonorthogonal eigenvalue problem, Eq.(14), is then transformed to an equivalent orthogonal one:

$$det|H_{\alpha \beta }^{(L)}(𝐤)E\delta _{\alpha \beta }(𝐤)|=0.$$
(25)
The orthogonal hopping parameters are obtained by substituting $`H_{\alpha \beta }^{(L)}(𝐤)`$ in Eq.(19). These are the parameters that we use in the recursion calculation which is described in the rest of this paper.

## III Recursion Calculation

The recursion method, initially introduced to electronic-structure calculation by Haydock, is one of the most efficient and stable ways of calculating the electronic Green’s function from the TB hamiltonian describing the system:

$$G_{uv}(E)=u|(EH)^1|v,$$
(26)
where $`|u`$ and $`|v`$ are two arbitrary vectors in the Hilbert space for the problem. The local density of states (LDOS) is a projection of the total density of states to a local orbital:

$$\rho _\alpha (E)=\underset{n}{}|\alpha |\mathrm{\Psi }_n|^2\delta (EE_n),$$
(27)
where $`|\mathrm{\Psi }_n`$ is an eigenvector with eigenvalue $`E_n`$. This can be obtained from the Green’s function:

$$\rho _\alpha (E)=\frac{1}{\pi }\mathrm{Im}\left[G_{\alpha \alpha }(E+iϵ)\right],$$
(28)
with $`ϵ0^+`$. The token Im$`[\mathrm{}]`$ denotes the operation of taking the imaginary part of its argument. The recursion method outputs the Green’s function as a continued fraction:

$$G(E)=\frac{1}{(Ea_0)\frac{b_1^2}{\left(Ea_1\right){\scriptscriptstyle \frac{b_2^2}{\mathrm{}}}}}.$$
(29)
Note that there are two conventions used in the literature regarding the indexing of the $`a`$-parameters in the above continued fraction. One convention is to assign $`a_0`$ as the first $`a`$-parameter, while the other use $`a_1`$ as the first parameter. Here we follow the first convention. The calculation of the LDOS is therefore reduced to the computation of the strings of recursion coefficients $`\{a_n\}`$ and $`\{b_n\}`$. The procedure for performing this has been described extensively in the literature and will not be repeated here.

We have used the recursion method to calculate the LDOS for NiAl crystal. This alloy crystallizes in B2 structure with two atoms per (cubic) unit cell. One component (Ni or Al) sits at $`(0,0,0)`$, while the other is at the body-center of the cube, $`(\frac{1}{2},\frac{1}{2},\frac{1}{2})`$. We obtained the real-space orthogonal hopping parameters by direct extraction and Löwdin transformation as described in the preceding section. The parameters were extracted from a single-$`\kappa `$ FPLMTO calculation using 9 $`spd`$ orbitals and $`\kappa ^2=0.04`$ Ry for each Ni or Al atom. This rather small value of $`\kappa `$ produces basis functions with relatively long range (the envelope of the basis function decays roughly as $`e^{|\kappa |r}`$). Using a larger negative value of $`\kappa `$, however, produces a greater discrepancy between the output bands from this single-$`\kappa `$ calculation and the corresponding accurate bands from a multiple-$`\kappa `$ calculation. The chosen value of $`\kappa `$ that we used is obtained by compromising the need to have a short-ranged TB basis functions with the requirement to have a precision comparable to accurate multiple-$`\kappa `$ result.

All of the extracted parameters are saved in a database for use in the recursion calculation (crystal symmetry was used to reduce the size of the database). They were calculated for separation distance of up to 8 lattice constants. The obtained TB parameters fall roughly one order of magnitude for each increase in separation distance of one lattice constant. The recursion calculation was performed using a cubic cluster of $`(16)^3`$ unit cells and the LDOS’s were computed for Ni and Al orbitals within the unit cell at the center of the cluster. Since the TB parameters connecting the orbitals at the central unit cell to the ones at the surface of the cluster are practically zero, the calculated LDOS should not depend on the size of the cluster if we increase it further.

Figs. 2 and 3 show the calculated recursion parameters for Ni-$`d`$ and Al-$`p`$ states, respectively. Here the $`d`$-state is defined by:

$$|d=\frac{1}{\sqrt{5}}\underset{m=2}{\overset{2}{}}|l=2,m,$$
(30)
with a corresponding definition for the $`p`$-state. It is seen that the parameters tend to settle around certain constant values. Indeed it is known that as $`n\mathrm{}`$:

$$a_n(E_t+E_b)/2,b_nW/4,$$
(31)
where $`E_t`$ and $`E_b`$ are the values of energy at the top and bottom of the spectrum, respectively, and $`W=(E_tE_b)`$ is the total bandwidth. Furthermore, asymptotically the parameters should oscillate around their limit values with a decaying amplitude as $`n\mathrm{}`$. The frequency and the rate of decay are determined by specific features (the Van Hove singularities) of the spectrum.

The recursion parameters shown in Figs. 2 and 3 are seen to roughly follow this decayed oscillation prediction. However, for large values of $`n`$ (greater than $`n_m40`$), the amplitudes of the deviation tend to increase rather than decrease. We interpret this as the limit at which the machine precision noise start to interfere with the calculation. This interpretation is supported by our experience that incorporating the recursion parameters with $`n>n_m`$ in the reconstruction of LDOS in general does not help improve the agreement, and indeed it tends to reduce the agreement, with the accurate LDOS obtained from the FPLMTO method. Note that the main workhorse in the recursion method involves a matrix-vector multiplication with the dimension of the vector being the total number of orbitals in the cluster. In our case this dimension is equal to 73728 (= $`16^3`$ unit cells $`\times `$ 2 atoms/cell $`\times `$ 9 orbitals/atom). Compounded with this large dimension is the fact that the elements of the matrix and the vector in general have widely different orders of magnitude which makes the calculation rather prone to rounding errors.

Figs. 4 display the main results in this paper: the s,p,d LDOS of Ni and Al in NiAl as calculated using the TB recursion method. These LDOS’s are compared with the corresponding accurate spectra obtained from Brillouin zone integration using atom and angular-momentum projection of the eigenvalues obtained directly from the FPLMTO method. The agreement in general is good although the TB recursion method is unable to precisely reproduce the sharp peaks in the FPLMTO spectra. Not very surprisingly, the best fit is obtained for the tightly-bound Ni-$`d`$ state.

The limited number ($`n_m40`$) of reliable recursion parameters that we can use here is generally not sufficient for acceptable reproduction of the FPLMTO LDOS. Furthermore, construction of LDOS by directly taking the imaginary part of the Green’s function, Eq.(29), normally introduces an undesirable feature in the output LDOS in the form of spurious rapid oscillations near the edges of the spectrum, as can still be seen in, e.g., Fig. 4a. The simplest extrapolation method for the recursion parameters is the square-root terminator method: one simply replaces all $`a_n`$ and $`b_n`$ for $`n>n_m`$ with constants $`a_{\mathrm{}}`$, and $`b_{\mathrm{}}`$. The tail of the continued fraction for the Green’s function can then be summed exactly as a square-root expression. The best values for the terminating coefficients $`a_{\mathrm{}}`$ and $`b_{\mathrm{}}`$ can be obtained using the Beer-Pettifor method.. Other ways to construct a terminator for the continued fraction have been suggested, but for our calculation we have chosen to use the Beer-Pettifor method. Finally, further smoothing of the spectra shown in Fig. 4 have been performed using a Chebyshev-polynomial method.

## IV Summary

In this paper, we have described a method for directly extracting real-space TB parameters from a FPLMTO method. The basis functions used in the FPLMTO calculation and their construction are described in considerable details. Special emphasis has been placed on the fact that these basis functions do not have a well-defined angular-momentum quantum number, although to a good approximation they do. We believe this fact should be kept in mind in transcribing FPLMTO matrix elements to the Slater-Koster two-center form, since there it is implicitly assumed that the basis functions do have a good angular-momentum quantum number.

This direct extraction method has been applied to an intermetallic alloy NiAl. The TB parameters extracted have been used as input to a real-space calculation of local densities of states using the recursion method. We believe this is the first application of the direct extraction method to a real-space calculation of the electronic structure of an intermetallic alloy. The good agreement between the LDOS obtained using the TB extracted from the FPLMTO method and the accurate LDOS obtained directly from the FPLMTO method shows the feasibility of using this method to study more complicated cases, e.g., alloys with small concentration of impurity atoms.

This work was supported by AF-OSR Grant No. F49620-99-1-0274. DD would like to thank Dr. M.J. Mehl of the Naval Research Laboratory, Washington D.C., for useful discussions and insightful inputs regarding the NRL Tight-Binding code.