text:
# VLT/ISAAC H-band spectroscopy of embedded massive YSOs

## 1 Introduction

Little is known about the earliest stages of the formation of massive stars. It is unclear whether the general picture applicable to the formation of low-mass stars, i.e. an accretion disk and infall onto the central star (Shu et al. (1987)) can be scaled up to masses higher than $``$10M☉ (in the case of spherically symmetric accretion), as radiation pressure from the high mass stellar core on dust mixed with the infalling gas prevents further accretion. Alternative explanations involving the coalescence of intermediate mass protostars in the dense central regions of rich clusters have been proposed (Bonnell et al. (1998)).

Observationally, the earliest stages of massive star formation are only accessible by infrared and radio techniques, because of the presence of $``$10 – 100 magnitudes of visual extinction. Ultracompact (UC) H ii regions represent this earliest phase in the evolution of a very massive star. UC H ii regions are typically identified using radio surveys (e.g. Bronfman et al. (1996); Walsh et al. (1999)). Massive YSOs capable of ionizing a UC H ii region are recognisable by their extremely red IRAS colours (Wood & Churchwell (1989); Osterloh et al. (1997)), and the IRAS database has become a rich source to be mined for intermediate and high-mass YSOs using colour selection criteria (e.g. Persson & Campbell (1987); Wood & Churchwell (1989); Campbell et al. (1989); Chan et al. (1996)).

Despite this, measurements of the fundamental parameters of UC H ii regions, and more importantly those of the massive embedded stars which ionize the region, cannot easily be derived from the radio and far-IR observations. This is due principally to uncertainties in distances and in the indirect methods used to convert radio continuum or IRAS fluxes to the physical properties of the embedded, ionizing OB star (Walsh et al. 1997; Kaper et al. (2002) and references therein).

To circumvent these difficulties, observational efforts are now beginning to probe directly the photospheres of massive embedded YSOs using observations in the near-infrared. Such wavelengths are accessible from the ground (1 – 5 $`\mu `$m, JHKLM bands) and the extinction owing to the gas and dust in which the YSO is embedded is much less than in the optical. Moreover, at thermal infrared wavelengths, emission from the dust becomes too great, and the flux from the photosphere itself can no longer be detected. The first spectral classification of the ionizing star of a UC H ii region was made by Watson & Hanson (1997) in the K-band.

While an important aim of the present paper is to clarify the status of the three objects discussed here as massive embedded YSOs, a second no less important aim is to show that intermediate resolution spectroscopy in the H-band (where veiling from circumstellar emission is minimal) is a necessary and valuable tool to detect, classify and study the photospheres of such objects. This paper attempts to build on recent work on the classification of OB stars in the H-band, and to apply it to highly obscured objects.

### 1.1 Classification spectroscopy of OB stars in the H-band

Recent developments in classification spectroscopy of OB stars in the H- and K-bands suggest application to objects too deeply obscured to be optically visible. In the earliest work, both Lancon et al. (1992) and Dallier et al. (1996) presented libraries of stellar spectra in the H-band, but these were largely directed to stellar population analyses. More recently, the most important study directly relevant here has been that of Hanson et al. (1998), who have identified a number of H-band temperature and luminosity diagnostics for OB stars in R=2000 spectra, involving the ratio of Br 11 (1.681 $`\mu `$m) to He i (1.700 $`\mu `$m) and He ii (1.693 $`\mu `$m). Using a full grid of spectra representing both dwarf and supergiant luminosities, these authors find that the presence of all three lines definitively indicates a spectral type earlier than late O. With combinations of only two lines present, in various strengths, dwarfs can be divided into O9 – B1, B2 – B5, and late B – early A type.

This work overlaps with that of Meyer et al. (1998), who have observed 85 MK standards of near-solar metallicity, at R$``$3000 in the H-band. They find a number of strong Teff and luminosity-sensitive lines in the H-band spectra of stars of spectral type A, and later, all the way down to M5. These two studies show that the classification of heavily embedded YSOs, and the near-IR study of their photospheres, could be carried out in the H-band for all spectral types from O to M, using intermediate resolution, high signal-to-noise spectra. This could be advantageous, since the veiling of photospheric absorption lines by continuum emission from hot circumstellar dust is much less severe than in the K-band (Meyer et al. 1998; Ishii et al. (2001)). Taken together, these studies suggest that while K-band spectroscopy has been extremely successful in studies of low-mass YSOs (see Greene & Lada (1996) for one example among many) classification in the range B0 – G0 is best performed in the H-band, as pointed out by Luhman & Rieke (1998).

## 2 Observations and data reduction

### 2.1 Target selection

Targets for this study were selected from recent literature on UC H ii regions and IRAS-selected catalogues of massive YSOs (see Table 1). The main criteria were that such targets would be sufficiently bright in the H-band to be accessible for medium resolution ISAAC spectroscopy in short exposure times, and that they should be present in the publicly available 2MASS dataset, to avoid the need for prior imaging. A further important criterion was that the 2MASS point sources should have no optical counterparts in the GSC2.2 and USNO catalogues. The three selected targets are listed in Table 1 together with 2MASS magnitudes and colours, Av, and indirect spectral type estimates from the given references. The mid-infrared fluxes are shown in the IRAS colour-colour diagram (Fig. 1), where different locations occupied by sources of various types are indicated. IRAS 17175-3544 is too faint in the J-band for 2MASS to give an accurate magnitude in this band. However, this source shows a strong near infrared excess as shown in Fig. 2 where its JHK colours were derived instead from the deeper imaging of Tapia et al. (1996).

Fig. 3 shows the 2MASS H-band images, covering 152″ square (the FOV of ISAAC) for each target. Fig. 3 can be usefully compared to Fig. 19 of Walsh et al. (1999) for IRAS 17175-3544, who show a K-band image overlaid with the positions of methanol and OH maser sites and 8.64 GHz radio continuum emission, which are shown to be coincident with the near-infrared source. Furthermore, L-band imaging of IRAS 17175-3544 by the same authors indicates that the source is extremely bright and red (L=3.8, K-L=4.7). In the cases of IRAS 17441-2910 and IRAS 18079-1756, 2MASS H-band data indicate that the bright sources shown in Fig. 3 lie within 2–3″ of the IRAS positions, allowing unambiguous identification of the targets.

The presence of nebulosity associated with all targets is clearly shown in the H-band images and even more obvious in the K-band (see Fig. 4, which uses the same intensity scales as the equivalent H-band images). This strengthens the identification of the sources in the 2MASS data with the IRAS sources given in the literature. The very bright source IRAS 17441-2910 appears to be multiple, but manipulation of the image indicates that there is one single, major source for which the point spread function in the 2MASS H-band image is symmetric. Therefore it was included as a target under the assumption that this one source would contribute by far the greater flux to the H-band spectrum.

### 2.2 Data extraction

Observations were carried out at VLT/UT1 (Antu) with the ISAAC spectrograph, in SW mode at medium resolution (R$``$5000), in service mode, on 2002 May 9 and 10. The total integration time for each object is given in Table 1 and totals 150 min. A central wavelength of 1.71 $`\mu `$m was selected to provide a wavelength range of $``$1.67 – 1.75 $`\mu `$m, covering lines diagnostic of an OB spectral type, most importantly H i Br 10 1.7362 $`\mu `$m, Br 11 1.6811 $`\mu `$m, He i 1.7007 $`\mu `$m and He ii 1.693 $`\mu `$m. Data reduction was performed with the ISAAC pipeline package ECLIPSE (version 4.1.2). Firstly, all data and calibration frames were ”de-ghosted” using the ghost routine. Next, master flat fields were created for each night using the sw\_spflat routine. Co-added spectrum images were then created using the sw\_spjitter routine. Full details of the ISAAC reduction pipeline can be found in Amico et al. (2002).

Spectra were then extracted using the standard apall routine in iraf111iraf is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc. (AURA) under cooperative agreement with the National Science Foundation., and wavelength calibrated using the sky OH emission lines found throughout the spectra. In all cases $``$17 OH lines were used, yielding RMS values on the fits to the dispersion of $`<`$0.1 Å.

Telluric divisor stars were selected at the time the observations were made, and were typically observed to within 0.05 airmasses of the science objects. Telluric standards were reduced, extracted and wavelength calibrated in exactly the same way as the science data. As the divisor stars are all of spectral type B, Br 11, Br 10 and He i profiles were estimated by eye and directly fitted out of the standard spectra, and pure telluric spectra created in these regions by dividing the spectrum by the fit. These sections were then spliced into the rest of the standard star spectra (assumed free of photospheric lines) to create normalised pure telluric spectra, one derived from each standard star observation.

One problem is that in all standard star spectra He i is itself blended with a telluric line, too closely for the components to be separated. Therefore the whole He i \+ telluric profile had to be fit, and replaced by a flat ‘continuum’ region. Hence the region around He i remains uncorrected for telluric absorption in all target spectra shown in Fig. 5. The problem is dealt with further in Sect. 4, where we have used synthetic spectra to model the telluric standards (see Fig. 6).

Object spectra were divided by those of the appropriate normalised standard (after fitting out photospheric lines) using the iraf package telluric. This method retains the same flux scale (number of counts) and spectral slope as the uncorrected spectra. The signal-to-noise ratio (S/N) of the final spectra is $``$100 and in the case of the very bright target IRAS 17441-2910, it approaches 150.

## 3 Background for individual objects

The three target objects are located in heavily obscured regions in the Galactic plane. Strong extinction in the respective areas is especially obvious from 2MASS K-band images (Fig. 4). This results in a sharp drop in 2MASS K-band detections within a few arcminutes of the target objects. The obscuring material appears to have a filamentary structure, especially visible near IRAS 17441-2910 and IRAS 18079-1756. Neither of these two objects are located within a known massive star forming region or giant molecular cloud.

By contrast, IRAS 17175-3544 lies in the North-Eastern part of the extended molecular cloud complex NGC6334 (photometric distance 1.76 kpc, Neckel 1978). This star forming region is known as NGC 6334 I (Emerson et al. 1973; McBreen et al. 1979) and contains three IR sources IRS-I 1, IRS-I 2 and IRS-I 4 (Harvey & Gatley 1983). IRS-I 1 and IRS-I 2 are within a few arcsec of each other, and both have been proposed to be the driving source of a high-velocity ($`70\mathrm{km}\mathrm{s}^1`$) CO bipolar outflow (see Persi et al. 1996).

