text:
# Blue horizontal branch stars in metal-rich globular clusters. II. 47 Tuc and NGC 362

## 1 Introduction

The colour-magnitude diagrams (CMD’s) of globular clusters show a great variety of horizontal branch (HB) morphologies. It has long been known that the HB morphology tends to become bluer with decreasing cluster metallicity. However, the existence of pairs of globular clusters with the same metallicity but different HB morphologies shows that there must be at least one “second parameter”. Theoretically, an increase in age is expected to reduce the envelope mass and thus yield a bluer HB morphology, provided that no other parameters are changing. If age is the second parameter, then there are profound consequences for models of the formation of the galactic halo (e.g. Sarajedini et al. sara95 (1995)) and the origin of the ultraviolet upturn in elliptical galaxies (e.g. Park & Lee park97 (1997)). However, much recent work indicates that age may not be the only or even the dominant second parameter controlling the HB morphology (Stetson et al. stet96 (1996), Ferraro et al. 1997a , Catelan cate00 (2000)). Fusi Pecci et al. (fusi93 (1993)) suggested that the cluster density may influence the extension of the blue HB to rather faint magnitudes (“blue tail”) that is seen in the CMD’s of many globular clusters, due to enhanced mass loss triggered by stellar interactions (see also Bailyn bail95 (1995)).

The clusters 47 Tuc and NGC 362 are famous for their red HB morphologies. 47 Tuc is the prototype metal-rich (\[Fe/H\] = $``$0.71; Harris harr96 (1996)) red HB cluster, with only a single known RR Lyr star (V9; Carney et al. cast93 (1993)). NGC 362 \[Fe/H = $``$1.16; Harris harr96 (1996)) is the red member of a famous second-parameter pair of clusters (along with NGC 288). Thus it was a surprise when observations with the Ultraviolet Imaging Telescope (UIT) revealed the presence of 51 hot ($`\mathrm{T}_{\mathrm{eff}}`$ $`>`$ 8500 K) stars in the field of 47 Tuc (O’Connell et al. ocdo97 (1997)) and 84 hot stars within a radius of 8$`.\mathrm{}`$25 from the center of NGC 362 (Dorman et al. dosh97 (1997)). Although both 47 Tuc and NGC 362 are located in projection on the outer halo of the Small Magellanic Cloud (SMC), O’Connell et al. and Dorman et al. were able to show from the distribution with cluster radius that a significant fraction of the detected hot stars should be cluster members. Independently, Kaluzny et al. (kakr97 (1997), 1998a ) used optical photometry to identify eight candidate blue HB stars in 47 Tuc (all of which were detected with UIT), but suggested from the absence of a radial gradient that most of these candidates were actually SMC halo members.

If the hot stars in NGC 362 and 47 Tuc are HB stars and cluster members, then they should have the same age and metallicity as the majority red HB population, and might provide important clues to the other parameters that can yield blue HB stars in globular clusters. However, as noted above, not all the hot stars detected on the UIT images are cluster HB stars. From the distribution of young stars in the outer halo of the SMC mapped by Gardiner & Hatzidimitriou (gaha92 (1992)), we expect significant SMC contamination of the hot stars in NGC 362, whereas 47 Tuc is located in a region mostly devoid of young stars in the northwest SMC halo. Some of the hot stars in NGC 362 and 47 Tuc might also be extreme blue stragglers, analogous to those seen in the cores of NGC 6397 (Burgarella et al. bupa94 (1994)) and M3 (Ferraro et al. 1997b ). The temperatures of these blue stragglers can reach 12,000 K, and recent HST spectroscopy has confirmed that the mass of the two bright blue stragglers in the core of NGC 6397 exceeds twice the turnoff mass (Sepinsky et al. sesa00 (2000)) .

To reach any conclusions about these stars and their evolutionary status we need information about their atmospheric parameters (effective temperature and surface gravity) and also about their cluster membership. We therefore started a programme to obtain spectra of the hot HB star candidates in 47 Tuc and NGC 362. The radial velocities can be used to assess cluster membership, and the derived atmospheric parameters ($`\mathrm{T}_{\mathrm{eff}}`$ and $`\mathrm{log}\mathrm{g}`$) can be used to distinguish HB stars from blue stragglers and background stars.

We describe the observations and their reduction Sects. 2 and 3. The determination of atmospheric parameters is described in Sect. 4 and the results are discussed in Sect. 5.

## 2 Observations

We selected a total of 12 blue stars in 47 Tuc and NGC 362 that were sufficiently uncrowded to allow ground-based follow-up spectroscopy. These stars are listed in Table 1 and marked in Fig. 1 with triangles.

We observed with the ESO Multi Mode Instrument (EMMI) at the New Technology Telescope (NTT) on October 26 and 27, 1997, using only the blue channel as the use of the two-channel mode leads to destructive interference around the Hβ line. We obtained low resolution spectrophotometric data with large slit widths to analyse the flux distribution and medium resolution spectra to measure Balmer line profiles and helium line equivalent widths. Seeing values for these observations varied between 0$`.\mathrm{}`$7 and 1$`.\mathrm{}`$3, but the nights were not photometric. Since no EMMI grism allows observations below 3600 Å (necessary to measure the Balmer jump) we used grating #4 (72 Å/mm) and reduced the dispersion for the low resolution spectra by binning along the dispersion axis by a factor of 4. We used a slit width of 1$`.\mathrm{}`$0 (5$`.\mathrm{}`$0) for the medium (low) resolution spectra and the slit was kept at parallactic angle for all observations. The CCD was a Tek 1024$`\times `$1024 chip with (24$`\mu `$m)2 pixels, a read-out-noise of 5.7e- and a gain of 2.84 e-/count.

For calibration purposes we observed each night 10 bias frames and 10 dome flat-fields with a mean exposure level of about 10,000 counts each. In addition we obtained sky flat fields to correct for the slit illumination. We took wavelength calibration spectra before and after each science spectrum. As flux standard stars we used LTT 7987 and EG 21.

## 3 Data Reduction

We averaged the bias frames of the two nights and used their mean value instead of the whole frames as there was no spatial or temporal variation detectable. To correct the electronic offset we adjusted the mean bias by the difference between the mean overscan value of the science frame and that of the bias frame. The dark currents were determined from several long dark frames and turned out to be negligible (3$`\pm `$2 counts/hr/pixel).

The flat fields were averaged separately for each night, since there was a slight variation in the fringe patterns of the flat fields from one night to the next (below 5%). We averaged the sky flat fields and then condensed them along the dispersion axis to derive the spatial profile of the slit. The dome flat fields were treated the same way, divided by their spatial profile and then multiplied by the spatial profile derived from the sky flat fields. The spectral energy distribution of the flat field lamp was determined by averaging the dome flat fields along the spatial axis. This one-dimensional “flat field spectrum” was then heavily smoothed and was afterwards used to normalize the dome flats along the dispersion axis.

For the wavelength calibration we fitted 3rd order polynomials to the dispersion relations. We rebinned the frames two-dimensionally to constant wavelength steps. Before the sky fit the frames were median filtered along the spatial axis to erase cosmic ray hits in the background. To determine the sky background we had to find regions without any stellar spectra, which were sometimes not close to the place of the object’s spectrum. Nevertheless the flat field correction and wavelength calibration turned out to be good enough that a constant spatial distribution fit the sky light well enough to subtract the sky background at the object’s position with sufficient accuracy. This means in our case that we do not see any absorption lines caused by the predominantly red stars of the clusters or by moon light. The fitted sky background was then subtracted from the unsmoothed frame and the spectra were extracted using Horne’s algorithm (Horne horn86 (1986)) as implemented in MIDAS.

Finally the spectra were corrected for atmospheric extinction using the data of Tüg (tueg77 (1977)). The data for the flux standard stars were taken from Hamuy et al. (hamu92 (1992)). The response curves were fitted by 8th order polynomials for the low and medium resolution spectra. We took special care to correctly fit the response curves for the low resolution data in the region of the Balmer jump, since we planned to use this feature for the determination of the effective temperatures of the stars. The response curves showed a roughly constant shape, but varied in intensity, indicating non-photometric conditions during the observations.

We also used the medium resolution data to derive radial velocities, which are listed in Table 1 (corrected to heliocentric system). The error of the velocities is about 30 km/sec (estimated from the r.m.s. scatter of the velocities derived from individual lines). The normalized and velocity-corrected spectra are plotted in Fig. 2.

## 4 Atmospheric Parameters

### 4.1 $`\mathrm{T}_{\mathrm{eff}}`$ derived from energy distributions

To avoid any spurious results due to contamination from neighbouring stars we searched all stars for red neighbours. For this purpose we used the photometry of 47 Tuc and NGC 362 (Montgomery & Janes, priv. comm.), and extracted for each of our targets all stars within a radius of 20″ (down to the limiting magnitude of $`B`$ $`18\stackrel{\mathrm{m}}{.}5`$). We then estimated the straylight provided by these stars assuming that the seeing at larger distances is best described by a Lorentz profile and scaling the intensities with the V fluxes. We used 1$`.\mathrm{}`$5 as “intrinsic seeing” for all observations (which overestimates the seeing near the meridian) and took into account the elongation caused by atmospheric dispersion. Under these assumptions we get stray light levels of more than 3% in $`B`$ only for MJ3832 (36%) in NGC 362.

We used line blanketed LTE model atmospheres of Kurucz (kuru92 (1992), ATLAS9) for a metallicity \[M/H\] = $``$1.0 to derive effective temperatures from the low resolution spectra. To correct for interstellar reddening we applied the reddening law of Savage & Mathis (sama79 (1979)) and used an EB-V of $`0\stackrel{\mathrm{m}}{.}04`$ for both clusters.

The simultaneous fitting of the spectrophotometric data (Balmer jump and continuum) was generally possible except for MJ38529 in 47 Tuc and MJ6558 and MJ3832 in NGC 362. MJ3832 has a cool star close by (MJ3749, distance 3″, $`B`$ = $`14\stackrel{\mathrm{m}}{.}39`$, $`BV`$ = $`+1\stackrel{\mathrm{m}}{.}43`$) which contaminates the low resolution spectrum that was observed with a 5″ slit (see above). MJ6558 was observed at rather high airmass (1.8) and MJ38529 is probably a binary (see below). However, comparing our measured intensities to those derived from $`B`$ and $`V`$ magnitudes of the stars led to discrepancies of up to $`0\stackrel{\mathrm{m}}{.}3`$. In most of the cases the observed spectra were brighter than the corresponding photometric fluxes. The slope as derived from $`BV`$ was consistent with the slope of the spectrophotometric continuum in most cases except MJ38529 in 47 Tuc and MJ5381 in NGC 362 ($`BV`$ being more positive than indicated by the spectrophotometric continuum). These variations are consistent with non-photometric weather conditions, assuming that the varying atmospheric absorption was basically grey, i.e. independent of wavelength. If the absorption towards the flux standard stars were greater than towards the globular clusters, the cluster stars would be calibrated too bright. We found that the observed Hβ line in the low resolution spectrum of MJ8453 in NGC 362 was significantly shallower than the theoretical one (see also below).

### 4.2 Balmer line profile fits

To derive effective temperatures, surface gravities, and helium abundances we fitted the observed Balmer and helium lines with appropriate stellar model atmospheres. Beforehand we corrected the spectra for radial velocity shifts, derived from the positions of the Balmer lines. To establish the best fit we used the routines developed by Bergeron et al. (besa92 (1992)) and Saffer et al. (saff94 (1994)), which employ a $`\chi ^2`$ test. The $`\sigma `$ necessary for the calculation of $`\chi ^2`$ is estimated from the noise in the continuum regions of the spectra. The fit program normalizes model spectra and observed spectra using the same points for the continuum definition. We computed model atmospheres using ATLAS9 (Kurucz 1991, priv. comm.) and used Lemke’s version111For a description see http://a400.sternwarte.uni-erlangen.de/$``$ai26/linfit/linfor.html of the LINFOR program (developed originally by Holweger, Steffen, and Steenbock at Kiel university) to compute a grid of theoretical spectra, which include the Balmer lines Hα to H22 and He i lines. The grid covered the ranges 8000 K$``$ $`\mathrm{T}_{\mathrm{eff}}`$ $``$20,000 K, 2.5$``$ $`\mathrm{log}\mathrm{g}`$ $``$6.0, and $``$2.0, $``$1.0 in $`\mathrm{log}\frac{\mathrm{n}_{\mathrm{He}}}{\mathrm{n}_\mathrm{H}}`$ at metallicities of $``$1 and $``$0.75.

We fitted the Balmer lines from Hβ to H12 (excluding Hϵ because of the Ca ii H line) and the He i lines $`\lambda \lambda `$ 4026 Å, 4388 Å, 4472 Å, 4713 Å, and 4922 Å (for $`\mathrm{T}_{\mathrm{eff}}`$ $``$ 10,500 K). For cooler stars we fitted only the Balmer lines for a fixed helium abundance of $`\mathrm{log}\frac{\mathrm{n}_{\mathrm{He}}}{\mathrm{n}_\mathrm{H}}`$ = $``$1 (the helium abundance in these cooler stars should be close to solar as they should not be affected by diffusion). In the medium resolution spectrum of MJ8453 the observed Hβ line is significantly less deep than the theoretical one and also Hγ shows some evidence for filling. We therefore excluded those two lines from the fit for this star. Several stars above 10,500 K do not show He i lines. We nevertheless fitted those regions of the spectra, where the He i lines are expected, to obtain an upper limit of the He abundance.

To check for any effects of metallicity on the final results we fitted the spectra with models of \[M/H\] = $``$1.00 and \[M/H\] = $``$0.75. We found that in all cases except MJ3832 in NGC 362 the differences in effective temperature and surface gravity were below 1% resp. 0.05 dex (2% and 0.11 dex for MJ3832). As the differences are below our expected errors we decided to use only the results for \[M/H\] = $``$1 for the further discussion, which lies between the metallicity for 47 Tuc (\[M/H\] = $``$0.71) and that of NGC 362 (\[M/H\] = $``$1.16).

The fit program gives r.m.s. errors derived from $`\mathrm{\Delta }\chi ^2`$ = 2.71 ($`\mathrm{T}_{\mathrm{eff}}`$, $`\mathrm{log}\mathrm{g}`$) resp. 3.53 ($`\mathrm{T}_{\mathrm{eff}}`$, $`\mathrm{log}\mathrm{g}`$, $`\mathrm{log}\frac{\mathrm{n}_{\mathrm{He}}}{\mathrm{n}_\mathrm{H}}`$). However, if $`\chi ^2`$ is not close to 1 (see Table 2) these errors will most likely underestimate the true errors. As the errors causing $`\chi ^2>1`$ are most probably systematic (see below) they are hard to quantify. We decided to get an estimate of their size by assuming that the large $`\chi ^2`$ is solely caused by noise. Increasing the noise parameter $`\sigma `$ until $`\chi ^2`$ = 1 then yields new formal errors for $`\mathrm{\Delta }\chi ^2`$ = 2.71 and 3.53, respectively, which are given in Table 2. These errors still underestimate the internal errors as we oversampled the spectra by a factor of 2.5 when rebinning to constant wavelength steps. Thus only 40% of the wavelength points used for the fit are truly independent. We correct this by applying a factor of $`\sqrt{2.5}=1.6`$ to all error estimates. We believe that the $`\chi ^2`$ values above 1 are due to the low resolution of the data (5.7 Å): At $`\mathrm{T}_{\mathrm{eff}}`$ = 11,000 K and $`\mathrm{log}\mathrm{g}`$ = 4 the FWHM of the theoretical Balmer lines are around 12 Å. Thus the instrumental profile makes up a considerable part of the observed line profile. Therefore any deviations of the instrumental profile from a Gaussian (which is used to convolve the model spectra) will lead to a bad fit. Fortunately we could compare the effects for one hot HB star (B 3253) in NGC 6752, for which we have an NTT spectrum with the same setup as is used here and a spectrum from the ESO 1.52m telescope with a resolution of 2.6 Å. While the reduced $`\chi ^2`$ values for the two spectra differ by a factor of 2 (5.9 vs. 3.1) the resulting effective temperatures and surface gravities are rather similar (13,700 K/3.75 vs. 13,700 K/3.80, see Moehler et al. mosw00 (2000) for details).