IRAS 17175-3544 itself appears to be associated with IRS-I 1 (see Fig. 19 of Walsh et al. (1999) and Persi et al. 1996). High resolution near-IR imaging showed IRS-I 1 to be complex, consisting of at least 4 components (Persi et al. 1996). The nearby (6″) source IRS-I 2 is a mid-IR source (detected at 20 and 30 $`\mu `$m). Whether or not this source contributes to the H-band spectrum is determined by its extension and position with respect to the slit, but no source brighter than K $``$ 17.5 is detected at the position of IRS-I 2 by Persi et al. (1996), so we consider such a contribution unlikely. This source could however contribute to the IRAS fluxes.

Further evidence for the massive YSO nature of at least two of the objects discussed in this work, namely IRAS 17175-3544 and IRAS 18079-1756, is that both have IRAS low resolution (8 – 23 $`\mu `$m) spectra classified as type H by Kwok et al. (1997). These authors state that the vast majority of IRAS LRS spectra of this type, exhibiting very red continua, are associated with H ii regions.

Another noteworthy feature concerning IRAS 17175-3544 is that it is surrounded by an embedded stellar cluster (centred on IRS-I 2), which was detected in the K-band by Straw et al. (1989). The cluster is clearly visible in Fig. 4. Deeper imaging by Tapia et al. (1996) yielded a stellar density estimate of about $`1200\mathrm{pc}^3`$, and cluster size of 0.6 pc (for $`\mathrm{D}=1.74\mathrm{kpc}`$).

In addition to the cluster near IRAS 17175-3544, our acquisition images revealed that IRAS 18079-1756 is also accompanied by a embedded cluster, which the 2MASS image in Fig. 3 is too shallow to fully reveal. It is not clear whether there is an increase in stellar density near the third source IRAS 17441-2910.

Finally, apart from the location in an obscured region and the presence of a cluster (near at least two sources), a third common characteristic of the target sources is that all are associated with maser emission and/or are known to be outflow sources, as has already been mentioned for IRAS 17175-3544. IRAS 17441-2910 shows evidence for OH masering (Caswell 1998), while IRAS 18079-1756 displays $`\mathrm{H}_2\mathrm{O}`$ maser emission located at $``$12″ distance and is known to be a CO outflow source (Osterloh et al. 1997). To conclude, all three sources show all the characteristics of young, actively forming massive stars.

## 4 Results from individual spectra

### 4.1 IRAS 17175-3544

Fig. 5 shows the observed H-band spectrum of IRAS 17175-3544 (upper large panel). As has already been discussed, Walsh et al. (1999) find maser and radio continuum sources exactly coincident with the near-IR source discussed here. IRAS 17175-3544 is one of the few sources for which Walsh et al. (1999) can unambiguously derive a spectral type. As a result of a methanol maser emission survey to identify UC H ii regions from an IRAS-selected sample, Walsh et al. (1997) calculated total luminosities and spectral types directly from the IRAS fluxes and the kinematic distances given by the maser emission. However, the large size of the IRAS beam ($``$1′), may encompass a number of emitting regions, so the IRAS-based spectral types cannot be definitive.

In the case of IRAS 17175-3544, though, extrapolation of the near-IR SED to the far-IR shows that the near-IR counterpart will be the one dominant star in the field. Using a kinematic distance of 2.2 kpc, Walsh et al. (1999) find O6 – O7 (from a total luminosity of 12.1$`\times `$104 L☉).

The H-band spectrum does not slope as steeply to the red as the other two objects in the sample, perhaps suggesting it is less embedded; yet the infrared excess is by far the largest of this sample. Narrow, relatively weak emission from H i Br 10 and Br 11 is clearly seen. In Fig. 5, the region around the He i line at 1.7007 $`\mu `$m is not corrected for telluric lines, for reasons already discussed. However, there is certainly a suggestion of He i emission. This is clarified in Fig. 7, where in order to perform telluric correction in this region, we have fitted the complete telluric standard spectrum with a synthetic spectrum. For the standard for IRAS 17175-3544, we adopt Teff = 20000 K, log $`g`$ = 3.5 and vsin$`i`$ = 200 km s-1 (see also Fig. 6). The He i line for these parameters is shown in the left panel of Fig. 7 (dashed line) together with the blended He i \+ telluric profile in the standard star. After division, a pure telluric line is obtained (Fig. 7, right panel, dotted line). Using this profile, the He i region of IRAS 17175-3544 was corrected for telluric absorption, yielding a clear emission line rising to $``$ 20% above the continuum level. Using this line, together with the narrow H i emission, we derive an LSR radial velocity of +15.9 $`\pm `$ 1.5 km s-1 for this object. This is somewhat different from the methanol maser emission velocity given by Walsh et al. (1997) (–10 km s-1) and could indicate that the atomic and molecular emission do not arise in precisely the same circumstellar region. However we note that the velocity resolution of our data is insufficient to press this point. It is clear, however, from inspection of the 2D spectrum images for this object, that the H and He emission lines are extended spatially, and hence are of nebular origin.

By themselves, emission lines of H and He suggest an early-type object, but do not allow derivation of a spectral type. We may speculate that, as the emission is much weaker than in IRAS 17441-2910, either the exciting central object is not as luminous in this case, or perhaps the object is in an ‘intermediate’ stage of evolution, displaying nebular emission but not the powerful circumstellar emission of IRAS 17441-2910, nor resembling a main sequence type spectrum like IRAS 18079-1756 (see below).

As this object has previously been characterised as O6 – O7, we have searched the He ii region carefully. Fig. 8 shows this region in the telluric corrected spectrum, with the telluric standard star also plotted as a dashed line. The feature in the IRAS 17175-3544 spectrum at $``$1.693 $`\mu `$m does not appear to be an artifact of the telluric correction process, since the standard star spectrum is featureless across its profile.

We identify this feature as He ii absorption. The equivalent width is close to 300 mÅ, compatible with the measurements of Hanson et al. (1998): in their small sample, they measure 0.5$`\pm `$0.2 Å for their latest O-dwarf (HD48279, O8V) and 0.3$`\pm `$0.2 Å for the O7V HD 47839, while they do not detect He ii in their O9.5V object.

Moreover, Meyer et al. (1998) include the O7V star HR2546 (15 Mon) in their sample. Even their best spectrum of this object is rather noisy (S/N$``$50), but clearly exhibits He ii absorption. Since the spectrum is noisy, and the continuum level around He ii uncertain, we have not attempted to plot the profile in Fig. 8. However, the equivalent width is in the range 150 – 400 mÅ, close to both Hanson et al. (1998) and our measurement. Our He ii detection therefore supports the previous classification of IRAS 17175-3544 as O6 – O7. A caveat, however, is that according to Hanson et al. (1998), He i should still be strong in absorption ($``$ 1 Å) for an O7V type. Such absorption appears to be masked by nebular emission in IRAS 17175-3544.

That IRAS 17175-3544 is located at the centre of the complex star-forming region of NGC6334 I has already been discussed. In addition to the H and He features already discussed, it shows absorption features which can be identified as Mg i (17113 Å) and the Al i triplet (16723, 16755 and 16768 Å). Usually these features are associated with cooler atmospheres, increasing in strength to K-types (Meyer et al. 1998). Inspection of the telluric standard spectrum excludes a telluric nature for this lines (unlike the emission spikes near Br 11 in IRAS 18079-1756), nor a badly determined and subtracted sky spectrum (contamination by other, late-type stars in the slit), since the slit position angles were carefully chosen to avoid any other stars along its length. However, a close examination of the VLT acquisition image for IRAS 17175-3544 reveals the presence of a second object, which would have been in the slit at the time the spectrum was obtained. The two sources are not resolved in the 2MASS image.

These two sources are coincident with IRS1-I 1E and IRS-I 1W of Persi et al. (1996, see also Tapia et al. 1996) where the former source is shown by these authors to be the reddest (H = 15, H-K = 3.9), dominating the emission longward of 2 $`\mu `$m and coincident with the peak of the 30 $`\mu `$m emission. This source is the OB star responsible for the ionization of the compact H ii region and for pumping the nearby H2O and OH masers. The second source has H = 14.4 and J-H = 2.3, and clearly also contributes to the H-band spectrum, giving rise to the composite spectrum we observe. The presence of the Mg i and Al i lines, presumed to arise from IRS1-I 1W, shows that this source is cooler. We have performed experiments with synthetic spectra which show that the flux emitted by a Teff = 20000 K photosphere is only $``$ 5 times greater than one of Teff = 5000 K in the H-band, so only a very few cooler embedded stars might suffice to cause the appearance of Mg i and Al i. In this scenario, IRS-I 1W might be a system of less massive embedded objects associated with IRS-I 1E, the embedded OB star which gives rise to the H and He features in our spectrum.

Lastly, it might be argued that the observed spectrum could be that of a foreground late-type object, with emission lines from the nebular background. We cannot definitively rule out this possibility using the data presented here, but would point out that there is no optical counterpart to such a putative foreground object, even on the very deep SuperCosmos (UKST I-band, Hambly et al. (2001) and references therein) images.

### 4.2 IRAS 17441-2910

The very high signal-to-noise ($``$150) spectrum of this bright (H=10.29) source is shown in Fig. 5 (central large panel). The spectrum exhibits broad, strong emission lines of H i Br 10, Br 11 and numerous other species, notably He i and \[Fe ii\]. The emission lines we have been able to identify are labelled in Fig. 5. The object is known to exhibit very strong CO emission in the K-band (Porter et al. (1998)). The equivalent widths of the Br 11 and Br 10 emissions are $``$6 Å and $``$5.5 Å respectively, with FWHM $``$180 kms-1 and 160 kms-1. Porter et al. (1998) find Br$`\gamma `$ in emission with equivalent width 12 Å. Walsh et al. (1997) have derived a spectral type of O5.5 from the IRAS fluxes, together with the kinematic distance given by the methanol maser (6.67 GHz) emission profile (9.8 kpc), and assuming the source is a single star. However, the lack of photospheric absorption lines precludes a spectral type determination here.