### 4.3 Comparison of effective temperatures obtained with different methods

As can be seen from Table 2 the temperatures determined from the spectrophotometric data deviate from those derived from the Balmer lines, which may be an effect of strongly varying atmospheric extinction (as is already suggested by the discrepancies between spectrophotometric and photometric fluxes for the stars, see Sect. 4.1). Differences in the reddening and interstellar extinction law between galactic globular clusters and the SMC may also affect the flux distribution of the stars, which are by default only corrected for the galactic extinction towards the respective cluster. Sasselov et al. (sabe97 (1997)) report that the extinction law for the SMC is essentially the same as that for the Large Magellanic Cloud (LMC). We therefore dereddened the extinction corrected spectra once more, this time using $`\mathrm{E}_{\mathrm{B}\mathrm{V}}`$ = $`0\stackrel{\mathrm{m}}{.}09`$ and the reddening law towards the LMC (Howarth howa83 (1983)) as implemented in MIDAS. The resulting differences between spectra with and without this additional reddening correction are too small to explain the temperature differences between the results from the low resolution data and the line profiles. We thus conclude that the most probable explanation for the discrepancies lies with the non-photometric observing conditions.

### 4.4 Masses

Knowing $`\mathrm{T}_{\mathrm{eff}}`$, $`\mathrm{log}\mathrm{g}`$, and the distances of the stars we can derive the masses:

$`\mathrm{log}\frac{\mathrm{M}}{\mathrm{M}_{}}=const.+\mathrm{log}g+0.4((mM)_VV+V_{th})`$

(Vth denotes the theoretical brightness at the stellar surface as given by Kurucz kuru92 (1992).) The results are listed in Table 2. We determine two masses for each star, one assuming that it is a cluster member and one assuming it belongs to the SMC. We use the following reddening-free distance moduli $`(mM)_0`$: $`13\stackrel{\mathrm{m}}{.}31`$ for 47 Tuc, $`14\stackrel{\mathrm{m}}{.}67`$ for NGC 362 (both from Djorgovski djor93 (1993); Harris harr96 (1996) gives $`13\stackrel{\mathrm{m}}{.}25`$ resp. $`14\stackrel{\mathrm{m}}{.}65`$ in the June 22, 1999 version of his compilation), and $`18\stackrel{\mathrm{m}}{.}91`$ for the SMC (Sasselov et al. sabe97 (1997)). Using the Hipparcos based distance moduli to 47 Tuc and NGC 362 of $`13\stackrel{\mathrm{m}}{.}51`$ and $`14\stackrel{\mathrm{m}}{.}88`$ (Gratton et al. grfu97 (1997), Reid reid98 (1998)) would increase the masses by about 20% in both cases.

In addition to the error in $`\mathrm{log}\mathrm{g}`$ given in Table 2 errors in the absolute magnitude and the theoretical brightness at the stellar surface also enter the final error in $`\mathrm{log}\mathrm{M}`$. We assume an uncertainty in the absolute brightness (combining errors in the photometric data, uncertainties in the distance moduli and reddenings) of $`0\stackrel{\mathrm{m}}{.}3`$, and a mean error in $`\mathrm{log}\mathrm{g}`$ of 0.2 dex. The error in the theoretical $`V`$ brightness is dominated by the errors in $`\mathrm{T}_{\mathrm{eff}}`$. We estimate the total error in $`\mathrm{log}\mathrm{M}`$ to be about 0.24 dex, corresponding to an error of $`{}_{42\%}{}^{}{}_{}{}^{+73\%}`$ in mass (for a detailed discussion of error sources see Moehler et al. mohe97 (1997)).

## 5 Discussion

47 Tuc: As can be seen from Table 2 three of the four stars in 47 Tuc are probably members of the cluster, but MJ8279 according to its radial velocity belongs to the SMC. The comparison of the spectroscopically derived mass (3.8$`\mathrm{M}_{}`$) with the mass obtained from the evolutionary tracks (5$`\mathrm{M}_{}`$) further supports its SMC membership. Kaluzny et al. (kakr97 (1997),1998b ) reported observations of six stars in the field of 47 Tuc with similar magnitudes and colors to MJ 8279, that is $`18<`$ V $`<18.7`$, and $`BV0.1`$. The evidence of MJ 8279 supports their conclusion (based on a lack of a radial gradient) that most of these stars are probably SMC members.

Of the remaining three stars MJ38529 is very probably a binary: Its UV-visual colour suggests an effective temperature of 11,900 K, but $`BV`$ for this star is much too red222Kaluzny et al. (1998a ) give a relatively red $`VI`$ colour of $`+0.14`$ for MJ38529, which provides further evidence for the presence of a cool companion. for such a temperature (see also its position in Fig. 1 below the HB). The optical low resolution spectrum suggests effective temperatures of 12,500 K (Balmer jump only) resp. 8000 K (continuum only). The superposition of two different stellar spectra (with different Balmer lines) may explain the strange values obtained for $`\mathrm{log}\mathrm{g}`$ from the Balmer line profiles. If we assume that the temperature obtained from UV data is that of the hot component MJ38529hot, it would lie between MJ280 and MJ33410 in $`\mathrm{T}_{\mathrm{eff}}`$. If MJ38529hot were a ZAHB star, its $`V`$ magnitude should also lie between those two stars – in contradiction to the fact that the complete binary system is already fainter in $`V`$ than MJ280 and MJ33410. MJ38529hot could thus also be an extreme blue straggler, which would be fainter than a ZAHB star of the same $`\mathrm{T}_{\mathrm{eff}}`$. On the other hand, the UV-based temperature is a lower limit for the hot component, as optical data are also involved in its determination (UV-visual colour). If MJ38529hot were considerably hotter than MJ280 it would be much fainter visually and could therefore be a hot ZAHB star accompanied by a red subgiant. Thus the contamination by the cool star prevents any definite conclusions.

MJ280 is hot enough to be affected by diffusion processes which may result in an enrichment of heavy elements in its atmosphere (cf. Moehler et al. mosw00 (2000)). Fitting it with enriched model atmospheres as described by Moehler et al. (mosw00 (2000), \[Fe/H\] = $`+`$0.5) indeed yields $`\mathrm{T}_{\mathrm{eff}}`$ = 14,200 K, $`\mathrm{log}\mathrm{g}`$ = 4.29 and $`\mathrm{log}\frac{\mathrm{n}_{\mathrm{He}}}{\mathrm{n}_\mathrm{H}}`$ = $``$1.63, moving the star closer to the ZAHB.

MJ33410 finally is a normal cool HBB star. The masses of MJ280 and MJ33410 agree well with canonical predictions.