After correction to the LSR velocity, we derive a radial velocity of +44.9 $`\pm `$ 4 km s-1 using 4 lines labelled in Fig. 5 which we are able to deblend accurately: Br 11 16811.11, Br 10 17366.85, Fe ii 16877.81 and Fe ii 17115.95 222Rest wavelengths from www.pa.uky.edu/$``$peter/atomic/. This velocity is somewhat different than that cited by Walsh et al. (1997), and again may indicate that the atomic emission arises in a different region from that producing the methanol maser emission (centred near +23 km s-1). However, owing to the large FWHM of the emission lines measured, and consequent difficulties in accurately measuring the line centre, we do not consider this point further.

Kaper et al. (2002) have drawn attention to a significant population of objects whose K-band spectra also show broad, strong emission lines, up to 100 kms-1, two of which also show CO emission. IRAS 17441-2910 clearly belongs in this category also. Following the suggestion of Kaper et al. (2002), IRAS 17441-2910 appears to be a third object for which strong emission of CO and other species, together with a significant infared excess (from 2MASS colours), indicate that the exciting object is surrounded by dense circumstellar material. It is clear that the strong emission line objects represent massive stars in the earliest stages of evolution.

### 4.3 IRAS 18079-1756

Previous literature on this object gave Lbol $``$ 32000 L☉ and Av = 74 (Osterloh et al. 1997). This bolometric luminosity would correspond roughly to an early B spectral type. The entire H-band spectrum is shown in Fig. 5 (lower large panel). The plot shows that the continuum slope for IRAS 18079-1756 rises to the red across the length of the H-band spectrum by $``$20%, comparable to IRAS 17441-2910, suggesting it is indeed deeply embedded. However, the absence of significant emission, and the general appearance of the spectrum (similar to an unobscured object of the same type) suggests that it could be at a somewhat later evolutionary phase.

The spectrum clearly shows broad He i absorption, which is centred at $``$ 17007Å, as well as broad hydrogen absorption. In Fig. 9, the normalised spectrum of IRAS 18079-1756 is compared to an H-band spectrum of HR 5191 ($`\eta `$ UMa). This star is one observed by Meyer et al. (1998), who give a spectral type of B3V and vsini = 205 kms-1. The He i profile is well matched, but the H lines are weaker than in the standard star, presumably because the H profiles are filled in by emission. This is supported by inspection of the Br 10 profile, which clearly shows central emission. The He i profile is contaminated by a telluric line, since it was not possible to correct for telluric absorption over a 50Å region across the He i profile (see Sect. 2.2). Fig. 10 (left panel, dashed line) shows the blended He + telluric profile in the B-type standard observed for IRAS 18079-1756. In the same way as for IRAS 17175-3544, we have synthesised this spectrum (see Fig. 6) adopting parameters of Teff = 25000K, log $`g`$ = 4.0 and vsin$`i`$ = 100 km s-1 (dotted line). Division yields a pure telluric profile (solid line). In Fig. 10 (right panel), this resultant profile (dotted line) is divided through the uncorrected spectrum for IRAS 18079-1756 (dashed line) to yield the corrected He i profile (heavy line). This profile matches the broad feature in Fig. 9 well, supporting the case that the stellar He i line is the broader feature. Comparing this profile to the galactic standard yields a very close match between the two profiles. It is clear that this match strongly supports an early-B spectral type for IRAS 18079-1756, and furthermore suggests it is a rapid rotator, as expected for very young massive YSOs. However, a rigorous measurement of vsini would require spectra of rather higher resolution.

After telluric correction, the equivalent width of He i (heavy profile in Fig. 10) was measured as 600 $`\pm `$ 100 mÅ. Comparison with equivalent widths tabulated by Hanson et al. (1998) against optical spectral type, suggest IRAS 18079-1756 is very close to B2, assuming a dwarf luminosity. Comparison with HR 5191 suggests B3V. However, strong He i persists in supergiants to $``$B7. If IRAS 18079-1756 is of lower gravity, which would not be unexpected by analogy with low-mass embedded protostars (see Greene & Lada 1996), then a somewhat later spectral type is possible. The equivalent width of Br 11 is $``$2.0$`\pm `$0.1 Å, which also suggests early B, although because of the possibility of emission this is not such a reliable indicator.

Furthermore, He i 1.700 $`\mu `$m remains strong into the O-star regime, as far as O7, and the Br 11 widths given by Hanson et al. (1998) for their late O-stars show large scatter, varying non-monotonically from 1.6 to 3.0Å between O9.5V and O7V. If IRAS 18079-1756 is as early as O8V, however, it should show He ii absorption at 1.693 $`\mu `$m with a strength comparable to that of the observed He i profile. Such a feature is not seen. Hanson et al. (1998) do not detect He ii in their O9.5V object, so in principle IRAS 18079-1756 could be as early as this. Using slightly lower resolution (R=2000), comparable S/N spectra, they give equivalent widths as small as 0.3 Å reliably. Our upper limit on the strength of He ii in IRAS 18079-1756 is smaller than this ($``$150 mÅ: we marginally detect a feature of this strength close to the expected He ii wavelength, but it is compatible with the noise). Therefore, it is not possible to say definitively that IRAS 18079-1756 is not an embedded very late O-type photosphere. Nevertheless, since we do not detect He ii with certainty, we adopt a spectral type of B2-B3, provided the luminosity class is nearer that of a dwarf than a supergiant, which seems likely.

As in the case of IRAS 17175-3544, IRAS 18079-1756 also shows a composite spectrum, with the Al i and Mg i lines. The clear presence of a cluster around the latter object suggests that this spectrum may be composite for similar reasons to IRAS 17175-3544; namely the existence of a small cluster of unresolved low-mass embedded objects around or associated closely with the massive star. For IRAS 18079-1756, there is no suggestion of a second object in the slit which can be identified with a known object, unlike the case of IRAS 17175-3544 (IRS1-I 1W). The acquisition image PSF exhibits some asymmetry, but how many components it contains is unclear. Once again, we note that if the Mg i and Al i lines are to be explained by a chance superposition of a late-type foreground object, there is no optical counterpart at the correct position on the SuperCosmos UKST I-band image.

## 5 Discussion

### 5.1 The nature of the three sources

Our sample of three embedded massive YSOs in UC H ii regions has uncovered a variety of objects which may represent different stages of the earliest evolution of massive stars.

IRAS 17441-2910 shows broad, strong emission lines, and is known to show strong CO emission, indicating the presence of dense circumstellar material. While the intrinsic near-infrared excess emission is marginal compared to that of IRAS 17175-3544, the presence of emission lines suggests mass accretion and/or outflow processes. This object appears to be at the earliest evolutionary phase of the three candidate massive YSOs.

IRAS 17175-3544 is in some ways the most enigmatic source, since it has a flat continuum, but the largest near-IR excess of all. Its association with IRS-I 1 in the molecular cloud complex NGC 6334 (see Sect. 3) supports a young evolutionary status. The nebular emission detected (Br 10, Br 11, He i) is narrower and weaker than in IRAS 17441-2910, but still clear. If the rather flat continuum of IRAS 17175-3544 is taken to mean it is relatively less embedded than IRAS 17441-2910, the infrared excess is not explained. Instead, the weakness of the emission lines compared to IRAS 17441-2910 may be indirect evidence that IRAS 17175-3544 is less luminous than IRAS 17441-2910, and that the spectral type of the latter is earlier than O6 – O7, which is supported by the previous derivation of a spectral type of O5.5 (see Table 1). Note however, that the broad emission in IRAS 17441-2910 is most likely of circumstellar rather than nebular origin. We propose that IRAS 17175-3544 may be at an equally early stage of evolution as IRAS 17441-2910.

Finally, the spectrum of IRAS 18079-1756 resembles a main sequence B-star, lacking direct spectral evidence for circumstellar material. The absence of near-infrared emission line features suggests a later evolutionary stage than the other two objects. However, the continuum does slope to the red, rising by $``$ 20% over the 750 Å range covered by these spectra, similar to the slope exhibited by IRAS 17441-2910. Perhaps tellingly, IRAS 18079-1756 shows little near-infrared excess, while IRAS 17441-2910 does. This supports the hypothesis that IRAS 18079-1756 is more evolved but still obscured, possibly predominantly by intracluster dust.

While this would explain the slope of the continuum and the location of the source in Fig. 2, the association of the source with a UC H ii region and H2O masering, and its identification as a CO outflow source by Osterloh et al. (1997), lend weight to the suggestion that it is an embedded (that is, still surrounded by circumstellar material) YSO.

We note also that the very presence of emission lines does not confirm these sources are massive embedded YSOs; OB supergiants such as P Cygni also exhibit emission lines in the near-infrared (Hanson et al. 1998), but such objects are not generally optically obscured, nor associated with masering phenomena.

### 5.2 Spectral types from H-band data

Where H and He absorption lines are detected, it is possible to obtain unambiguous spectral types, by comparison to galactic standards. He ii absorption is detected in IRAS 17175-3544, in close quantitative agreement with the spectral type derived from IRAS fluxes, O6 – O7. From the presence of strong He i and the absence of detectable He ii in IRAS 18079-1756, we also confirm a spectral type of B3, with considerable accuracy ($`\pm `$ $``$3 subclasses). These close agreements strongly suggest that near-infrared spectral types could be derived for many of the known UC H ii regions for which far-infrared spectral types, derived from bolometric luminosities, are ambiguous (e.g. Walsh et al. 1997).

In a very recent paper, Hanson et al. (2002) have performed a large survey of UC H ii regions in the K-band, including spectroscopy at R=1200. Their main finding relevant here is that $``$ 50% of their radio-selected sample were detectable by Br$`\gamma `$ emission, but of those, only for 5 – 10% was it possible to detect photospheric features of the ionizing star, and thus derive spectral type and hence effective temperature directly. However, in the K-band spectroscopic survey of Kaper et al. (2002), at R$``$8000, these authors detect characteristic K-band OB star features (emission and absorption) in the majority of their sample of 75, of which $``$ $`\frac{1}{3}`$ are found to be late-type foreground contaminants, and around 20 are found to be extremely strong Br$`\gamma `$ emitters with obvious resemblances to IRAS 17441-2910, as discussed above. This suggests that if sufficiently high resolution and high signal-to-noise spectra (R$``$5000, S/N$``$100) are used in either the H or K band, the photospheres of massive YSOs, at very early stages of evolution, could be readily detected, and absorption lines analysed. Such spectra would be amenable to similar techniques as those common in optical studies, in particular the derivation of gravities log g and projected rotational velocities vsini, both of which may be expected to differ systematically from that observed on the main sequence. Finally, given suitable line(s) present in the H-band, the detection of a surface enhancement of 14N in young massive stars suggests an observational test of the coalescence formation scenario, by analogy with O-type blue stragglers (Kendall et al. (1995)).

## 6 Conclusions

1. We have used high signal-to-noise, intermediate resolution H-band spectroscopy to probe the central, ionizing stars of three UC H ii regions. In two of the objects, IRAS 18079-1756 and IRAS 17175-3544, we detect the photosphere directly via absorption lines, allowing the derivation of spectral types which are in close agreement with those found by far-infrared observations.
2. In a third object, we detect no absorption, but a rich emission line spectrum. We identify this object, IRAS 17441-2910, as an additional member of the small subset of massive YSOs displaying strong emission (including CO) and which are likely to be at an extremely early evolutionary stage, surrounded by dense circumstellar material. The large infrared excess of IRAS 17175-3544, coupled with the observation of emission lines, also places it in this category.
3. IRAS 17175-3544 and IRAS 18079-1756 exhibit composite spectra, including absorption lines of neutral metals normally observed in cool stars, as well as H and He lines arising in the OB photospheres of or nebulosity associated with the targets themselves. Both these objects are associated with small clusters, as observed in our acquisition images. For IRAS 17175-3544, we identify the source of the metal lines as IRS-I 1W (Persi et al. 1996) and suggest that this object is a system of low mass embedded objects associated with the massive source itself. For IRAS 18079-1756 we tentatively suggest a similar scenario, although there is no direct evidence from imaging. In neither case can the spectra presented here rule out a chance superposition of a foreground late-type object causing the composite spectra, but if such a foreground star exists, it does not appear on SuperCosmos UKST I-band images of the two near-infrared sources.
4. H or K-band spectra at a comparable or slightly higher resolution to those in this study, with high signal-to-noise ratio, may be amenable to quantitative analysis; e.g. derivation of projected rotational velocities and perhaps surface gravities, in the same way as for optical spectra. This will open up a new method to study massive YSOs at the earliest stages of their formation. Both bandpasses are likely to be useful: sources will be brightest in the K-band, where CO is also observed, but the H-band promises greatly reduced veiling and is likely to be a better wavelength region to seek weak metal absorption lines.

###### Acknowledgements.

TRK and WJdeW acknowledge financial assistance from the European Union Research Training Network ‘The Formation and Evolution of Young Stellar Clusters’ (RTN1-1999-00436). The authors thank the referee, Lex Kaper, for exhaustive comments and discussions. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.


format: text

text:
# Toward one-band superconductivity in MgB2

## Abstract

The two-gap model for superconductivity in MgB2 predicts that interband impurity scattering should be pair breaking, reducing the critical temperature. This is perhaps the only prediction of the model that has not been confirmed experimentally. It was previously shown theoretically that common substitutional impurities lead to negligible interband scattering—if the lattice is assumed not to distort. Here we report theoretical results showing that certain impurities can indeed produce lattice distortions sufficiently large to create measurable interband scattering. On this basis, we predict that isoelectronic codoping with Al and Na will provide a decisive test of the two-gap model.

It is now widely accepted that MgB2 is a two-gap superconductor: its Fermi surface consists of two distinct sheets characterized by strong and weak electron-phonon coupling, respectively (see Ref. mazin03,  for a review). This view is supported by numerous experiments probing either the larger or smaller gap, or both simultaneously. Experimental observation of the merging of the two gaps would constitute even stronger evidence. Such a merging is expected, for example, from interband scattering by impurities. Within the theory of multiband superconductivity, interband scattering mixes the “weak” and “strong” Cooper pairs, averaging the order parameter and reducing $`T_c`$ golubov97 . At small defect concentrations the suppression of $`T_c`$ should be linear, with the larger gap decreasing and the smaller gap increasing. The effect should be pronounced in samples with high defect concentrations, but despite the relatively low quality of many samples such an effect has not been observed. Indeed, some samples with high resistivity have nearly the same critical temperatures as clean single crystals. On the other hand, samples doped with carbon show both gaps decreasing but not merging, despite a considerable reduction in $`T_c`$. Other types of intentional defects—whether from doping or irradiation—have also failed to produce a merged gap. Even efforts to introduce defects into MgB2 for the explicit purpose of inducing interband scattering and merging the two gaps have failed to observe this effect wang02 ; schmidt03 ; samuely03 ; balaselvi03 .

This apparent lack of evidence for a central prediction of the two-gap model is disturbing. In Ref. mazin02  it was shown that interband scattering from substitutional impurities is inherently weak, if the lattice is assumed not to distort. This partially explains the null results of current experimental efforts to induce pair breaking, but it does not address the possibility of pair breaking from impurities specifically chosen to maximize the interband scattering due to large lattice distortions. In this Letter we show that this strategy is likely to succeed and, by identifying a simple relationship between impurity atoms and the resulting lattice distortions in MgB2, suggest an impurity-doping protocol that will produce measurable pair breaking—and thereby provide the final “smoking-gun” evidence for the two-band model.

Hampering such investigations is the currently limited insight into which defects are most effective in creating interband scattering. Since the states at the Fermi level of MgB2 are formed by the boron orbitals, one might expect impurities (such as carbon) in the B planes to produce large interband scattering. This is not borne out experimentally: substitutional C impurities have only a weak effect on the interband scattering samuely03 . This finding had been anticipated theoretically as a consequence of the special symmetry properties of the electronic states within the $`\sigma `$ band near the $`\mathrm{\Gamma }`$ point mazin02 . The crucial point is that although impurities in the B plane do have a strong effect on the electronic structure, they do not change the local point symmetry and therefore do not lead to significant $`\sigma \pi `$ scattering.

The situation can be quite different for substitutions in the Mg plane, which may create out-of-plane distortions of B atoms in neighboring planes. Such relaxations do change the local point symmetry of nearby B atoms, mixing the in-plane $`p_{x,y}`$ and out-of-plane $`p_z`$ orbitals, and for sufficiently large disturbances can lead to significant $`\sigma \pi `$ scattering. Here we demonstrate by first-principles calculations that this is indeed the case for certain substitutional impurities but—surprisingly—not for Mg vacancies. We predict that the interband scattering effects will be most pronounced for isoelectronic co-doping with Na and Al, and that the effect on the superconducting properties should be detectable for impurity concentrations above 2%.

We used density-functional theory (DFT) to study the lattice distortion created by Mg-plane substitutional impurities from Groups I, II, and III, by a Mg vacancy, and by B-plane C substitution. To model the distortion induced by single defects, we used 2$`\times `$2$`\times `$2 supercells of bulk MgB2 (and 3$`\times `$3$`\times `$3 supercells for convergence checks). Total energies and forces were calculated using projector-augmented-wave potentials and the generalized-gradient approximation kresse93a ; kresse96a . All atomic positions were relaxed within the constraint of fixed (theoretical) bulk lattice parameters. The resulting displacements of the nearest-neighbor B atoms are given in Table I. Given the strong intraplanar covalent bonding, it is not surprising that these displacements are dominated by the out-of-plane component $`\delta z`$, which ranges from $``$0.04 Å (for Be) to $`+`$0.09 Å (for K); for C it is zero by symmetry.

The out-of-plane displacements for seven different impurities on the Mg site are plotted in Fig. 1. In a previous related study, changes in interlayer spacing induced by the complete substitution of Al for Mg in a single plane were ascribed to electrostatic effects barabash02 . We find no such correlation between $`\delta z`$ and the formal valence of the impurity—for example, Be and Ca give displacements of opposite sign despite having identical valence—and thus infer that electrostatic effects are not important. On the other hand, Fig. 1 shows that there is an excellent correlation between $`\delta z`$ and the ionic radius of the impurity atom. Hence, we conclude that the out-of-plane displacement of B atoms by impurity atoms substituting for Mg is mostly a size effect.

In light of this, one might anticipate a Mg vacancy to produce a large inward displacement. Our DFT results reveal very different behavior: the vacancy creates a negligibly small displacement, $`\delta z<0.01`$ Å. This result is consistent with the experimental fact that pair breaking is not observed in low quality samples, which presumably contain many vacancies. However, it is very different from the trend shown in Fig. 1. Indeed, if one naively considers the vacancy as an impurity of zero size, the predicted displacement is almost twice that found for Be and Al, in sharp distinction to our DFT result.

A simple model explains the surprisingly small displacement created by the Mg vacancy. We consider a MgB2 crystal containing a single substitutional defect $`D`$ (either an impurity atom or a vacancy) in the Mg plane. Such a defect has 12 nearest-neighbor B atoms, consisting of two hexagonal rings. We consider the out-of-plane displacement of the B atoms in these rings to arise from two opposing effects. The first of these represents the change in covalent bonding between B planes. In the absence of any defects, the MgB2 interlayer spacing $`c_{\mathrm{Mg}}`$ is primarily determined by assisted hopping between B $`p_z`$ orbitals through Mg orbitals; for each B atom there are three such hopping paths through nearest-neighbor Mg atoms. With the defect present, one of these three paths now passes through the defect site. This new hopping path results in an out-of-plane force on the B atom. We assume the magnitude of this effect to be one-third of that found for a fully substituted Mg plane, which we approximate using the energy vs. layer spacing, $`E_D(c)`$, for fully substituted $`D`$B2. Thus, we consider the change in spacing between the two displaced hexagonal rings, $`2\delta z`$, to contribute an energy per B given by $`\frac{1}{3}E_D(c_{\mathrm{Mg}}+2\delta z)`$. We have calculated the binding-energy curves $`E_D(c)`$ within DFT for AlB2, NaB2, and “vacancy-substituted” B2; the results are shown in Fig. 2. All can be accurately represented by a Morse potential, which is the form we will use in the discussion below.