NGC 362: According to Table 1 three of the eight targets in NGC 362 are radial velocity members. MJ6558 is not a proper motion member, but the error of its proper motion is rather large compared to the measurements of the other stars. MJ6558 and MJ8241 lie close to the canonical ZAHB, whereas MJ3832 has evolved towards lower $`\mathrm{log}\mathrm{g}`$. The mean mass of all three stars is 0.50$`{}_{0.23}{}^{}{}_{}{}^{+0.44}`$$`\mathrm{M}_{}`$ (r.m.s. errors only) and thus also in good agreement with canonical HB theory. Note that MJ6558 was one of the three candidate blue HB stars (= H1382) in NGC 362 identified by Harris et al. (harr82 (1982)); the other two Harris candidates were not detected with UIT and radial velocity studies indicate that they are unlikely to be cluster members (R. Rood, priv. comm.). None of the radial velocity members are good candidates for being extreme blue stragglers, in the sense of having both a mass and $`\mathrm{T}_{\mathrm{eff}}`$ consistent with being on an extension of the NGC 362 main sequence.

The three radial velocity member stars are all located within 2$`.\mathrm{}`$5 of the cluster center, while the remaining five stars are all more than 3$`.\mathrm{}`$5 from the center. Three of these five stars may belong to the SMC according to their radial velocities: MJ65, MJ94, and MJ5381. They lie, however, farther from the SMC main sequence in Fig. 1 than the remaining two stars MJ2341 and MJ8453. Comparing the spectroscopically determined masses to those derived from evolutionary sequences would put MJ65 (9.5$`\mathrm{M}_{}`$ vs. 7$`\mathrm{M}_{}`$), MJ94 (7.6$`\mathrm{M}_{}`$ vs. 4$`\mathrm{M}_{}`$), and MJ2341 (5.9$`\mathrm{M}_{}`$ vs. 5$`\mathrm{M}_{}`$) to the SMC. The radial velocity of MJ2341 differs by 74 km s-1 from that of the SMC, and so it might be a high-velocity halo star of the SMC (see, e.g., Keenan keen92 (1992), keen97 (1997) for a discussion of main-sequence B stars in the halos of galaxies). While two of the three stars lie away from the main sequence one should keep in mind that our selection is heavily biased towards bright (= evolved) hot SMC stars. Note that MJ65 and MJ94 appear to be normal HB stars in the UV CMD of Dorman et al. (dosh97 (1997)), which emphasizes the need for spectroscopic observations to verify the evolutionary status of stars in fields with high contamination.

The evolutionary status of MJ5381 is unclear: In the $`m_{162},m_{162}V`$ diagram of Dorman et al. (dosh97 (1997)) it lies rather isolated – being fainter than the theoretical ZAHB for NGC 362, but much cooler than the SMC main sequence.

The filled Hβ line of MJ8453 reminds us of the Be stars analysed by Mazzali et al. (male96 (1996)), but it is considerably cooler than those objects. A spectrum of the Hα region of this star would be helpful to decide whether it indeed shows a Be-like spectrum.

## 6 Conclusions