The second effect is the restoring force experienced by the displaced B atoms, due to the strong covalent bonding within the B planes. For the small displacements we are considering it is reasonable to take this effect to be harmonic in $`\delta z`$, again weighted by 1/3. Thus we take the total energy change per B to be

$$E(\delta z)=\frac{1}{3}K(\delta z)^2+\frac{1}{3}kw^2[1e^{(c_{\mathrm{Mg}}+2\delta zc_D)/w}]^2,$$
(1)
where $`w`$ is the width of the Morse potential, $`kw^2`$ its depth, and $`c_D`$ its equilibrium interlayer spacing. For small $`c_{\mathrm{Mg}}c_D`$, it is easy to show that this energy is minimized for $`2\delta z=(c_Dc_{\mathrm{Mg}})/(1+K/2k)`$. In other words, for impurities whose size is comparable to Mg, the displacement is linear in the size mismatch and reduced by the factor $`1+K/2k`$ (which is typically in the range 5–10); this is consistent with the DFT results shown in Fig. 1. Qualitatively, when $`c_D`$ is close to $`c_{\mathrm{Mg}}`$, as it is for AlB2 and NaB2, the equilibrium displacement represents a balance between the harmonic restoring force $`\frac{2}{3}K|\delta z|`$ and a nearly harmonic (attractive or repulsive) force $`\frac{2}{3}k(c_{\mathrm{Mg}}c_D+2\delta z)`$ from the interlayer bonding.

For the vacancy there is a very large size mismatch: $`c_{\mathrm{vac}}`$ is over 40% smaller than $`c_{\mathrm{Mg}}`$. Thus for any reasonable displacement, the hexagons experience only a weak attractive force from the tail of the Morse potential. Hence, the energy is minimized for a very small displacement, $`\delta z0.007`$ Å, in agreement with the negligible displacement given by DFT. Qualitatively, the result of a large size mismatch is to largely preempt the mechanism of interlayer binding, leading to very small displacements strongly suppressed by the penalty for perturbing the planarity of the B layer.

The origin of the large mismatch between $`c_{\mathrm{Mg}}`$ and $`c_{\mathrm{vac}}`$ can best be understood by comparing the band structures of the fully substituted $`D`$B2 compounds at their equilibrium interlayer spacings. For these pure compounds (as well as the parent material) interlayer bonding arises primarily from the interaction between B $`p_z`$ orbitals in different layers. This interaction depends on the assisted hopping through $`s`$ and $`p_z`$ orbitals located on the $`D`$ site. However, for the fully “vacancy substituted” compound the occupancy of the $`\pi `$ bands is so much reduced that their contribution to bonding becomes quite small. At the same time, the $`\sigma `$ bands acquire substantial $`z`$-dispersion from $`pp\pi `$ hopping, which contributes to bonding. Even for high vacancy concentrations (without complete removal of a Mg plane) it is impossible to engage the $`\sigma `$ bands in interlayer bonding by any reasonable dimpling of the B planes, because the planes remain too far apart. As a result, hardly any distortion occurs at all: indeed, even for 50% vacancies within a single Mg plane, the interlayer spacing changes by less than 0.05 Å.

We have established that the defects leading to the largest B displacements are Mg-plane substitutional impurities with a large size mismatch (but not Mg vacancies). We now estimate the magnitude of the interband scattering associated with displacements from such impurities. We assume that the only relevant scattering is that due to the out-of-plane distortion, $`\delta z`$, and use the analogy between the formulas allen78  for the impurity-induced scattering rate,

$$\gamma _{\mathrm{imp}}=\pi n_{\mathrm{imp}}\frac{_{\mathrm{𝐤𝐤}^{}}\delta (\epsilon _𝐤)\delta (\epsilon _𝐤^{})|\delta V_{\mathrm{𝐤𝐤}^{}}|^2}{_𝐤\delta (\epsilon _𝐤)},$$
(2)
and the electron-phonon coupling constant,

$$\lambda =\frac{_{\nu ,\mathrm{𝐤𝐤}^{}}\delta (\epsilon _𝐤)\delta (\epsilon _𝐤^{})|M_{\nu ,\mathrm{𝐤𝐤}^{}}|^2/\mathrm{}\omega _{\nu ,𝐤𝐤^{}}}{_𝐤\delta (\epsilon _𝐤)}.$$
(3)
Here $`n_{\mathrm{imp}}`$ is the impurity concentration; $`\delta V_{\mathrm{𝐤𝐤}^{}}`$ is the matrix element of the impurity perturbation potential (defined as the difference between the full crystal potentials with and without an impurity); $`\epsilon _𝐤`$ is the electron energy with respect to the Fermi level; $`\omega _{\nu 𝐤}`$ is the phonon frequency; $`M_{𝐤,𝐤^{}}`$ is the electron-ion matrix element $`𝐤|dV/dQ_\nu |𝐤^{}`$, where $`dV/dQ_\nu `$ is the derivative of the crystal potential with respect to the phonon normal coordinate $`Q=\sqrt{2m\omega /\mathrm{}}x`$; and the summations are over all electron states and all phonon branches allen78 .

To proceed we make three approximations, all qualitatively reasonable if not quantitatively reliable. First, we assume that in-plane phonons, including the well-known $`E_{2g}`$ modes, contribute little to interband electron-phonon coupling; this follows from the same symmetry arguments given in Ref. mazin02, . Second, we assume by the same reasoning that interband impurity scattering comes only from the out-of-plane relaxation of B atoms. Finally, we approximate Eq. 2 as

$$\gamma _{\mathrm{imp}}12\pi n_{\mathrm{imp}}\left(\underset{𝐤}{}\delta (\epsilon _𝐤)\right)\left|dV/du_z\right|^2(\delta z)^2,$$
(4)
and likewise Eq. 3 as

$$\lambda 2\left(\underset{𝐤}{}\delta (\epsilon _𝐤)\right)\left|dV/du_z\right|^2/2m\omega ^2.$$
(5)
Here we have assumed that the average of the crystal potential with respect to the vertical displacement of the B atom is the same in both cases. The numerical factors 12 and 2 are the coordination of Mg and the number of B atoms in the unit cell, respectively.

First-principles calculations give the interband part of the electron-phonon coupling as $``$0.2 mazin03 . The phonon frequencies for the out-of-plane modes are about 400 cm$`^1.`$ Using these values, the scattering rate is given by $`\gamma _{\mathrm{imp}}50n_{\mathrm{imp}}(\delta z)^2`$ eV, where $`\delta z`$ is in angstrom. Hence, for 2% Al doping we find $`\gamma _{\mathrm{imp}}1.1`$ meV. For 2% Na doping we estimate a similar scattering rate, $`\gamma _{\mathrm{imp}}1.3`$ meV. These scattering rates are small, but should still have a measurable effect on the superconducting gaps and temperature. The effect on the gaps is difficult to estimate without full Eliashberg calculations. The reduction of $`T_c`$ can be easily estimated using Eq. 13 from Ref. golubov97 , which gives 2.0–2.5 K (in addition to any suppression due to the electron doping of the $`\sigma `$ band). While this is a small reduction compared to the changes observed in heavily electron-doped samples, the underlying mechanism is quite different. In particular, pair breaking from interband scattering is unique in that it reduces $`T_c`$ and the gap ratio while simultaneously increasing the smaller gap golubov97 . This distinctive behavior should facilitate the separation of interband scattering from other sources of $`T_c`$ reduction.

Finally, we suggest that an especially attractive test of these predictions would be simultaneous codoping by equal parts Al and Na. This would effectively be an isoelectronic substitution, and any effect on the superconducting properties could then be ascribed to impurity-induced interband scattering. Moreover, Na and Al induce distortions of the same magnitude but of the opposite sign, which should mitigate the effects of a reduction in scattering due to possible short-range ordering of the impurities.

In conclusion, we have performed first-principles calculations of the lattice distortion in MgB2 induced by several common substitutional (for Mg) impurities, and a Mg vacancy. We find out-of-plane displacements as large as 0.04 Å for common impurities such as Al and Na. The magnitude and sign of the displacement are mainly determined by the ionic size of the impurity. On the other hand, for the Mg vacancy we find an essentially negligible displacement of nearby B atoms. The different behavior of vacancies and impurities is explained by a simple physical model representing the competition between interlayer binding and intralayer planarity. We estimate the interband scattering rate due to the Na and Al impurities to be of order 1 meV, sufficiently large to give a detectable change in the superconducting transition temperature. Finally, we propose that the codoped material Mg1-xNaxAlxB2, which is isoelectronic with MgB2, should provide an excellent test of these predictions.

We are grateful to O. K. Andersen for many stimulating discussions. This work was supported by the Office of Naval Research.


text:
# A Stochastic Model for the Evolution of the Web Allowing Link Deletion

## 1 Introduction

Power-law distributions taking the form

$$f(i)=Ci^\tau ,$$
(1)
where $`C`$ and $`\tau `$ are positive constants, are abundant in nature \[Sch91\]. The constant $`\tau `$ is called the exponent of the distribution. Examples of such distributions are: Zipf’s law, which states that the relative frequency of words in a text is inversely proportional to their rank, Pareto’s law, which states that the number of people whose personal income is above a certain level follows a power-law distribution with an exponent between 1.5 and 2, Lotka’s law, which states that the number of authors publishing a prescribed number of papers is inversely proportional to the square of the number of publications, and Gutenberg-Richter’s law, which states that the number of earthquakes over a period of time having a certain magnitude is roughly inversely proportional to the magnitude.

Recently several researchers have detected asymptotic power-law distributions in the topology of the World-Wide-Web \[BKM+00, DKM+02\] and, in parallel, researchers from a variety of fields are trying the explain the emergence of these power laws in terms of the evolution of complex networks. See \[AB02, DM02, New03\] for comprehensive surveys of recent work in this area, detailing various stochastic models that can explain the evolution of the web and other networks such as the internet, citation networks, collaboration networks and biological networks. A common theme in these models is that of preferential attachment, which results in the “rich get richer” phenomenon, for example, when new links to web pages are added in proportion to the number of currently existing links to these pages. Related approaches are the general theoretical model covering both directed and undirected web graphs \[CF03\], the stochastic multiplicative process in which nodes appear at different times and the rate of addition of links varies between nodes \[AH01\], and the statistical physics approach that uses the rate equation technique \[KR02\].

To explain situations where pure preferential attachment models fail, we \[LFLW02\] and others \[PFL+02\] have previously proposed extensions of the stochastic model for the web’s evolution in which the addition of links is prescribed by a mixture of preferential and non-preferential mechanisms. In \[LFLW02\], we devised a general stochastic model involving the transfer of balls between urns; this also naturally models quantities such as the numbers of web pages in and visitors to a web site, which are not naturally described in graph-theoretic terms. We note that our urn model is an extension of the stochastic model proposed by Simon in his visionary paper published in 1955 \[Sim55\], which was couched in terms of word frequencies in a text. We also considered an alternative extension of Simon’s model in \[FLL02\] by adding a preferential mechanism for discarding balls from urns (corresponding to deleting web pages); this results in an exponential cutoff in the power-law distribution.

Our urn transfer model is a stochastic process, in which at each step with probability $`p`$ a new ball (which might represent a web page) is added to the first urn and with probability $`1p`$ a ball in some urn is moved along to the next urn. We assume that a ball in the $`i`$th urn has $`i`$ pins attached to it (which might represent web links). It is known that the steady-state distribution of this model is a power law, with exponent $`\tau =1+1/(1p)`$ \[Sim55, LFLW02\]. As mentioned above, the power-law distribution breaks down when balls may be discarded, resulting in a power-law distribution with an exponential cutoff. So the question arises as to whether it also breaks down when removal of pins is allowed. We answer this question here by showing that the power-law distribution does not break down, under the constraint that more balls are added to the system than are removed. (When the only remaining pin is removed from a ball in the first urn that ball is removed from the system.) This model gives a more realistic explanation for the emergence of power laws in complex networks than the basic model without link deletion (i.e. pin removal). Considering, for example, the web graph, the modification we make to the basic model is that after a web page is chosen preferentially, say according to the number of its inlinks, there is a small probability that some link to this page will be deleted. (This is equivalent to deleting a link chosen uniformly at random.) A possible reason for this may be that popular web pages compete for inlinks, so that the number of inlinks a page has acquired will fluctuate with its popularity. This is evident, for example, in the rise and fall in the popularity of several search engines. Link deletion has also been considered by Dorogovtsev and Mendes \[DM01\] in the context of a different model in which, at each time step, a new node is added with just one link, which preferentially attaches itself to an existing node. In addition, at each time step $`m`$ links between existing nodes are added to preferentially chosen nodes and $`c`$ links between existing nodes are deleted, again preferentially. Their conclusion that link deletion increases the power-law exponent is also obtained here for our stochastic model.

Consider a power law such as Lotka’s law \[Nic89\]. If this holds, then a plot on a log-log scale of the number of authors (the frequency) against the number of publications (the value) should reveal a straight line with a negative slope of around $`2`$. There is an obvious problem with a log-log transformation if any of the frequencies are zero; that is to say, when there are values $`v_1`$ and $`v_2`$, with $`v_1<v_2`$, for which no author published $`v_1`$ papers but at least one author published $`v_2`$ papers (i.e. there is a gap in the values for the number of publications).

In general, we expect such gaps to occur in any data set, mainly for large values. (This is due to the stochastic variation and the fact that the frequencies have to be integral.) This observation is consistent with the finding that Lotka’s law does not fit well for large values (i.e. authors who published a large number of papers) and this tail region is characterised by the presence of gaps \[Nic89\].

One way of dealing with the problem of gaps is simply to ignore all value-frequency pairs where the frequency is zero. However, ignoring gaps in this way seems questionable, since the zero frequency values do give relevant information about the data set and should not be treated as missing values. Our unease with this approach is reinforced by the fact that computing the exponent in this way results in a much lower value; for example, for the web inlink data \[BKM+00\], it gives a value around $`1.5`$, rather than the generally accepted figure of $`2.1`$.

An alternative approach is to squash the non-zero frequencies up towards the y-axis (i.e. the frequency axis), after ignoring the zero frequency values. This is equivalent to omitting from the data set those values for which there are no authors having that number of publications and then ranking the remaining values in increasing order, i.e. renumbering the values starting from one. This method also seems somewhat dubious, since the power-law relationship should really involve the actual values rather than their ranks.

The standard technique for fitting a power-law distribution uses linear regression on log-log transformed data, which is not possible if gaps are present. Another approach for handling gaps is to consider only the values up until the first gap; this approach is only reasonable if the first gap does not occur at too low a value. We call this the unranked approach, and note that this approach ignores the large values in the tail of the distribution. The unranked approach suggests one possible solution, but other approaches, such as smoothing or using the Hill plot \[DdR00\], are possible. Ad-hoc approaches have the disadvantage that they are hard to justify, and the Hill plot, as originally defined, is only applicable after the data has been sorted, and this also seems difficult to justify.

Another solution for handling gaps is the second method described above of ranking the values with non-zero frequencies and squashing up; we call this the ranked approach. In general, we should not expect the ranked and unranked approaches to yield the same power-law exponent. We would expect the exponent computed by linear regression on the log-log transformed data from the unranked approach to be somewhat higher. One of the findings of this paper is that for inlinks in the web graph the ranked and unranked approaches lead to a small but noticeable difference between the exponents. (Our previous comments indicate why we have misgivings with the ranked approach.)

The rest of the paper is organised as follows. In Section 2 we present an urn transfer model that extends Simon’s model by allowing a pin, chosen uniformly at random, to be discarded with some specified probability. In Section 3 we derive the steady state distribution of the model, which, as stated, follows an asymptotic power law like (1). In Section 4 we show that, by examining a snapshot of the distribution of balls in urns at the steady state, we are able to tell whether removal of pins has taken place and estimate the pin removal probability. In Section 5 we utilise our urn model to describe a discrete stochastic process that simulates the evolution of the degree distribution of inlinks in the web graph. As a proof of concept, we analyse the May 1999 data set for inlinks presented in \[BKM+00\], which we obtained from Ravi Kumar at IBM. We are able to show that our model is consistent with the data and determine the extent to which link deletion has occurred. We also investigate the discrepancy between the ranked and unranked approaches. With the ranked approach, using linear regression on the log-log transformed data, an exponent of $`2.1`$ is obtained, as in \[BKM+00\]. However, the unranked approach results in an exponent of $`2.15`$, which is shown to be consistent with our stochastic model. Although the difference between $`2.1`$ and $`2.15`$ may not seem significant, it has been remarked in \[BKM+00\] that “2.1 is in remarkable agreement with earlier studies” and in \[DKM+02\] that the exponent is “reliably around 2.1 (with little variation)”, which justifies us in making a distinction between the exponents obtained using the ranked and unranked approaches. Finally, in Section 6 we give our concluding remarks.

## 2 An Urn Transfer Model

We now present an urn transfer model for a stochastic process with urns containing balls (which might represent web pages) that have pins (which might represent either inlinks or outlinks) attached to them. Our model allows for pins to be discarded with a small probability. This model can be viewed as an extension of Simon’s model \[Sim55\]. We note that there is a correspondence between the Barabási and Albert model \[BA99\], defined in terms of nodes and links, and Simon’s model, defined in terms of balls and pins, as was established in \[BE01\]. Essentially, the correspondence is obtained by noting that the balls in an urn can be viewed as an equivalence class of nodes all having the same indegree (or outdegree).

We assume a countable number of urns, $`urn_1,urn_2,urn_3,\mathrm{}`$, where each ball in $`urn_i`$ is assumed to have $`i`$ pins attached to it. Initially all the urns are empty except $`urn_1`$, which has one ball in it. Let $`F_i(k)`$ be the number of balls in $`urn_i`$ at stage $`k`$ of the stochastic process, so $`F_1(1)=1`$. Then, for $`k1`$, at stage $`k+1`$ of the stochastic process one of two things may occur:

1. with probability $`p`$, $`0<p<1`$, a new ball is inserted into $`urn_1`$, or
2. with probability $`1p`$ an urn is selected, with $`urn_i`$ being selected with probability proportional to $`iF_i(k)`$, and a ball is chosen from $`urn_i`$; then,

1. with probability $`r`$, $`0<r1`$, the chosen ball is transferred to $`urn_{i+1}`$, (this is equivalent to attaching an additional pin to the ball chosen from $`urn_i`$), or
    2. with probability $`1r`$, the ball is transferred to $`urn_{i1}`$ if $`i>1`$, otherwise, if $`i=1`$, the ball is discarded (this is equivalent to removing and discarding a pin from the ball chosen from $`urn_i`$).

In the special case when $`r=1`$, the process reduces to Simon’s original model.

We note that choosing a ball preferentially (in proportion to the number of pins) is equivalent to selecting a pin uniformly at random and choosing the ball it is attached to. Thus, with probability $`(1p)(1r)`$, a pin chosen uniformly at random is discarded.

Since $`iF_i(k)`$ is the total number of pins attached to balls in $`urn_i`$, the expected total number of pins in the urns at stage $`k`$ is given by

$`E\left({\displaystyle \underset{i=1}{\overset{k}{}}}iF_i(k)\right)`$ $`=`$ $`1+(k1)\left(p+(1p)r(1p)(1r)\right)`$ (2)
$`=`$ $`1+(k1)\left(12(1p)(1r)\right).`$
Correspondingly, the expected total number of balls in the urns is given by

$`E\left({\displaystyle \underset{i=1}{\overset{k}{}}}F_i(k)\right)`$ $`=`$ $`1+(k1)p(1p)(1r){\displaystyle \underset{j=1}{\overset{k1}{}}}\varphi _j,`$ (3)
where $`\varphi _j`$, $`1jk1`$, is the expected probability of choosing $`urn_1`$ at step (ii) of stage $`j+1`$, i.e.

$$\varphi _j=E\left(\frac{F_1(j)}{_{i=1}^jiF_i(j)}\right).$$
(4)
Now let