There has been increasing interest in understanding the origin of the hot stars that can appear in old, metal-rich systems. The metal-rich open clusters NGC 6791 (\[Fe/H\] = +0.4) and NGC 188 (\[Fe/H\] = 0.0) contain hot HB populations sufficiently large that the UV-turnup in their integrated cluster spectra would be as strong as that observed in elliptical galaxies (Landsman et al. labo98 (1998)). At these high metallicities, hot HB stars can be a natural product of single-star HB models (e.g. Dorman et al. dorm95 (1993)), although a significant fraction of the hot stars in the two clusters are binaries (Green et al. grli97 (1997)). The metal-rich bulge globular clusters, NGC 6388 (\[Fe/H\] = $``$0.60, Harris harr96 (1996)) and NGC 6441 (\[Fe/H = $``$0.53, Harris harr96 (1996)) contain a sizable population of blue HB stars. A possible scenario for the blue HB in these clusters is a merger of two systems with different ages, with the blue HB arising from the minority old population (Yoon et al. yole00 (2000)). However, this scenario does not explain several other peculiar aspects of the blue HB population in NGC 6388 and NGC 6441, including the slope of the HB, and the gravities which are higher than predicted by canonical models (Moehler et al. mosw99 (1999)).

Here we have shown that at least some hot HB stars exist in the classical red HB clusters 47 Tuc and NGC 362, and – unlike the case of NGC 6388 and NGC 6441 – they perfectly fit the expectations of canonical HB evolution. The three blue members in 47 Tuc all have 10,000 K $`<`$ $`\mathrm{T}_{\mathrm{eff}}`$$`<`$ 15,000 K, and thus are well separated from the rest of the HB population, which is (except for the single RR Lyr V9) entirely redward of the instability strip. (A fourth probable hot star member of 47 Tuc is UIT-14, which is only $`1.7^{}`$ from the cluster center, and for which the IUE spectrum obtained by O’Connell et al. ocdo97 (1997) indicates $`\mathrm{T}_{\mathrm{eff}}`$$``$ 50,000 K.) The small number of hot HB stars in 47 Tuc, and their high temperatures, point to a scenario in which they have a different physical origin (such as binary interactions) than the dominant red HB population. It is suggestive that one of the hot stars in 47 Tuc (MJ38529) appears to be a photometric binary, whose hot component could, however, also be an extreme blue straggler.

Interestingly, whereas NGC 362 has 14 bright blue stars within 14″ radius of its center (Dorman et al. dosh97 (1997)), the inner 1′ of 47 Tuc probably does not contain any hot HB stars, as shown by observations with IUE (Rich et al. rimi93 (1993)), UIT (O’Connell et al. ocdo97 (1997)), and WFPC2 (Rich et al. riso97 (1997)). Rose & Deng (rode99 (1999)) found that only about 7% of the mid-UV light in the core of 47 Tuc comes from stars hotter than about 7,500 K (most of which are probably blue stragglers).

The three hot HB stars in NGC 362 are somewhat cooler than those in 47 Tuc, and the dominant HB in NGC 362 is not as red as that of 47 Tuc. Thus it is plausible that the blue HB stars in NGC 362 arise from a small percentage of red giants with unusually high mass loss, and no special mechanism is required for their production. Unfortunately, the SMC contamination is larger in NGC 362 than in 47 Tuc, and five out of our eight targets turned out to be nonmembers. It would be interesting to study the stellar parameters of the hot stars in the core region, which we could not observe due to high crowding, but where the relative SMC contamination should be much lower.

###### Acknowledgements.

We want to thank the staff of the ESO La Silla observatory for their support during our observations and an anonymous referee for valuable comments. SM acknowledges support by the DLR (grant 50 OR 96029-ZA).


format: text

text:
# 1 Introduction

## 1 Introduction

Among the hydrogen-bonded ferroelectrics, the compound KH2PO4 (KDP) has been the most intensively studied over the second half of the last century. However, a detailed knowledge about its ferroelectric phase transition and some related phenomena is still lacking. In particular, its transition temperature (Tc $``$ 122K), shows an increase of about 107K upon deuteration. The origin of this huge isotopic effect is still controversial. In the sixties, Blinc and Zecks proposed a quantum tunnelling model to explain it, which focuses purely on mass-dependent effects. Improvements of this model include coupling to vibrations of the K-PO4 complex. Experimental studies over the past fifteen years, especially high-resolution neutron difraction measurements, provided increasing evidence that the geometrical modification of the hydrogen bonds upon deuteration, the well-known Ubbelohde effect, is intimately connected with the mechanism of the phase transition. McMahon et al. have shown that the distance $`\delta `$ between the two collective equilibrium positions of the protons is remarkably correlated with Tc. They proposed that $`\delta `$ might be the only relevant parameter for Tc and, in addition, no concluding evidence of tunnelling was established. It was then suggested the necessity of a revised theory for this class of transitions. Recently, a phenomenological model was proposed, which couples, via geometrical details of the H-bond, the proton coordinate with a nonlinear dynamics of the K-PO4 system, where the spontaneous polarization actually develops. This model, which includes both tunnelling and geometrical isotope effects, also allows to explain the Tc behaviour upon deuteration.

Summarizing, it appears well-established that, in this class of systems, proton ordering and geometrical and electronic properties of the host cage are connected in a non trivial way. So far, to the best of our knowledge, no ab initio calculations have been reported for KDP. This type of approach has the advantage of allowing a confident parameter-free analysis of the microscopic electronic and structural changes upon proton ordering that drive the ferroelectric transition. To this end, we have performed Density Functional (DFT) electronic structure calculations  for KDP, within a gradient corrected approach to exchange and correlation, both in the tetragonal paraelectric phase (forcing the protons to be in the center of the H-bonds) and in the orthorhombic ferroelectric phase, where the protons order cooperatively at off-center positions in the H-bond. We thus analyzed in detail the structural parameters and the microscopic electronic density redistributions associated to the proton ordering phenomenon, which eventually provoke the ferroelectric instability.

This work is organized as follows: in Section 2 we describe the methods and the details of the calculations. Section 3 is devoted to the analysis of the geometry results. In Section 4, we study and discuss the polarization mechanism produced by the off-centering of the hydrogen atoms along the H-bonds, and its relation to the ferroelectric instability. Finally, in Section 5 we elaborate our conclusions.

## 2 Ab initio methods

We have performed two different types of ab initio calculations. First, we have used a recently developed, very efficient localized basis (LB) set approach (SIESTA). This is a fully self-consistent, pseudopotential DFT method that employs a basis set of pseudoatomic orbitals (PAO). We have chosen a double-zeta basis set with polarization functions (DZP). The exchange-correlation energy terms were computed using the Perdew-Burke-Ernzerhof gradient-corrected functional. Angular-dependent, norm-conserving Troullier-Martins pseudopotentials  were employed to represent the interaction between ionic cores and valence electrons. We have also included nonlinear core corrections (NLCC) for a proper description of the K ion, due to the large overlap of core and valence charge densities. We used a 125 Ry energy cutoff for the grid computation of numerical integrals , and a value of $`E_c=50`$ meV for the orbital confinement energy, which gives total energies and geometries of sufficient accuracy for this system. Due to the approximate character of these confined basis orbitals, we have checked the overall accuracy by performing standard pseudopotential plane wave (PW) calculations. The parallel version of the Car-Parrinello SISSA FPMD code , and the same set of pseudopotentials, exchange-correlation functional, and NLCC were used. The energy cutoff for the plane wave expansion was set to 150 Ry. In both approaches, the electronic Brillouin zone sampling reduced to the $`\mathrm{\Gamma }`$-point, which is a reasonable approximation for such a wide band gap insulator. We have, however, checked that this sampling is sufficient for the large supercells considered.

In the paraelectric phase, the KDP structure (space group I$`\overline{4}`$2d) has a body-centered tetragonal primitive lattice. For our calculations, we used instead the conventional face-centered tetragonal cell (space group F$`\overline{4}`$d2), which includes 64 atoms (8 formula units). With this choice, the transformation to the orthorhombic structure of the ferroelectric phase is simply described by a change in the ratio of the basal plane lattice parameters, in addition to changes in the internal structure parameters.

## 3 Geometry results

In order to obtain the equilibrium configurations relevant to ferroelectricity in KDP, we performed a series of computational experiments, which are described below.

First we considered a static view of the paraelectric phase by fixing a tetragonal cell with parameters at their experimental values at 293 K  and relaxed the atomic positions under the constraint that the H atoms remain in the middle of the H-bonds. In this way we found the transition state (TS) which connects the two oppositely polarized realizations of the ferroelectric phase. In a second time we still kept fixed the tetragonal lattice, and allowed the hydrogen atoms to relax from their centered positions. When the unstable TS state was slightly perturbed, two H atoms approach each PO4 tetrahedron. The two PO4 oxygens which were approached by the hydrogens (in the following named O2) lie on one side of the (xy)-plane containing the P ion. The other two oxygens (O1) lie on the other side. The geometry relaxation is carried out also for all the other atoms in the unit cell. As a consequence of proton ordering, also the K ions move away from the centered position between P planes, thus leading to a ferroelectric phase with tetragonal lattice, which we denote FPt. Finally, we fully optimized atomic positions and cell parameters. This leads to an orthorhombic ground state struture GS, which corresponds with the actual experimental ferroelectric structure.

The resulting structural parameters are shown in Table 1, together with experimental data for the paraelectric tetragonal and the ferroelectric orthorhombic phases. An excellent overall agreement is observed for the TS and FPt configurations obtained with the PW method, as compared to experimental data. The good quality of the LB method can be assessed by comparison with the PW results in columns 4 and 5. In the FPt calculation, the off-center ordering of the hydrogen atoms along the H-bonds leads to inequivalent O1 and O2 oxygens. Consequently, (a) there is a distortion of the PO4 tetrahedra (see the different P-O1 and P-O2 distances on Table 1) and (b) the K ions abandon the centered positions between the phosphate planes (see $`\mathrm{\Delta }z_K`$ in Table 1), thus inducing a spontaneous polarization which leads to ferroelectricity. These effects will be explained in the next Section. We also verified that, when the hydrogen atoms are constrained to remain in the middle of the H-bonds, then the centered positions are stable for the K ions. This provides a strong evidence that the origin of ferroelectricity is in the off-centering of the protons in the H-bonds.

It is worth mentioning that, in this fixed tetragonal FPt calculation, a substantial anisotropic stress appears in the $`xy`$-plane. Note that the hydrogen off-centering in the H-bond network, which lies nearly on the $`xy`$-plane, leads to an asymmetry of the $`x`$ and $`y`$ directions. Further relaxation, obtained by optimizing also the cell parameters in order to minimize the stress with respect to ambient conditions, releases the tetragonal symmetry in favor of the orthorhombic structure shown in column 3 of Table 1.

## 4 Polarization mechanism for the ferroelectric instability

In order to understand the origin of the ferroelectric instability, we have focused our analysis in the electronic charge redistributions that happen in going from the centered (TS) to the off-centered hydrogen case (FPt). To this purpose we compute the orbital Mulliken populations, which are shown on Table 2 for the different atoms in both cases. It is known that Mulliken populations depend strongly on the choice of the basis set. Differences, however, are much less sensitive, and thus meaningful. In going from TS to FPt, it can be clearly observed an increase of the charge localized around O1, the main contribution ($`70`$%) being provided by the O2 charge decrease. The analysis of bond and self atom-orbitals overlap populations, which are shown in Table 3, also helps to discriminate the origin of this charge redistribution.

The trends observed in Tables 2 and 3 are confirmed by the charge density difference $`\mathrm{\Delta }\rho (𝐫)=\rho _{FP_t}(𝐫)\rho _{TS}(𝐫)`$, which is plotted in the planes determined by the atoms P-O2-H (Fig. 1a) and P-O1$`\mathrm{}`$H (Fig. 1b). Analysing both, Table 3 and Fig. 1a, we observe a significant enhancement of the O1-O1 overlap population, accompanied by a smaller increment of the same sign in the O1-P orbitals. This happens at the expenses of the population of the overlap orbitals in the O1$`\mathrm{}`$H bond. This behaviour around the O1 atom indicates that, as the hydrogen atoms move away from O1, the strongly covalent O1-H bond weakens, its charge being redistributed. A substantial amount of this charge localizes around O1, while part of it localizes in the O1-P orbitals, consistently with the shortening of the O1-P bonds reported in Table 1. The contrary occurs in the vicinity of the O2 atom: as H moves towards O2 along the O1$`\mathrm{}`$H-O2 line, the O2-H distance decreases, thus increasing the degree of covalency and producing a charge localization along the corresponding O2-H bond. This is indicated by the bond overlap population in Table 3, and also by the contours in Fig. 1b. Finally, the O2-O2 and O2-P orbitals loose a significant amount of charge. We name this process of charge redistribution, produced by the collective hydrogen off-centering along the H-bond, the polarization mechanism (PM). The PM can be described in a simplified way as a charge accumulation and localization around O1 at the expenses of a charge depletion and delocalization in the vicinity of O2.

This charge accumulation around the O1 atom can also be viewed as a negative charge defect, which electrostatically attracts the positive K+ ion towards the PO4 cage, along the $`z`$-axis. This electrostatic force is counterbalanced by repulsion forces from orbital overlaps, which stabilize a new equilibrium position of the K+ ion, as can also be seen in Table 1. Our calculation then supports the idea that proton off-centering and ferroelectricity are very correlated phenomena.

Cell relaxation – from FPt to GS – also affects the internal structure (Table 1). As H approaches O2, $`\delta `$ increases, the O1$`\mathrm{}`$H bond weakens and the H-O2 bond becomes stronger. As a consequence, a larger amount of charge localizes around O1, and delocalizes away from O2. This behavior is related to the PM and will be useful to establish a connection between $`\delta `$ and the charge difference $`\mathrm{\Delta }q_O=q(O1)q(O2)`$. Moreover, we observe a striking correlation between $`\mathrm{\Delta }q_O`$, the distortion of the tetrahedra, $`\mathrm{\Delta }`$(P-O) $`d(PO2)d(PO1)`$, and the displacement of the K atoms, $`\mathrm{\Delta }z_K`$. These trends are shown in the theoretical curves in Fig. 2 for two different values of $`\delta `$, obtained from the cell-unrelaxed FPt and the GS calculations.

The crucial experimental observation that the effect of isotopic substitution can be almost exactly reproduced by applying pressure , supports the idea that the main effect of deuteration is to modify the internal geometry. In this sense, our calculations for FPt and GS, corresponding to two different pressures, can be used to analyse isotopic effects. It is known from experiment that H(D) ordering, $`\mathrm{\Delta }z_K`$, $`\mathrm{\Delta }`$(P-O), and spontaneous polarization Ps, show very similar trends just below Tc, in the ferroelectric phase. Moreover, these trends are quite well explained by the present PM (left curves in Fig. 2). Therefore, it is natural to attempt to establish a connection between $`\mathrm{\Delta }q_O`$ and structural data measured in KDP and in its deuterated analog, DKDP. This is done in Fig. 2, by including the experimental values of $`\mathrm{\Delta }z_K`$, $`\mathrm{\Delta }`$(P-O) and Ps for KDP and deuterated DKDP. The comparison between theoretical and experimental curves is very favorable, especially for what concerns the slopes of the curves.

In connection with the proposed PM, the above indicates that the enhancement of the saturated polarization due to deuteration, experimentally observed below Tc,  can be explained merely by the increase in the charge unbalance between O1 and O2 ($`\mathrm{\Delta }q_O`$), which in turn is induced by an increase in $`\delta `$. This picture is in line with the ideas developed in Ref., where the experimental behaviour of Ps and the inverse dielectric constant were rationalized in terms of a proton-dipole interaction model. In fact, at temperatures well below Tc, the protons are strongly bound to the neighboring PO4 dipoles, in a stable minimum of an asymmetric potential. In this asymmetric configuration, the protons localize into the deepest well because the energy needed to be promoted to the excited state is too high, thus ruling out the possibility of tunnelling. In this case, deuteration can only induce the indirect effect of enhancing $`\delta `$ , i.e. the so-called geometric effect. Preliminary ab initio calculations indicate, even more strikingly, that there is actually no double minimum structure in the potential energy surface for single protons. Appart from supporting this scenario, the current ab initio-derived PM establishes the importance of the short-range proton – PO4-dipole interaction for the transition, as was also speculated in .

## 5 Conclusions

We presented the results of ab initio DFT calculations for the H-bonded ferroelectric material KDP, which reproduce very well the structure of the tetragonal and orthorhombic (ferroelectric) phases in KDP. On this basis we identified the microscopic polarization mechanism that leads to the ferroelectric instability. This arises from an electronic charge redistribution between the oxygen atoms involved in the H-bonds, occurring when the protons move away from the center of the bonds. We indirectly showed that the interplay between proton ordering and ferroelectric polarization qualitatively explains the enhancement of the saturated Ps upon deuteration. In addition we provided, for the first time, fully ab initio evidence supporting the idea that isotope effects reduce, at least at low temperatures in the ferroelectric phase, to the indirect geometric effect of enhancing $`\delta `$, as speculated in recent works .

## Acknowledgments

We thank E. Tosatti for his invaluable collaboration in the initial stages of this work. JK also acknowledges very useful discussions with R.J. Nelmes and M.I. McMahon. We also thank P. Ordejón for his assistance in the use of the SIESTA package and C. Cavazzoni for providing us the FPMD SISSA code. RLM and SK acknowledge support from ICTP (Italy) and CONICET (Argentina).


text:
# Iterated fully coordinated percolation on a square lattice

## I INTRODUCTION

Iterated fully coordinated percolation (referred to as IFC percolation from this point on) is an extension to the ordinary percolation and the fully coordinated percolation (referred to as FC percolation hereon) models. These are well-known problems of geometrical phase transition where global connectivity is singular as a function of local connectivity (see, for a review, Stauffer and Aharony ). In ordinary site percolation on a square lattice, sites are randomly occupied with an independent occupation probability $`p`$. Two nearest neighbor sites are considered connected if they are both occupied. A collection of connected sites, i.e., every site in the group is connected to every other site either directly or through the paths of nearest neighbors, is termed a cluster. When $`p`$ is increased until it reaches a critical value $`p_c`$, an infinitely spanning cluster will form for the first time. The square lattice shown in Fig. 1(a) is an example of ordinary site percolation on a square lattice. Occupied sites are shown as circles while unoccupied sites are not shown.

To perform IFC percolation, first perform ordinary percolation. After this is done, determine which of the occupied sites have all four of their nearest neighbors also occupied. Occupied sites having this property are called fully coordinated sites and will remain occupied while the other sites will be removed. This process of winnowing occupied sites was described in  and we called the resulting model the fully coordinated percolation model (FC percolation). The square lattice shown in Fig. 1(b) is an example of fully coordinated percolation constructed from the ordinary percolation sites shown in Fig. 1(a). The gray-filled circles are fully coordinated sites. The final step in the process is to again choose fully coordinated sites but this time from the resulting system constructed from FC percolation. As in the FC percolation case, we also remove the sites that are not fully coordinated. The solid circles in the square lattice shown in Fig. 1(c) are examples of IFC percolation sites constructed from Figs. 1(a) and 1(b).

Notice that this problem is different from bootstrap percolation (see, for a review, Adler ) where sites of less connectivity are iteratively removed until no more sites have insufficient connectivities. This is in contrast to our problem where we only perform the culling of non-fully coordinated sites twice.

In our earlier work , we found that the sizes of the fully-connected subregions (with no internal holes) in an incipient infinite cluster in FC percolation are much larger than those found in ordinary percolation. In this study, we find this same increase in the sizes of the fully-connected subregions to continue when the full coordination requirement is iterated one more time. The regions of high local connectivity are generally interesting because they would seriously affect the vibrational modes in an equivalent disordered network of springs. In particular, it was found previously that they enter into the scaling of the normal mode spectrum .

In the same previous work , we also concluded that the dynamic critical exponents of FC percolation were close to but different from those of ordinary percolation, drawing this conclusion from numerical results on the FC percolation clusters which were grown from a seed site until prescribed sizes were reached. In this work, on the other hand, we obtain numerical results on the corresponding dynamical critical exponents of FC and IFC percolation from the clusters generated statically on square grids of prescribed sizes. Surprisingly, our new results indicate that the numerical values of these exponents are very similar for ordinary, FC, and IFC percolation, and in fact most likely to be the same. Since the actual numerical estimates for FC cases are compatible within the error estimates of both works, this might simply mean that we drew an overzealous conclusion in our previous work. However, we cannot rule out the possibility that the ensembles of clusters used in  and the current work may contain subtle differences which in fact alter dynamic universality classes.

The rest of this paper is divided into four sections. In Sec. II we present numerical results and calculations of the static and dynamic critical exponents of IFC percolation. The geometry of clusters is studied in Sec. III with a particular focus on the sizes and distributions of the highly connected subregions. The difference between the methods of generating the cluster ensemble in  and the current work is discussed in Sec. IV. The final section gives a summary of our results.

## II STATIC AND DYNAMIC CRITICAL BEHAVIOR

We have determined the static and dynamic critical exponents of IFC percolation by simulations on a square grid. The culling of non-fully coordinated sites occurs in two stages. After sites in the grid are randomly occupied with an independent probability $`p`$, all the fully coordinated sites are marked. After this marking process is done, all the occupied but unmarked sites are removed. From the resulting system of remaining sites, we again determine which sites are fully coordinated and remove those which are not. Figs. 1 depict this process. To determine the connectivities among the sites, we have used the standard Hoshen-Kopelman algorithm .

The critical occupation probability, $`p_c`$, is determined from the plot of susceptibility, $`\chi `$, against $`p`$. This plot is shown in Fig. 2. Notice that because of finite size effects, instead of a diverging $`\chi `$ at $`p_c`$, we see clearly defined peaks at $`p`$ values near the true $`p_c`$. The lattice sizes, $`L^2`$, we used in our simulations were $`L=256`$, $`512`$, $`1024,`$ and $`2048`$, and each data point in Fig. 2 is an average over $`1000`$ realizations. According to the plot, the susceptibility peaks are pegged at $`p_{peak}=0.946`$ for the larger lattices. Since the larger the lattices the more our simulations mimic the asymptotic limit, we can conservatively assign the value for the critical occupation probability to be at $`p_c=0.946\pm 0.001`$.

The ratio between critical exponents $`\gamma `$ and $`\nu `$ can be determined from the scaling of susceptibility peak values against the corresponding lattice sizes, i.e., $`\chi (p_{peak},L)L^{\gamma /\nu }`$. Using this relationship we find $`\gamma /\nu `$ for IFC percolation to be $`1.805\pm 0.018`$. This value agrees nicely with the ordinary percolation value which is known exactly to be $`43/241.7917`$. The ratio $`\gamma /\nu `$ can also be used to indirectly determine the fractal dimension, $`d_f`$, from the known scaling relationship, $`d_f=1+\gamma /2\nu `$. The result for IFC percolation is $`d_f=1.903\pm 0.009`$ while the known exact value for ordinary percolation is $`91/481.8958`$. The critical exponent $`\gamma `$ can also be determined independently from the scaling of susceptibility against the occupation probability $`p`$, which goes as $`\chi (p,L=\mathrm{})|pp_c|^\gamma `$, where the $`L=\mathrm{}`$ indicates the scaling to be ideally satisfied in the asymptotic limit of $`L\mathrm{}`$. As in the case of FC percolation, the susceptibility peaks are very close to $`p=1.0`$. This provides data to the right of the peaks to reside in a small interval. We therefore only use data to the left of the peaks to determine the scaling implied. The value we found for $`\gamma `$ in IFC percolation using $`p_c=0.946`$ is $`\gamma =2.398\pm 0.039`$. This value also agrees very well with the value in ordinary percolation, known to be exactly $`43/182.3889`$.

The fractal dimension, $`d_f`$, of IFC percolation can also be determined directly. At $`p_c`$, the size of the largest cluster, $`S`$, scales with $`L`$ as $`SL^{d_f}`$, where $`d_f=1.857\pm 0.007`$. This value is close to the known exact value for ordinary percolation as well as to the value obtainable from our measured $`\gamma /\nu `$.

Listed in Table I are the values for the ratio $`\gamma /\nu `$ and for $`d_f`$ we found from our simulations of ordinary, FC, and IFC percolation. As a reference, the previously determined values of the corresponding critical exponents for ordinary percolation are listed at the bottom row of the table.

The values for the spectral dimension, $`d_s`$, and the walk dimension, $`d_w`$, can be determined from the dynamic critical behavior of IFC percolation. We first construct the transition probability matrix W. The elements, $`W_{ij}`$, of this matrix are the hopping probabilities per step of a random walker to hop from site $`j`$ to site $`i`$. The diagonal elements $`W_{ii}`$ are the probabilities for the walker to stay at site $`i`$ at a time step. In the myopic ant rule, all of the $`W_{ii}`$ are zero. In contrast, in the blind ant rule, the $`W_{ii}`$ are nonzero, their values dependent on how many occupied neighboring sites the walker can hop into (see, for a review,). We are using the blind ant rule for this work.

Once the W matrix is determined, the value of $`d_w`$ can be estimated directly from its largest non-trivial eigenvalue, $`\lambda _1`$. This eigenvalue satisfies the finite size scaling law

$$|ln\lambda _1|1\lambda _1S^{2/d_s}.$$
(1)
Making use of the Alexander-Orbach scaling law , $`d_s=2d_f/d_w`$, and $`SL^{d_f}`$, we obtain

$$1\lambda _1L^{d_w}.$$
(2)
Eq. (2) gives a rather direct way of estimating $`d_w`$ from finite size scaling, in a manner that is considerably simpler compared to methods wherein a detailed knowledge of the eigenvectors is required (see, e.g., ).

Shown in Fig. 3 is our data plotted in a way that is consistent with Eq. (2). For each percolation type, we used $`13`$ different edge lengths of the square lattices that range from $`L=94`$ to $`L=198`$. Furthermore, for each data point, about $`1000`$ random realizations of the W matrix had been generated by Monte Carlo simulations. To numerically determine $`\lambda _1`$ from each W matrix, we used the Arnoldi-Saad method .

Listed in the second column of Table II are the $`d_w`$ values we found following Eq. (2) for ordinary, FC and IFC percolation. Notice that there is an excellent agreement among the values for all three cases and that they agree with the previously known value for ordinary percolation.

In the third column of Table II are listed the $`d_s`$ values calculated using the Alexander-Orbach scaling law. Although all three estimates from this work agree with the previously known value for ordinary percolation within the aggregate error estimates of our work and those of the previous work, there is a slight internal discrepancy between our own estimate of $`d_s`$ for ordinary percolation and that for IFC percolation. Since $`d_s`$ is derived from $`d_f`$ and $`d_w`$, the slight discrepancy may be attributed to the $`d_f`$ value we found for IFC percolation by direct simulation (refer to the third column of Table I). We suspect an unaccounted systematic error occuring when we perform each iteration process. In particular, as we will discuss in more details later, the length scales associated with fully-coordinated subregions appear to increase at each iteration. This suggests that we may have to increase the overall lattice length scale correspondingly in order to maintain a good estimate of the fractal dimension by direct simulation. Since this was not done in this work, a systematic error may have crept in this way.

## III GEOMETRY OF CLUSTERS

Previously , it was found that the number of interior sites in the largest cluster in FC percolation is much greater than that in ordinary percolation at their respective thresholds. In this work we find that this trend continues further for IFC percolation. As an illustration, consider Figs. 4(a) and 4(b). Shown in Fig. 4(a) is the largest cluster from a simulation of ordinary percolation on a $`200`$x$`200`$ square lattice at its percolation threshold. The darker dotted sites are fully-coordinated (FC) sites while the lighter gray dotted sites are partially-coordinated (PC) sites. FC sites are those all of whose neighbors are also occupied (thus these are similar to, though not identical with, the interior sites considered in ), while PC sites have at least one nearest neighbor not occupied. The same shading tone is used in Fig. 4(b) except that what is shown is the largest cluster from a simulation of IFC percolation at its own percolation threshold. Notice the difference in the number of FC sites in these figures. Partially coordinated sites dominate in ordinary percolation while FC sites dominate in IFC percolation.

The average number of FC and PC sites in ordinary, FC, and IFC percolation are plotted in Fig. 5. Linear least squares fits indicate that in ordinary percolation the PC sites are more than $`6`$ times as many as the FC sites, while the ratio in FC percolation is about $`1.5`$, and in IFC percolation it is only about $`0.86`$. A trend that can easily be seen here is that the number of FC sites increases as we progress from ordinary to FC and then to IFC percolation. Although only a local full-coordination correlation is imposed, it clearly results in a change on the global scale. This may be reflected in some qualitative differences in a physical property if it is sensitive to the size of connected clusters of FC sites. For example, if a percolation model is used to describe a disordered substrate on which diffusive physics takes place (such as in catalysis), what is expected within a given time scale (or frequency range) would be directly affected by the size of the high-connectivity or FC regions of the substrate where fast diffusion occurs compared to less well connected areas. In such cases, it is expected that there would be large differences depending on whether the appropriate percolation model to use is ordinary, FC, or IFC percolation. The latter cases could be a better representation if the culling of partially coordinated sites naturally occurs in the particular problem, say, due to desorption.

The distribution of cluster sizes (of all sites, i.e., both FC and PC sites) is shown in Fig. 6. Because of the finite size effects, only the intermediate size ranges (about $`10^2`$ to $`10^4`$ sites) are consistent with the cluster size scaling $`N_ss^\tau `$ with the ordinary percolation value of $`\tau =187/912.055`$, where $`s`$ is the cluster size and $`N_s`$ is the number of clusters having size $`s`$. Very small $`s`$ values are clearly out of the asymptotic scaling regime while the region of $`s𝒪(L^2)`$ is also limited by the finite size of the grid. Nonetheless, it is important to note that there is no significant difference between ordinary, FC, and IFC percolation for all $`s`$ above about $`10^2`$. This suggests strongly that the asymptotic cluster size scaling (which is known to occur for ordinary percolation) also occurs for both FC and IFC percolation.

In making the plot for Fig. 6, we didn’t distinguish between PC and FC sites in the clusters. If we instead just consider FC sites and make a semi-log plot of the number of clusters formed by these sites against the corresponding cluster sizes, we should be able to roughly estimate the average cluster sizes formed by FC sites alone. Shown in Fig. 7(a) is $`N_s`$ against $`s`$ of clusters formed by FC sites in ordinary percolation. The line shown is not a best fitting line for all of the available points but instead one that best fits only the four middle points and is shown merely to obtain a crude idea of the typical size scale of FC clusters. For ordinary percolation, we found the average FC cluster size to be around $`s_02`$ FC sites. Shown in Figs. 7(b) and 7(c) are data for FC and IFC percolation, respectively. We found that around $`s_118`$ FC sites comprise the average cluster in FC percolation while around $`s_272`$ FC sites make up the average FC cluster in IFC percolation. This means that in order to map the FC and IFC percolation problems to the same length scale as ordinary percolation, a length rescaling by $`\sqrt{s_1/s_0}3`$ and $`\sqrt{s_2/s_0}6`$ must be performed, respectively. Thus, although both static and dynamic critical exponents seem to remain the same as the full-coordination rules are iteratively applied, the local connectivity greatly changes or, equivalently, the length scale which would preserve the local connectivity changes.

## IV COMPARISON WITH PREVIOUS WORK ON FC PERCOLATION

In this work, we follow the prescription described in the Introduction closely in order to generate the ensemble of FC and then IFC clusters. The only limitation in the Monte Carlo portion of our work are that the grid we use to generate the clusters must be finite and that a particular boundary condition must be selected. Our largest grid is $`2048`$x$`2048`$ and clearly, this is not very large by today’s standards. Since we use finite size scaling in the determination of both the static and dynamic critical exponents, the finiteness of the grid (and cluster) is in fact necessary; yet, there is always a possibility that an asymptotic regime of large sizes may not have been reached. Also, we use the free boundary condition with no direct connectivity across the opposite edges of the grid, which may accentuate the finite size limitations. This type of ensemble is the most straightforward extension of the standard Monte Carlo simulation of ordinary percolation to the FC and IFC percolation models, and may be described as static since there is no dependence on the order in which sites are generated or the connectivities checked. Indeed, the static exponents $`\gamma `$, $`\gamma /\nu `$, and the like, agree closely with the corresponding ordinary static percolation values, and, if the length scale is identified with the grid edge length, then the cluster’s fractal dimensions should also agree closely with that of ordinary percolation.

In our previous work , we studied the static and dynamic critical behavior of FC percolation using the ensemble of clusters of specified number of sites, say, $`S^{}`$. Those clusters were grown from a seed site by applying a breadth-first algorithm twice; i.e., first generating an ordinary percolation cluster in nearest neighbor shells to a sufficiently large size (say, $`5S^{}`$), and then identifying and joining only the fully coordinated sites in a second searching process from the same seed and then stopping the process when the size of the resulting connected cluster reaches the predetermined value $`S^{}`$. Now, if we make the intermediate cluster of the first search process very large and if we do not terminate the second search process arbitrarily when $`S^{}`$ is reached, but rather continue the search until the natural boundary of the FC cluster is found, then the resulting cluster ensemble should be exactly the same as that which would result in the static ensemble of FC percolation clusters generated on an infinite grid. However, since the first stage is stopped at at most $`5S^{}`$ and since the second stage is forcibly stopped exactly at $`S^{}`$, which is a predefined value, we introduced a kinetic nature to the cluster ensemble. This may be somewhat similar to the differences in the degree of intrinsic anisotropy in static and kinetic ensemble of percolation clusters studied previously  (called equlibrium versus growing clusters there).

An illustration of the cluster geometry differences arising between a static cluster and a growing one is shown in Fig. 8. In the static ensemble, every site of the spanning, connected component of the ordinary, FC, or IFC percolation models’ realizations within a given length scale $`L`$ (defined to be the grid edge length) is a part of the cluster. In the growing cluster, on the other hand, only those sites which are within a certain path length along the nearest neighbor connections from the seed site are counted, independent of how close they are by straight line distance to the seed. Also, since cluster size $`S^{}`$ is the control parameter for this ensemble, the length scale $`L`$ is not precisely fixed but rather obtained as an average for a given value of $`S^{}`$. Exactly which sites are counted in a cluster of length scale $`L`$ and which ones excluded do clearly depend on the particular search order that is employed. Thus, typically, for the same $`L`$, the number of sites of the percolating cluster is smaller in the growing cluster ensemble than in the static one. This tends to reduce both the fractal dimension $`d_f`$ and the walk dimension $`d_w`$. If the effect is more pronounced for $`d_w`$, then it may explain why an apparent measurement of the spectral dimension $`d_s=2d_f/d_w`$ was larger for FC in the growing cluster ensemble. Since the search algorithm is repeated to obtain an IFC cluster, the effect will be even larger in IFC than in FC percolation, as well as it is larger in FC than in ordinary percolation.

## V SUMMARY

In summary, we have introduced and studied an extension to ordinary site percolation on a square lattice which we call iterated fully coordinated percolation. The static and dynamic critical behaviors of this model, as well as that of the ordinary and fully coordinated percolation models, were determined using methods such as Monte Carlo simulations, finite-size scaling, and normal mode analyses. All clusters used in our simulations were generated statically in grids of predetermined sizes on a square lattice, as opposed to growing them dynamically from a seed according to a dynamic rule. We found that the sizes of subregions of full connectivities in an incipient infinite cluster greatly increases as we simulate from ordinary to fully coordinated and then to iterated fully coordinated percolation. The length scales of such subregions increase by a factor of two to three with each iteration. Although this implies a large increase in the local connectivities, our results also indicate that the ordinary, fully coordinated, and iterated fully coordinated percolation models have identical static and dynamic critical exponents within our statistical errors and thus we believe that they belong to the same universality class, both statically and dynamically.