$$\varphi ^{(k)}=\frac{1}{k}\underset{j=1}{\overset{k}{}}\varphi _j.$$
(5)
In order to ensure that there are at least as many pins in the system as there are balls and that, on average, more balls are added to the system than are removed, we require the following constraint, derived from (2) and (3),

$$\frac{12r}{1r}\varphi ^{(k)}\frac{p}{(1p)(1r)}.$$
(6)
This implies

$$2(1p)(1r)1,$$
(7)
which obviously holds for $`r1/2`$. Inequality (7) expresses the fact that the expected number of pins should not be negative, and follows from (2).

We note that we could modify the initial conditions so that, for example, $`urn_1`$ initially contained $`\delta `$ balls, $`\delta >1`$, instead of just one ball. It can be shown, from the development of the model below, that any change in the initial conditions will have no effect on the asymptotic distribution of the balls in the urns as $`k`$ tends to infinity, provided the process does not terminate with all of the urns empty. More specifically, the probability that the process will not terminate with all the urns empty is given by

$$1\left(\frac{(1p)(1r)}{1(1p)(1r)}\right)^\delta .$$
(8)
This is exactly the probability that the gambler’s fortune will increase forever \[Ros83\]. (We note that (8) only makes sense if (7) holds.) Since, in practice, $`r`$ will be quite close to one, once the process has survived for a few steps the probability that it will subsequently terminate is small. From now on we will therefore assume that the process does not terminate.

## 3 Derivation of the Steady State Distribution

Following Simon \[Sim55\], we now state the mean-field equations for the urn transfer model. For $`i>1`$ the expected number of balls in $`urn_i`$ is given by

$$E_k(F_i(k+1))=F_i(k)+\beta _k\left(r(i1)F_{i1}(k)+(1r)(i+1)F_{i+1}(k)iF_i(k)\right),$$
(9)
where $`E_k(F_i(k+1))`$ is the expected value of $`F_i(k+1)`$, given the state of the model at stage $`k`$, and

$$\beta _k=\frac{1p}{_{i=1}^kiF_i(k)}$$
(10)
is the required normalising factor.

In the boundary case, when $`i=1`$, we have

$$E_k(F_1(k+1))=F_1(k)+p+\beta _k\left((1r)2F_2(k)F_1(k)\right),$$
(11)
for the expected number of balls in $`urn_1`$, given the state at stage $`k`$.

In order to obtain a solution for the model, we assume that, for large $`k`$, the random variable $`\beta _k`$ can be approximated by a constant value $`\widehat{\beta }_k`$ depending only on $`k`$. This is defined by

$$\widehat{\beta }_k=\frac{1p}{k\left(12(1p)(1r)\right)}.$$

The motivation for this approximation is that the denominator in the definition of $`\beta _k`$, the total number of pins, has been replaced by its expectation given in (2). This is a reasonable assumption since the number of pins is the difference between two binomial random variables, and with high probability this will be close to its expected value. We observe that using $`\widehat{\beta }_k`$ as the normalising factor instead of $`\beta _k`$ results in an approximation similar to that of the “$`p_k`$ model” in \[LFLW02\].

Replacing $`\beta _k`$ by $`\widehat{\beta }_k`$ and taking the expectations of (9) and (11), we obtain

$$E(F_i(k+1))=E(F_i(k))+\widehat{\beta }_k\left(r(i1)E(F_{i1}(k))+(1r)(i+1)E(F_{i+1}(k))iE(F_i(k))\right)$$
(12)
and

$$E(F_1(k+1))=E(F_1(k))+p+\widehat{\beta }_k\left((1r)2E(F_2(k))E(F_1(k))\right),$$
(13)
respectively.

In order to solve (12) and (13), we would like to show that $`E(F_i(k))/k`$ tends to a limit $`f_i`$ as $`k`$ tends to infinity. Suppose for the moment that this is the case, then, provided the convergence is fast enough, $`E(F_i(k+1))E(F_i(k))`$ tends to $`f_i`$. By “fast enough” we mean that $`ϵ_{i,k+1}ϵ_{i,k}`$ is $`o(1/k)`$ for large $`k`$, where

$$E(F_i(k))=k(f_i+ϵ_{i,k}).$$

Now, letting

$$\beta =k\widehat{\beta _k}=\frac{1p}{12(1p)(1r)},$$
(14)
we see that $`\widehat{\beta _k}E(F_i(k))`$ tends to $`\beta f_i`$ as $`k`$ tends to infinity. Thus, letting $`k`$ tend to infinity in (12) and (13), these yield

$$f_i=\beta \left(r(i1)f_{i1}+(1r)(i+1)f_{i+1}if_i\right),$$
(15)
and

$$f_1=p+\beta \left((1r)2f_2f_1\right),$$
(16)
respectively.

We now investigate the asymptotic behaviour of $`f_i`$ for large $`i`$. If we define $`g_i`$ to be $`if_i`$, we can rewrite (15) and (16) as

$$\frac{rg_{i1}}{g_i}+\frac{(1r)g_{i+1}}{g_i}=1+\frac{1}{i\beta }$$
(17)
and

$$\frac{p}{\beta g_1}+\frac{(1r)g_2}{g_1}=1+\frac{1}{\beta }.$$
(18)
Suppose we can write $`g_i`$ as

$$g_i=\frac{C\mathrm{\Gamma }(i+\alpha +1)}{\mathrm{\Gamma }(i+\alpha +1+\rho )}\left(1+\frac{\eta }{i^2}+O\left(\frac{1}{i^3}\right)\right),$$
(19)
where $`\alpha `$, $`\eta `$, $`\rho `$ and $`C`$ are constants.

Using the binomial theorem and the fact that $`\mathrm{\Gamma }(x+1)=x\mathrm{\Gamma }(x)`$, we get

$`{\displaystyle \frac{g_{i1}}{g_i}}`$ $``$ $`{\displaystyle \frac{i+\alpha +\rho }{i+\alpha }}\left(1+{\displaystyle \frac{\eta }{i^2}}+O\left({\displaystyle \frac{1}{i^3}}\right)\right)\left(1{\displaystyle \frac{\eta }{i^2}}+O\left({\displaystyle \frac{1}{i^3}}\right)\right)`$ (20)
$`=`$ $`{\displaystyle \frac{i+\alpha +\rho }{i+\alpha }}\left(1+O\left({\displaystyle \frac{1}{i^3}}\right)\right)`$
$`=`$ $`1+{\displaystyle \frac{\rho }{i}}{\displaystyle \frac{\rho \alpha }{i^2}}+O\left({\displaystyle \frac{1}{i^3}}\right).`$
Similarly,

$`{\displaystyle \frac{g_{i+1}}{g_i}}`$ $``$ $`{\displaystyle \frac{i+\alpha +1}{i+\alpha +\rho +1}}\left(1+O\left({\displaystyle \frac{1}{i^3}}\right)\right)`$ (21)
$`=`$ $`1{\displaystyle \frac{\rho }{i}}+{\displaystyle \frac{\rho (\alpha +\rho +1)}{i^2}}+O\left({\displaystyle \frac{1}{i^3}}\right).`$
Substituting (20) and (21) into the left-hand side of (17), we obtain

$$r\left(1+\frac{\rho }{i}\frac{\rho \alpha }{i^2}+O\left(\frac{1}{i^3}\right)\right)+(1r)\left(1\frac{\rho }{i}+\frac{\rho (\alpha +\rho +1)}{i^2}+O\left(\frac{1}{i^3}\right)\right)=1+\frac{1}{i\beta },$$

which simplifies to

$$1+\frac{\rho }{i}(2r1)+\frac{\rho }{i^2}\left((1r)(\rho +1)(2r1)\alpha \right)+O\left(\frac{1}{i^3}\right)=1+\frac{1}{i\beta }.$$
(22)
For this to hold for all large enough values of $`i`$, we require

$$\rho =\frac{1}{\beta (2r1)},$$
(23)
and

$$\alpha =\frac{(1r)(\rho +1)}{2r1}.$$
(24)
It is straightforward to obtain more terms in the expansion of (19) by a more detailed analysis, for example

$$\eta =\frac{\rho \left((1r)(\alpha +\rho +1)^2r\alpha ^2\right)}{2(2r1)}.$$

So, with $`\rho `$ and $`\alpha `$ as in (23) and (24), for large $`i`$, (19) gives an approximate solution to the recurrence defined by (15) and (16). Thus,

$$f_iCi^{(1+\rho )},$$
(25)
where $`C`$ is independent of $`i`$ and $``$ means is asymptotic to.

From (23) and (14),

$$\rho =1+\left(\frac{p}{1p}\right)\left(\frac{1}{2r1}\right).$$
(26)
If, as is usually the case, we assume that $`\rho `$ is positive, then $`r>1/2`$, and it follows that increasing $`p`$ increases $`\rho `$, whereas increasing $`r`$ decreases $`\rho `$. In particular, we observe from (26) that with pin removal, i.e. for $`r<1`$, the exponent of the power law is greater than it would be for $`r=1`$, i.e.

$$\rho \frac{1}{1p}.$$
(27)
## 4 Recognising Pin Removal

In this section we show that, using the mean-field equations derived in Section 3, we are able to detect whether pin removal has taken place by inspecting a static snapshot of the system at the steady state.

Based on (2) and (3) we have

$$\frac{pins}{k}12(1p)(1r),$$
(28)
$$\frac{balls}{k}p(1p)(1r)\varphi ,$$
(29)
where $`pins`$ and $`balls`$ stand for the expected numbers of pins and balls at stage $`k`$, and $`\varphi `$ is the asymptotic value of $`\varphi ^{(k)}`$, defined by (5), i.e. the expected proportion of pins in the first urn, $`|urn_1|/pins`$. The right-hand sides of (28) and (29) give asymptotic values of $`pins`$ and $`balls`$ as $`k`$ tends to infinity. It follows that the asymptotic value of the ratio $`balls/pins`$, which we denote by $`\mathrm{\Delta }`$, is given by

$$\mathrm{\Delta }=\frac{p(1p)(1r)\varphi }{12(1p)(1r)}.$$
(30)
From this and the fact that $`\mathrm{\Delta }\varphi `$, it follows that

$$\varphi \frac{p}{1(1p)(1r)}\mathrm{\Delta }\frac{p}{12(1p)(1r)}.$$
(31)
Using (7) and (26), this implies

$$\frac{\varphi }{2}p\mathrm{\Delta }\frac{\rho 1}{\rho }.$$
(32)
Now suppose we are given $`\rho `$, $`\mathrm{\Delta }`$ and $`\varphi `$. Then from (26) and (30) we can derive the following equations:

$$p=\frac{\varphi (\rho 1)}{2\rho (1\mathrm{\Delta })2+\rho \varphi },$$
(33)
$$r=\frac{1}{2}\left(1+\frac{\varphi }{2\rho (1\mathrm{\Delta })2+\varphi }\right).$$
(34)
Now suppose we observe an urn process with $`1/2<r1`$ and $`0<p<1`$, and we take a snapshot of the system at the steady state, obtaining empirical estimates for $`\rho `$, $`\mathrm{\Delta }`$ and $`\varphi `$ from the distribution of the balls in the urns. We would like to check whether it is possible that these values could have arisen from Simon’s process, i.e. with $`r=1`$, with the probability of inserting a ball into the first urn equal to, say, $`p^{}`$. So from (26) and (30) we would require $`\rho =1/(1p^{})`$ and $`\mathrm{\Delta }=p^{}`$. Thus, from (26), we would obtain

$$p^{}=\frac{\rho 1}{\rho }=\frac{p}{12(1p)(1r)},$$
(35)
and, from (30),

$$p^{}=\mathrm{\Delta }=\frac{p(1p)(1r)\varphi }{12(1p)(1r)}.$$
(36)
It is evident that (35) and (36) can only be consistent if $`r1`$. Thus, simulating the urn process with $`p^{}=(\rho 1)/\rho `$ and $`r=1`$ would result in the same value for $`\rho `$ as that obtained from $`p`$ and the original value of $`r`$. However, in the presence of pin removal, the value of $`\mathrm{\Delta }`$ (i.e. the asymptotic value of $`balls/pins`$) would be less than $`p^{}`$, the probability of inserting a new ball into the first urn. This provides a discriminator between the two processes. As a result, by examining the said snapshot we are able to ascertain whether the process is consistent with the urn model in which pins may be discarded, i.e. with $`r<1`$. We apply this analysis in the next section to the degree distribution of inlinks in the web graph.

## 5 A Model for the Evolution of the Web Graph

We now describe a discrete stochastic process for simulating the evolution of the degree distribution of inlinks in the web graph. In this model, balls correspond to web pages and pins correspond to inlinks. At each time step the state of the web graph is a directed graph $`G=(N,E)`$, where $`N`$ is its node set and $`E`$ is its arc set. In this scenario $`F_i(k)`$, $`i1`$, is the number of pages (nodes) in the web graph having $`i`$ inlinks. We note that, although we have chosen $`i`$ to denote the number of inlinks, $`i`$ could alternatively denote the number of outlinks, the number of pages in a web site, or any other reasonable parameter we would like to investigate; see \[LFLW02\] for further details.

Consider the evolution of the web graph with respect to the number of pages having $`i`$ inlinks at the $`k`$th step of the process. Initially $`G`$ contains just a single page. At each step one of three things can happen. First, with probability $`p`$ a new page having one incoming link is added to $`G`$; this is equivalent to placing a new ball in $`urn_1`$. Alternatively, with probability $`1p`$ a page is chosen, the probability of choosing a given page being proportional to $`i`$, the number of inlinks the page currently has; this is equivalent to preferentially choosing a ball from $`urn_i`$. Then, with probability $`r`$ the chosen page receives a new inlink; this is equivalent to transferring the chosen ball from $`urn_i`$ to $`urn_{i+1}`$. Alternatively, with probability $`1r`$ an inlink to the page is removed, and if the page has no inlinks remaining it is removed from the graph; this is equivalent to transferring the ball to $`urn_{i1}`$ when $`i>1`$, and discarding the ball from the system if $`i=1`$.

As a proof of concept, we use the inlinks data from a large crawl of the web performed during May 1999 \[BKM+00\]; the data set we used in this analysis was obtained from Ravi Kumar at IBM. After removing web pages with zero inlinks from the data set, there remain approximately 177 million pages with a total of 1.466 billion inlinks. Using the ranked approach, linear regression on the log-log transformed data gives an exponent of $`2.1052`$, which is consistent with the value of $`2.1`$ reported in \[BKM+00\] and confirmed subsequently in \[AB02\]; see also \[DKM+02\].

As discussed in the introduction, we have some reservations about the use of the ranked approach, as it uses ranks rather than values, thereby ignoring the gaps. In Broder et al.’s May 1999 data set, the first gap appears at 3121, i.e. there were no pages found with 3121 inlinks. Two possible reasons for such a gap are: (i) there were no pages in the web having 3121 inlinks, or (ii) there were one or more such pages but the crawl did not cover these pages. In any case, there is no reason to believe that the web will not contain such a page in the future. Moreover, such gaps, which are inherent in preferential attachment models, change over time due to the growth of the web graph and the stochastic nature of the evolutionary process.

One way to avoid this issue is to use the unranked approach and carry out the regression only on the data values preceding the first gap. Using the first 3120 data values of the inlinks data set, the regression yields an exponent of $`2.1535`$. In Figure 1 we show the two regression lines, which have negative slopes of $`2.1052`$, when using the ranked approach, and $`2.1535`$, when using the unranked approach. (Recall that the exponent of the power law is $`1+\rho `$.)

From the inlinks data we then estimated $`\mathrm{\Delta }`$ and $`\varphi `$, obtaining $`\mathrm{\Delta }balls/pins=0.1205`$ and $`\varphi |urn_1|/pins=0.0553`$. Using these estimates for $`\mathrm{\Delta }`$ and $`\varphi `$, and $`\rho =1.1535`$, we computed $`p`$ and $`r`$ from (33) and (34), respectively, obtaining $`p=0.0915`$ and $`r=0.8280`$. From the analysis in Section 4, it is evident that link deletion is taking place in the web graph, since $`r<1`$. The extent of link deletion indicated as a proportion of insertions and deletions is $`(1p)(1r)=0.156`$. It would be interesting to check this value empirically by looking at sequential snapshots of the web.

To validate our mean-field analysis we conducted 10 simulation runs each for $`k=10^7`$ steps with the above values of $`p`$ and $`r`$. Over the 10 runs, the mean for $`\mathrm{\Delta }`$ was $`0.1197`$ with standard deviation $`1.1\times 10^4`$, and computing $`\varphi `$ as $`|urn_1|/pins`$ gave a mean of $`0.0617`$ with standard deviation $`9.2\times 10^5`$. However, if we compute $`\varphi `$ from (29) the mean is $`0.0589`$ with standard deviation $`4.9\times 10^4`$, which is closer to the empirical value of $`0.0553`$. We suggest that the difference between the two estimates for $`\varphi `$ is mainly due to the slow convergence of $`|urn_1|/pins`$.

Lastly, we computed $`\rho `$ from the simulation results as

$$\frac{pins}{pinskp},$$

which follows from (26) and (28). This gave a mean value for $`\rho `$ of $`1.1534`$, with standard deviation $`5.6\times 10^5`$, compared to the value of $`1.1535`$ computed from the mean-field equation (26).

For regression purposes 10 million simulation steps are insufficient to get close to the asymptotic value of $`\rho `$, so we ran two additional simulations of 1 billion steps to compare the exponents obtained using the ranked and unranked approaches. For efficiency reasons, we modified the simulation by only considering the first $`7000`$ urns, approximating the effect at the upper boundary. The simulation was implemented in Matlab, running on a desktop PC with a 1485Mhz Intel Pentium 4 processor and 500MB of RAM, on a Windows 2000 platform. These 1 billion step simulations, with $`7000`$ urns, each took over 300 hours. (We conducted several more runs of 1 billion steps, varying the numbers of urns, to validate the robustness of the approximation; a single run of 2 billion steps with $`15000`$ urns gave similar results to those we report below.)

For the first (second) run with $`7000`$ urns, we found that the first empty urn was urn number $`2403`$ ($`2395`$) and that overall there were $`5859`$ ($`5914`$) non-empty urns. For the unranked approach, linear regression on the log-log transformed data for the first $`2402`$ ($`2394`$) urns gave an exponent of $`2.1603`$ ($`2.1558`$). On the other hand, for the ranked approach, regression on all the non-empty urns gave an exponent of $`2.1199`$ ($`2.1095`$). In summary, the simulations of our stochastic model are consistent with the inlinks data set, highlighting a small but noticeable difference in the exponent depending on whether the ranked or unranked approach is used. In fact, with our evolutionary model of the web graph, in common with all others, there is no easy way of handling the gaps and, moreover, the concept of gaps has no meaning in the context of an asymptotic mean field analysis.

## 6 Concluding Remarks

We have presented an extension of Simon’s classical stochastic process that allows for pins (which might represent web links) to be discarded, and have shown that asymptotically it still follows a power-law distribution. Given a snapshot of a system, the mean-field equations that we have derived give estimates of the parameters $`p`$ and $`r`$, which can then be input to our stochastic process simulating the evolution of the system. We applied our analysis to the May 1999 web crawl data \[BKM+00\] to detect the extent to which link deletion had taken place. The values of $`p`$ and $`r`$ that we have obtained indicate that approximately 15% of all link operations are deletions. We also ran a number of simulations to validate the mean-field analysis.

An interesting finding that came to light when analysing the data was that there is good evidence that the exponent of the power law for inlinks is in fact close to $`2.15`$ rather than to the widely published value of $`2.1`$. Although the difference between these exponents is small, we consider it to be significant because it suggests a more justifiable way to use linear regression to obtain exponents from power law data.

It would be interesting to study link deletion through historical data, such as that provided by the wayback machine \[Not02\] (http://www.waybackmachine.org), in order to gain more insight into the dynamic aspects of the evolution of the web graph. In particular, it would be desirable to determine whether and how the exponent, and the parameters $`p`$ and $`r`$, are changing over time.

Acknowledgements. The authors would like to thank Ravi Kumar at IBM for providing us with the inlinks data set. We would also like to thank the referees for their constructive comments, which have helped us to improve the presentation of the results.