text:
# Electronic structure of NiS1-xSex

## I Introduction

Electronic structure of 3$`d`$-transition metal compounds has been an active area of research for many decades, owing to their rich and diverse physical properties, often exhibiting exotic transport and magnetic properties. In particular, the electronic structure of the hexagonal nickel sulphide (NiS) has attracted a lot of interest  due to the fact that it has a first-order phase transition as a function of temperature at about 260 K. The high-temperature phase is a Pauli paramagnetic metal, transforming isostructurally to an antiferromagnetic state with nearly temperature independent resistivity . The high temperature phase is metallic with a resistivity in the order of 10-5 Ohm-cm  which is comparable with other highly metallic inorganic compounds of 3$`d`$ transition elements . The metallic state can be stabilized by the application of pressure or substitution with Se for S. A pressure of more than 20 kbar  or a Se percentage of more than 13%  suppresses the electronic phase transition, stabilizing the metallic phase down to the lowest temperature. At room temperature, the hexagonal NiS crystallizes in the NiAs structure having space group $`P6_3/mmc`$ ($`D_{6h}^4`$) and lattice parameters, $`a`$ = 3.440 Å and $`c`$ = 5.351 Å . The Ni2+ is in an octahedral environment of sulphur atoms in NiS and these NiS6 octahedra are edge-shared within the $`ab`$-plane and face-shared along the $`c`$-axis, forming a three-dimensional network of Ni-S connectivity in this system.

In spite of its interesting properties, there have been few studies to understand the electronic structure of these compounds. Band structure calculations  indeed yield a metallic ground state for the high-temperature phase. However, the details of the electronic structure, as shown in the present work, cannot be entirely described within the band theory due to the presence of substantial correlation effects even in the metallic phase. Different configuration interaction models have been employed in order to explain the valence band  and the transition metal core level spectra ; however, there has not been any attempt to describe both using a single model and with a single set of parameters. In the present paper, we investigate the electronic structure of metallic NiS1-xSex system using $`x`$-ray photoemission (XP), ultra-violet photoemission (UP), and bremsstrahlung isochromat (BI) spectroscopic measurements in conjunction with ab-initio band structure calculations as well as parameterized many-body calculations, in order to provide a consistent and quantitative description of the electronic structure and the electronic parameter strengths controlling the physical properties of this system.

## II EXPERIMENTAL

Polycrystalline samples of NiS1-xSex with $`x`$ = 0.0 and 0.15 used for the present study were prepared following the standard solid state reaction techniques reported in the literature ; $`x`$-ray diffraction patterns as well as resistivities of the samples were found to be in agreement with the reported data in the literature . Spectroscopic measurements were carried out in a combined VSW spectrometer with a base pressure of 2$`\times 10^{10}`$ mbar equipped with a monochromatized Al K$`\alpha `$ $`x`$-ray source, He discharge lamp and an electron gun. XP, BI, and UP spectroscopic measurements were performed on the samples with an overall instrumental resolution of better than 0.8 eV, 0.8 eV, and 90 meV, respectively. BI spectra have been recorded by monochromatizing the emitted photons at 1486.6 eV. The sample surface was cleaned in situ periodically during the experiments by scraping with an alumina file and the surface cleanliness was monitored by recording the carbon 1$`s`$ and oxygen 1$`s`$ XP signals. The reproducibility of the spectral features was confirmed in each case. The binding energy was calibrated to the instrumental Fermi-level which was determined by recording the Fermi-edge region from a clean silver sample. Experimental spectra of NiS reported here were collected for the highly conducting state at room temperature ($`295`$ K), while those of NiS0.85Se0.15, which does not exhibit any transition, were collected at $`120`$ K. Using high resolution photoemission spectroscopy, it has been established  that the electronic structure close to $`E_F`$ changes across the electronic transition at $`260`$ K. However, our low-resolution spectroscopic results show hardly any change in the spectra across the transition either in the valence band or in the core levels, in agreement with previous studies  establishing similar gross electronic structures with comparable electronic parameter strengths for both phases across the electronic transition.

Scalar relativistic linearized muffin-tin orbital (LMTO) band structure calculations have been performed within the atomic sphere approximation (ASA) for the paramagnetic phase of NiS with the real crystal structure. The unit cell has two formula units and the sphere radii for Ni and S were 2.4 and 3.12 au, respectively. In this calculation, convergence was obtained self-consistently using $`s`$, $`p`$ and $`d`$ basis at each atomic sphere with 186 $`k`$ points in the irreducible part of the Brillouin Zone. In order to understand and compare the experimental electronic structure with that obtained from the band structure calculations, it is necessary to include the effects of photoemission cross-sections of various levels. Thus, the theoretical spectra based on the band structure calculations were constructed by multiplying the various partial densities of states (DOS) with the corresponding cross-sections computed within the formalism of ref. . The cross-section weighted DOS was then convoluted with a Lorentzian with energy dependent full-width at half maximum, to account for the lifetime broadening and a Gaussian to describe the instrumental resolution broadening.

Core level and valence band (VB) spectra were calculated for a NiS6 cluster with the geometry in accordance with the bulk crystal structure of this compound, within a parameterized many-body multi-band model including orbital dependent electron-electron (multiplet) interactions . The multiplet interactions within the $`3d`$ manifold were expressed in terms of Slater integrals $`F_{dd}^0`$, $`F_{dd}^2`$ and $`F_{dd}^4`$ and those within the $`2p3d`$ manifold were expressed in terms of $`F_{pd}^0`$, $`F_{pd}^2`$, $`G_{pd}^1`$ and $`G_{pd}^3`$. The values $`F_{dd}^2`$ (=9.79 eV), $`F_{dd}^4`$ (=6.08 eV), $`F_{pd}^2`$ (=6.68 eV), $`G_{pd}^1`$ (=5.07 eV) and $`G_{pd}^3`$ (=2.88 eV) were taken to be 80% of the Ni atomic Hartree-Fock values. The monopole terms, $`F_{dd}^0`$ and $`F_{pd}^0`$ were treated as parameters and therefore, varied to fit the experimental spectrum. It is often more convenient to discuss the electron-electron interaction strengths in terms of the multiplet averaged quantities, $`U_{dd}`$ and $`U_{pd}`$, where

$$U_{dd}=F_{dd}^0\frac{2}{63}(F_{dd}^2+F_{dd}^4),$$
(1)
$$U_{pd}=F_{pd}^0\frac{1}{15}G_{pd}^1\frac{3}{70}G_{pd}^3.$$
(2)
The multiplet averaged $`U_{pd}`$ term was set to have a value 1.2 times that of $`U_{dd}`$ , thus reducing the number of freely adjustable parameters in the model. The hopping interaction strengths between Ni 3$`d`$ and S 3$`p`$ were expressed in terms of the Slater-Koster parameters $`(pd\sigma )`$ and $`(pd\pi )`$ and that between the different bonding S 3$`p`$ orbitals were expressed in terms of $`(pp\sigma )`$ and $`(pp\pi )`$ . We fixed the ratio between $`(pd\pi )`$ and $`(pd\sigma )`$ to be -0.5, while that between $`(pp\pi )`$ and $`(pp\sigma )`$ to be -0.2, in these calculations.

The calculations were performed in the symmetry adapted ($`t_{2g}`$ and $`e_g`$) basis including the transition metal 3$`d`$ and the bonding sulphur 3$`p`$ orbitals. In the calculation for the valence band spectrum, the S 3$`p`$ spectral contribution was also calculated within the same parameterization. As the size of the Hamiltonian is large, the Lanczos method was used to evaluate the spectral function and the calculated one-electron removal spectra were appropriately broadened to simulate the experimental spectra. In the Ni 2$`p`$ core level calculation, Doniach-Sunjic line shape function was used for broadening the discrete energy spectrum of the cluster model, in order to represent the asymmetric line shape of core levels from these highly metallic compounds, consistent with other core levels in the system. In the case of VB calculations, energy dependent Lorentzian function was used for the lifetime broadening. Other broadening effects such as the resolution broadening and solid state effects were taken into account by convoluting the spectra with a Gaussian function. The broadening parameters were found to be consistent with values used for similar systems . Taking into account the different atomic cross-sections for the Ni 3$`d`$ and S 3$`p`$ states, it is necessary to calculate a weighted average of these two contributions to the valence band. The atomic cross-section ratio  between S $`3p`$ and Ni $`3d`$ states ($``$ 0.17) is not appropriate in this context, since solid-state effects alter this ratio significantly. It was found that S 3$`p`$/Ni 3$`d`$ cross-section ratio of approximately 5.5 times that obtained from the atomic calculations gives the best result for the cluster calculations.

## III RESULTS & DISCUSSIONS

The total density of states (DOS) for NiS along with the Ni $`d`$ and S $`p`$ partial DOS obtained from the band structure calculations for the paramagnetic phase are shown in Fig. 1. The calculated DOS is in good agreement with previous reports . The Fermi level lies close to a minimum in the DOS, and there is a finite DOS at the $`E_F`$, consistent with the metallic nature of the system. The total DOS between -3 and 1 eV is dominated by Ni $`d`$ contributions with the energy region around the most intense peak in the DOS at about -1.65 eV dominated by the Ni $`t_{2g}`$ states and the features at higher energies ($`1`$ eV) arising primarily from the $`e_g`$ states. As a consequence of the high covalency of the system, the Ni 3$`d`$ ($`e_g`$) states have substantial contribution from S 3$`p`$ states. These states constitute the antibonding states due to Ni $`d`$-S $`p`$ interactions. Their bonding counterparts are located between -5 and -3 eV and have dominant S 3$`p`$ contributions. The non-bonding states of the system, centered around -7 eV is weaker in intensity and is stabilized at a higher binding energy compared to Ni $`d`$-S $`p`$ bonding states, which is a generic feature as seen in the case of the transition metal monosulphides .

We compare the experimental XP valence band and BI spectra corresponding to the occupied and unoccupied parts of the electronic structures for NiS and NiS0.85Se0.15 in Fig. 2(a) and (b), respectively. The calculated spectra for NiS, obtained using various partial DOS as explained in the previous section, are also shown for comparison. The XP and BI spectra of NiS and NiS0.85Se0.15 show finite spectral intensity at $`E_F`$, indicating their highly metallic conductivity. It is clear from the experimental spectra in Figs. 2(a) and (b) that the substitution of Se in place of S does not cause any perceptible change in the XP valence band spectra as well as in the BI spectra. This indicates that the occupied as well as the unoccupied parts of the electronic structure in these two compounds are similar, in agreement with the similar magnetic and transport properties of the two compounds.

Fig. 2(a) shows that the peak in the experimental spectra around 1.3 eV binding energy (marked A) essentially consists of Ni $`3d`$ states, while that at $`4`$ eV binding energy (marked B) is due to the S 3$`p`$ contribution. However, the calculated spectrum does not describe the experimental spectra accurately, as the energy positions of the different features, A and B are not correctly predicted by the calculation, but are shifted to higher binding energies compared to the experiment. Additionally, the experimental feature A is narrower compared to the calculated one. Moreover, the structure around 7 eV binding energy (marked C) in the experimental spectra is absent in the theoretical spectrum. In the unoccupied part of the electronic structure, the peak at around -0.7 eV in the experimental BI spectra is contributed mainly by the unoccupied eg states of the Ni $`3d`$ manifold, strongly hybridized with the S $`3p`$ bands, as seen from Fig. 1. The strong suppression in the S 3$`p`$ contribution in the calculated spectrum (Fig. 2(b)) compared to the partial DOS (Fig. 1) in the vicinity of $`E_F`$ arises from a substantially lower photoemission cross-section for the S 3$`p`$ states compared to Ni 3$`d`$ at Al K$`\alpha `$ $`x`$-ray energies. The broad spectral features above 4 eV are due to the various high energy unoccupied states in the system, such as Ni $`4s,4p`$ and S $`3d`$. Similar to the case in XP spectra, the band structure calculation, while providing a reasonable description of the overall spectral features observed in the BI spectra, fail to reproduce the energy position of the feature A in Fig. 2(b) by approximately 0.3 eV. These failures of the band theory to explain the electronic structure of this system indicates that electron correlation effects are important even in the paramagnetic metallic phase.

The valence band spectra of NiS and NiS0.85Se0.15 recorded using different photon energies, namely 1486.6 eV (XPS), 40.8 eV (He II) and 21.2 eV (He I) are compared in Fig. 3. The various features are marked as A, B and C and represent the same features as in Fig. 2(a). At all the photon energies, there is a finite spectral intensity at $`E_F`$. Besides providing a much higher resolution compared to the XPS, UP spectra with He I and II radiations allow for a wider variation in the relative photoemission cross-sections of Ni 3$`d`$ and S 3$`p`$, leading to an easy identification of their relative contributions to the valence band spectra. The marked loss in the intensity of the feature B for a moderate increase of the photon energy between He I and He II radiations arises from the Cooper minimum in the S $`3p`$ photoionization cross-section at around 50 eV, clearly establishing the dominance of S 3$`p`$ contributions in the energy range of the feature B, with feature A having predominantly Ni 3$`d`$ contributions, in agreement with previous works . The valence band spectra for NiS and NiS0.85Se0.15 taken at various photon energies are very similar, showing the similarity in their electronic structure within the solid solution, NiS1-xSex for $`x0.15`$.

The transition metal $`2p`$ core level spectra in transition metal compounds often reveal important information on their electronic structures. The Ni $`2p`$ core level spectra for NiS and NiS0.85Se0.15 are shown in the inset of Fig. 4; the two spectra are very similar to each other. Each spectrum consists of spin-orbit split, $`2p_{3/2}`$ and $`2p_{1/2}`$ peaks at 853.4 eV and 870.7 eV binding energies, respectively with pronounced satellite features at 860 eV and 876 eV, suggesting that electron correlations are important in the system. The satellite intensity relative to the main peak is considerably more in the $`2p_{1/2}`$ region compared to that in the $`2p_{3/2}`$ region. While it is indeed possible to have different shapes and intensity for the two satellite features accompanying $`2p_{3/2}`$ and $`2p_{1/2}`$ features due to details of the multiplet interactions, such a strong variation between the two satellites as observed here is unexpected. In order to explore the possibility of strong variations in the inelastic background function, normally assumed to vary smoothly, we have performed electron energy loss spectroscopy (EELS) on these samples, with the same primary energy as that of the Ni $`2p`$ core level peak. The inelastic loss background appropriate for the Ni 2$`p`$ photoemission spectrum was constructed from the experimental EELS spectrum by convoluting it with a Gaussian to obtain similar width for the elastic peak as that in the core level main peak. The region near the elastic peak was replaced by an integral background function. The total inelastic background was then generated by adding the loss functions corresponding to the $`2p_{3/2}`$ and $`2p_{1/2}`$ features with the intensities given by the degeneracy ratio (2:1) of the core levels. The resulting background function is shown in the inset of Fig. 4 as a solid line. We find that there is an intense and structured contribution to the background function arising from a plasmon band overlapping the $`2p_{1/2}`$ satellite region and resulting in the anomalously large satellite intensity in the $`2p_{1/2}`$ region compared to that in the $`2p_{3/2}`$ region. It is also seen that at around 856 eV, there is a peak-like structure in the inelastic background; this appears at the same position as that of a weak intensity feature at about the same energy in the core level experimental spectrum as shown by the vertical dashed lines in the inset of Fig. 4. This structure in the inelastic scattering background could have its origin from the inter-band $`pd`$ transitions. Thus, it appears that this weak intensity shoulder in the experimental Ni 2$`p`$ spectrum is largely due to this extrinsic loss feature in the background function. Our explanation is contrary to the earlier interpretation  of this feature as arising from intrinsic processes.

In order to obtain a quantitative description of the electronic structure of this system, we have calculated the Ni $`2p`$ core level and valence band (VB) spectra within a single model involving a NiS6 cluster, retaining the local geometry of the real solid. The calculated Ni $`2p`$ spectrum (solid line) for the parameter set S-Ni $`(pd\sigma )`$ = -1.4 eV, S-S $`(pp\sigma )`$ = 0.7 eV, $`\mathrm{\Delta }`$ = 2.5 eV and $`U_{dd}`$ = 4.0 eV, including the experimentally determined inelastic background is superposed on the experimental spectrum (open circles) in the main panel of Fig. 4, exhibiting a good agreement; the calculated spectrum without any broadening is presented as bar diagrams. There is an underestimation of the intensity around the shoulder to the main peak near 857 eV and around 883 eV binding energies. These small discrepancies are possibly due to the differences between the experimental EELS spectrum used for generating the background function and the actual inelastic loss background for the photoemission spectrum. The previous estimates  of the various parameter strengths in NiS obtained from a model that included a fictitious conduction band in addition to Ni 3$`d`$-S 3$`p`$ basis within the cluster model used for the core level calculation, are $`(pd\sigma )`$ = -1.47 eV, $`\mathrm{\Delta }`$ = 2.5 eV, and $`U_{dd}`$ = 5.5 eV. Thus, the present estimates differ significantly only for $`U_{dd}`$, for which we have obtained a smaller value.

The $`(pd\sigma )`$ values estimated for NiS is similar to that estimated for other divalent nickel sulphides, for example, $`(pd\sigma )`$ = -1.5 eV for NiS2 and BaNiS2 , and for trivalent Ni oxide systems, such as LaNiO3 ($`(pd\sigma )`$ = -1.57 eV ). However, divalent Ni oxide systems, such as NiO, typically exhibit a Ni 3$`d`$-O 2$`p`$ $`(pd\sigma )`$ = -1.25 eV . This decrease in the strength of $`(pd\sigma )`$ for the divalent oxides compared to trivalent ones is due to an increase in the Ni-O distances in the divalent oxides. The similarity in the $`(pd\sigma )`$ of nickel sulphides and the trivalent Ni oxides, in spite of the differences in their Ni-ligand bond length ($`d_{NiO}`$ 1.93 Å for LaNiO3 compared to $`d_{NiS}`$ 2.39 Å in NiS) arises from the fact that S 3$`p`$ orbitals are much more spatially extended in nature than O 2$`p`$ and hence have a larger overlap integral with Ni 3$`d`$ orbitals. In general, $`\mathrm{\Delta }`$ is expected to be smaller for sulphides compared to oxides, since the O $`2p`$ levels are more stable than the S $`3p`$ levels; for example, the estimated $`\mathrm{\Delta }`$ for NiO is 5.5 eV  compared to 2.5 eV in NiS. The value of $`U_{dd}`$ for NiS (4.0 eV) is much smaller than that in NiO (6.5 eV) possibly due to a more efficient screening in NiS, arising from the smaller charge transfer energy, $`\mathrm{\Delta }`$ and larger hopping strength, $`(pd\sigma )`$, leading to stronger covalency effects.

We have analyzed the ground state wavefunction of NiS corresponding to the estimated parameter strengths in terms of the various electron configurations. The ground state of the system was found to consist of 60.7%, 35.8% and 3.4% of $`d^8`$, $`d^9\underset{¯}{L}^1`$, and $`d^{10}\underset{¯}{L}^2`$ configurations with a high-spin configuration ($`S=1`$). The average value of the $`d`$occupancy $`(n_d)`$ is found to be 8.43, showing a highly covalent ground state of the system. In the case of NiO, a typical divalent oxide of Ni, $`n_d`$ is found to be $`8.16`$ . We have analyzed the characters of the final states of the system responsible for different features in the experimental spectrum, in order to understand their origins. The analysis was carried out for the final state energies marked 1-10 in Fig. 4. The different contributions to the final states from various electron configurations ($`d^8`$, $`d^9\underset{¯}{L}^1`$, $`d^{10}\underset{¯}{L}^2`$) are listed in Table. I. These features can be grouped into three different regions; the main peak region, 852-856 eV (labeled 1-3 in Fig. 4); intense satellite region, 858-861 eV (labeled 4-6), and weak satellites in the region 863-866 eV (labeled 7-10). The first group of features in the main peak region has a dominant $`d^9\underset{¯}{L}^1`$ character as seen from the table, which are the “well-screened” states of the system and correspond to one ligand (sulphur) electron being transferred to the Ni site to screen the Ni $`2p`$ core-hole potential created by the photoemission process. The second group of features, has a mixed character with significant contributions from all the configurations. This is in contrast to the case of NiO where the intense satellite structure results primarily from the “poorly-screened” $`d^8`$ configuration. The third group of features has a dominant $`d^{10}\underset{¯}{L}^2`$ contribution for the lower energy region (labeled 7 and 8) and a dominant $`d^8`$ character for the higher energy region (labeled 9 and 10 in Fig. 4), with a considerable amount of $`d^9\underset{¯}{L}^1`$ character for the entire region.

The XP valence band spectrum of NiS has also been calculated within the same model. An important feature of the present calculation is that the S $`3p`$ contribution to the VB spectrum has been calculated on the same footing as the contribution from the Ni 3$`d`$ within the same model without the need to adjust either the shape or the energy position of the ligand contribution to the valence band region, thereby enhancing the reliability of the estimated parameter strengths. The calculated spectrum (solid line) along with the Ni $`3d`$ (dashed line) and S $`3p`$ (dot-dash line) contributions for NiS are shown superimposed on the experimental data (open circles) in Fig. 5. An inelastic scattering background function (dotted line) is also included in the total calculated spectrum. The calculated contributions from the Ni $`3d`$ to the total spectrum for NiS is also shown without any broadening effect as an energy stick diagram in Fig. 5. The parameter set used for the VB calculation is similar to that used for the core level calculation, but with a Ni 3$`d`$-S 3$`p`$ $`(pd\sigma )`$ of -1.2 eV instead of -1.4 eV. Such minor adjustments of parameter strengths between the simulations for the core and valence level photoemission spectra, are of common occurrence. The present estimates of the electronic structure parameters are similar to that obtained from the previous valence band calculations  within the cluster-model. The agreement between the experimental spectrum and the calculated spectrum in the present study (see Fig. 5) is remarkable over the entire energy range. The main peak in the valence band spectrum in Fig. 5 at about 1.3 eV arises essentially from Ni 3$`d`$ photoemission contribution, though there is a small contribution arising from S 3$`p`$ states also due to the covalent mixing of states. It is possible to delineate such contributions only in the present approach where S 3$`p`$ contribution is also calculated within the same model. Within the traditional approach, S 3$`p`$ contribution would be approximated by a single Gaussian with its peak position, width and intensity being freely adjusted to obtain the best possible description. Clearly the present approach is more reliable and satisfactory. A second, and more intense, S 3$`p`$ contribution is responsible for the experimental spectral feature near 4 eV binding energy. The shoulder appearing between about 6 and 9 eV is once again dominated by the Ni 3$`d`$ contributions. This observation along with the fact that this feature is not in agreement with the ab-initio band structure results (see Fig. 2), clearly suggests that this feature is driven by correlation effects. The analysis of the ground state wavefunction for the parameter set used for the VB spectrum shows that the ground state has 64.9%, 32.6% and 2.5% of $`d^8`$, $`d^9\underset{¯}{L}^1`$ and $`d^{10}\underset{¯}{L}^2`$ characters giving rise to an average $`d`$occupancy ($`n_d`$) of 8.38 electrons with a high-spin configuration for the Ni2+, similar to that obtained from the core level calculation. The parameter strengths estimated here place NiS1-xSex samples in $`pd`$-metal region of the ZSA phase diagram .

The results of the character analysis of the final states labeled 1-11 in the figure corresponding to the Ni $`3d`$ contribution in the VB spectrum of NiS are shown in Table. II. On the basis of this analysis, the spectral features can be grouped into three regions, the main peak region (0-6 eV, labeled 1-6), the strong satellites in the 6-8.5 eV range (labeled 7-9) and weaker satellites beyond 8.5 eV (marked 10 and 11). The final states in the main peak region predominantly consist of $`d^8\underset{¯}{L}^1`$ states with a non-negligible contributions from $`d^7`$ and $`d^9\underset{¯}{L}^2`$ configurations. This is similar to the scenario as observered in the case of other charge-transfer systems, like NiO. In the satellite region (6-8.5 eV), the final states are dominated by $`d^7`$ and $`d^9\underset{¯}{L}^2`$ configurations. The weak intensity features at higher energies are dominated by $`d^9\underset{¯}{L}^2`$ configuration with a substantial contribution also from $`d^8\underset{¯}{L}^1`$ states. It is to be noted that the first ionization state, marked 1 in the figure, is not at the Fermi energy. This is the well known artifact of finite cluster-model calculations which invariably exhibit a gap in the charge excitation spectrum due to the discrete eigen-spectrum of a finite cluster. It is believed that this ionization state acquires a band width in the thermodynamic limit closing the discrete gap and overlaping the $`E_F`$. A dispersion width of less than 1 eV is seen to be sufficient to bring this about; such a dispersion is consistent with our results obtained from LMTO band structure calculations for the infinite solid. While we have not performed an analysis of the orbital symmetry, our analysis for the first ionization state shows it to be a spin-doublet state.

It is evident that the shoulder-like structure at around 7 eV in the experimental VB spectra is primarily due to the satellite structure of Ni $`3d`$ states and is expected to have some non-bonding S $`3p`$ contribution as seen from its photon energy dependence (see Fig. 3). This explains the inability of the band structure calculations to explain this feature, as the origin of this feature is primarily due to the electron correlations present in the system, giving rise to strong satellite structures. The same correlation effects are also responsible for shifting the main peak to a lower energy and making it narrower compared to the band structure results (see Fig. 3). This satellite structure, with a dominant $`d^7`$ character, is the spectral signature of the lower Hubbard band. This observation along with the highly conducting Pauli paramagnetic property of these compounds suggests that the Hubbard side-bands exist deep inside the metallic regime. This is in accordance with the predictions of the dynamical mean field theory (DMFT) calculations based on the Hubbard model , that well defined Hubbard bands continue to exist in the spectral functions even far away from the Mott transition in the metallic region, though the Mott-Hubbard gap collapses by the transfer of spectral weight from the Hubbard side-bands to the coherent spectral features near $`E_F`$ region.

In conclusion, we have investigated the electronic structure of the system, NiS1-xSex, using various electron spectroscopic techniques to elucidate the underlying electronic structure. The substitution of Se in place of S does not show appreciable changes in the electronic structure of the system, indicating similar electronic parameters for this solid solution. The band structure results do not describe well the details of the spectral features in the paramagnetic phase of the system, while a parameterized many-body multi-band model is found to be successful in describing the core level and valence band spectrum of the NiS. The core and valence band spectral calculations of paramagnetic metallic phase of NiS show that it is a correlated metal with a highly covalent character. The asymmetric line shape as well as the weak feature in the core level spectra have been ascribed to extrinsic loss processes in the system. The valence band calculations show that lower Hubbard band exists well inside the $`pd`$ metallic regime, as predicted by the DMFT calculations. The electronic parameter strengths obtained for the compounds of this solid solution show that they belong to the $`pd`$ metallic regime of the ZSA phase diagram.

## IV Acknowledgments

The authors thank Professor C. N. R. Rao for continued support and the Department of Science and Technology, and the Board of Research in Nuclear Sciences, Government of India, for financial support. SRK thanks the Council of Scientific and Industrial Research, Government of India, for financial assistance. DDS thanks Dr. M. Methfessel, Dr. A. T. Paxton, and Dr. M. van Schiljgaarde for making the LMTO-ASA band structure program available. The authors also thank Professor S. Ramasesha and the Supercomputer Education and Research Center, Indian Institute of Science, for providing the computational facility.

## V figure captions

Fig. 1. Density of states (DOS) of NiS from the band structure calculation for the paramagnetic state. The thick solid line represents the total DOS. The thin solid and dashed lines are the Ni 3$`d`$ and S 3$`p`$ partial DOS, respectively.

Fig. 2. (a) Experimental valence band spectra using Al K$`\alpha `$ for NiS and NiS0.85Se0.15 along with the calculated spectrum obtained from band structure calculations. The thick solid line represents the total spectrum, while the thin solid and dashed lines correspond to the contributions to the calculated spectrum from Ni 3$`d`$ and S 3$`p`$ partial DOS, respectively. (b) Experimental BI spectra for NiS and NiS0.85Se0.15 along with calculated spectrum obtained from the band structure calculations. Various contributions to the total spectrum are represented by different line types.

Fig. 3. Experimental valence band spectra for NiS and NiS0.85Se0.15 using 1486.6 eV (XPS), 40.81 eV (He II) and 21.2 eV (He I) photon energies. Various features in the spectra are shown by the vertical lines and are labeled as A, B, and C (see text).

Fig. 4. Experimental Ni $`2p`$ spectrum (open circles) along with the calculated spectrum (solid line) for NiS obtained from the cluster calculation. Various final states of the cluster calculation and the corresponding intensity contributions without any broadening are shown as the bar diagram. Inset shows Ni $`2p`$ core level spectra for NiS and NiS0.85Se0.15 along with the inelastic scattering background function obtained from EELS for NiS. The dashed vertical line shows the structure in the background function responsible for the feature in the core level spectra around 856 eV.

Fig. 5. The experimental VB spectrum (open circles) along with the calculated spectrum (solid line), Ni $`3d`$ component (dashed line), S $`3p`$ component (dot-dashed line) and the integral background (doted line) are shown for NiS. The final states of the calculation and the corresponding intensities without any broadening are shown as the energy stick diagram.


format: text

text:
# A Galactic Bar to Beyond the Solar Circle and its Relevance for Microlensing

## 1 Introduction

A somewhat unexpected discovery of microlensing experiments towards the Galactic Bulge was the need to take into account a bar-like structure there (see, e.g. Paczyński 1996). The precise nature of this bar - its composition, extent and evolution - is only sketchily known. This is an area where progress is likely to come from the combination of microlensing results with other types of observations. In the present paper we summarize the results of an analysis of the Galactic kinematics of Mira variables based on infrared photometry, radial velocities, and Hipparcos parallaxes and proper motions (Whitelock, Feast & Marang 2000, Whitelock & Feast 2000, Feast & Whitelock 2000b). This provides evidence for a rather extended bar-like structure in the Galaxy.

## 2 Background

Mira variables have a number of properties which makes them particularly interesting for Galactic structure studies 111Note that the present paper refers entirely to oxygen-rich Miras..

Bolometrically and in the near infrared, e.g. at $`K`$ (2.2$`\mu `$m), they are the brightest members of the old populations in which they are found (see for instance fig 1 of Feast & Whitelock 1987). These variables therefore define the tip of the AGB in these populations. Because of their brightness they can be studied with relative ease to large distances.

Observations in the LMC show that Miras follow a well defined, narrow period \- luminosity (PL) relation at $`K`$ or in $`M_{\mathrm{bol}}`$ (Feast et al. 1989). Thus Mira distances can be estimated from infrared photometry with good accuracy.

It has long been known that the kinematics of Galactic Miras are a function of period (see, e.g. fig 12 of Feast 1963). Their asymmetric drift gets numerically larger as one moves from the longer period Miras to the shorter period ones with a maximum of $`100\mathrm{km}\mathrm{s}^1`$ at a period of about 175 days. This is intermediate between the value usually quoted for the thick disc ($`3040\mathrm{km}\mathrm{s}^1`$, e.g. Freeman 1987) and that appropriate to halo objects ( $`230\mathrm{km}\mathrm{s}^1`$). This apparently intermediate population has generally been neglected in discussions of Galactic structure and kinematics. Judging from Miras in globular clusters this corresponds to a metallicity, \[Fe/H\], range from about –0.8 to about –1.3. As one moves to Miras of even shorter periods the numerical value of the asymmetric drift decreases, making the shortest period Miras similar in kinematics to Miras of much longer period. The relationship of this shortest period group to the other Miras was long a puzzle.

Another important property of the Miras is shown by those in globular clusters. Miras occur in relatively metal rich clusters. If there are more than one Mira in any cluster, their periods are all close together. Such clusters also contain semiregular (SR) variables. In any one cluster these have a range in periods, and in the period - luminosity plane they lie on an evolutionary track which terminates on the Mira PL relation (see Whitelock 1986, Feast 1989). The clusters containing Miras show that for these stars there is a good relation between period and metallicity (see fig 1 of Feast & Whitelock 2000a). Since the change of mass with metallicity is presumably small for these stars, the period - metallicity relation is primarily indicating a relation between metallicity and stellar radius at the tip of the AGB.

These various properties have been important in using Miras to study the Galactic Bulge itself. The work of Lloyd Evans (1976) and others showed that there was a wide distribution of periods amongst Miras in the Bulge. This, then, implies a wide range in metallicities. Indeed the metallicity distribution inferred from the Mira period distribution (Feast & Whitelock 2000a) is similar to that derived by Sadler et al. (1996) from Bulge K-type giants. The PL relation allows one to estimate the distances to individual stars in the Bulge (see, e.g. Glass et al. 1995). In this way Whitelock & Catchpole (1992) found that Bulge Miras at positive Galactic longitudes were nearer to us than those at negative longitudes. This was one of the first pieces of evidence for a barred structure of the Bulge.

## 3 Present Work

Our present study consists of extensive near infrared, $`JHKL`$, photometry of Galactic Mira-like variables222The sample contains a few stars classified as semiregular (SR) variables which have Mira-like properties (see Whitelock et al. 2000)., made at SAAO Sutherland, combined with Hipparcos astrometric and photometric data and with radial velocities from the literature.

Combining the infrared observations with the Hipparcos photometry allows us to see clearly the cause of the anomalous result for the asymmetric drift of the shortest period Miras discussed above. In the $`(HpK)_0`$ \- $`\mathrm{log}P`$ plane (where $`Hp`$ is the mean Hipparcos magnitude) there are two near parallel sequences (see Fig 1). A main (blue) sequence and a less populated, redder, sequence which is largely confined to the shorter periods and dominates at the very shortest periods. In the period range where the two sequences overlap their kinematics are quite different. The red sequence stars have kinematics which associates them with much longer period Miras. It is tempting therefore to identify some at least of the stars on this red sequence as analogous to the SRs in metal-rich globular clusters, which were discussed above, but at somewhat higher luminosities, and to suggest that they are evolving into longer period Miras. If this is the case we would expect the red sequence stars to be somewhat brighter than the blue sequence ones at a given period. This is because evolutionary tracks (at least in metal-rich globular clusters) lie above the Mira PL relation in the period - luminosity plane. Indeed the Hipparcos parallaxes do indicate that at a given period the red sequence stars are brighter than the blue sequence ones. It is also possible that some of the red sequence stars are pulsating in a higher mode than the blue sequence stars. The following discussion refers only to the main (blue) sequence stars. These stars may also be recognized from their infrared colours and these show that the Miras that have been found in globular clusters and the LMC belong to the blue sequence.

Hipparcos parallaxes of blue sequence stars were used to derive the zero-point of the PL relation. Using this calibration one can determine a distance modulus for the LMC. This is found to be $`18.64\pm 0.14`$ mag or slightly greater if a correction is necessary for a metallicity difference between the local and LMC Miras. This modulus agrees well with that determined from Cepheids and some other distance indicators (see the summary in Feast 1999).

We have studied the Galactic kinematics of the blue sequence Miras using Hipparcos proper motions and radial velocities from the literature. These radial velocities include radio observations of OH and other molecules. Distances were derived from the mean $`K`$ magnitudes and the Hipparcos calibrated Mira PL relation. However, the results are rather insensitive to the distance scale adopted.

The three components of the space motion of each star were derived in non-rotating cylindrical co-ordinates centred at the Galactic centre. These are $`V_R`$ radially outwards and parallel to the plane, $`V_\theta `$ at right angles to $`V_R`$ in the direction of Galactic rotation, and $`w`$ in the direction of the North Galactic Pole. The data were analysed by dividing the Miras into six (blue sequence) groups, according to period.

The results for $`V_\theta `$ confirm previous work. In the longest period group (mean period 453 days), $`V_\theta =223\pm 4\mathrm{km}\mathrm{s}^1`$. This is only slightly smaller than the circular velocity implied by the kinematics of Cepheids ($`231\mathrm{km}\mathrm{s}^1`$, Feast & Whitelock 1997). The value of $`V_\theta `$ drops as the period decreases reaching $`147\pm 14\mathrm{km}\mathrm{s}^1`$ in the group of mean period 175 days (17 stars). In this calculation one Mira in the period range 145 to 200 days (S Car) has been omitted, as it is moving on a highly eccentric retrograde orbit.

The main surprise of the analysis was in the values of $`V_R`$. This is $`+67\pm 17\mathrm{km}\mathrm{s}^1`$ in the group of mean period 175 days, whilst for the other groups there is a small net outward velocity.

There seem to be three possible explanations of this large mean outward radial motion in the 145 to 200 day group.

If the Galaxy is axially symmetric the result would imply that there is a general, axially symmetric, outward, radial motion of Miras in the period range 145 - 200 days. This seems unlikely.

Another alternative is that Miras in this period range belong to some galactic interloper. Whilst this cannot be ruled out we might ask why there are, in that case, no local stars of this type belonging to our Galaxy (such stars exist in the Galactic Bulge itself and in Galactic globular clusters). Also the mean motion perpendicular to the plane is small for this period group ($`w=12\pm 12\mathrm{km}\mathrm{s}^1`$) indicating that any interloper must be moving very closely parallel to the Galactic plane.

The most likely explanation would seem to be that there is a Galactic asymmetry leading to a deficit of incoming orbits in the solar neighbourhood for Miras in this period range. The values of $`V_\theta `$ and $`V_R`$ for this group then show that the major axes of their Galactic orbits are concentrated in the first quadrant of galactic longitude. All workers place the major axis of the Bulge-bar in this quadrant. A simple first-order calculation shows that the mean orbit has a major axis making an angle of $`17_4^{+11}`$ degrees with the Sun - Centre line. This agrees with the position angle of the Bulge-bar derived by Binney et al. (1991) from gas dynamics in the central region (16 degrees). It would be consistent with a number of other estimates which are in the 20 to 30 degree range, though values near 45 degrees derived by some workers agree less well. Rough estimates suggest that the orbits of some of these short period Miras are sufficiently eccentric to penetrate into the Bulge region itself at perigalacticum.

Figure 2 shows a plot of $`V_\theta `$ against $`V_R`$ for stars in the 145 to 200 day group. The curve is the relation that would be expected, on a simple model, if the major axes of the orbits of all the stars were aligned at an angle of 17 degrees to the Sun - Centre line.

The number of Miras in the short period group is small, 17 stars (omitting S Car for reasons mentioned above). However, the lifetime of a Mira is short ($`2\times 10^5`$ years). Thus these stars are tracers of a much larger population. This can be seen from the metal-rich globular clusters, where one such cluster may contain just one Mira. It would be particularly interesting to identify other objects associated with the 145 to 200 day Miras. The globular cluster results suggest that such objects would have metallicities in the range $`0.8[Fe/H]1.3`$ and the age of globular clusters in this metallicity range. However, the situation may be complex. Minniti et al. (1997) suggest that the RR Lyrae variables in the Bulge do not belong to a bar-like distribution, despite the fact that the ones in the NGC 6522 field, at least, have $`[Fe/H]1`$ (Walker & Terndrup 1991). It may be that choosing Miras by period allows one to sort out a more homogeneous population than is possible in other ways.

So far as microlensing is concerned this work would seem relevant for at least two reasons.

1. It gives hope of defining further the composition and general character of the Galactic bar.
2. It warns that one should not necessarily consider Galactic structure and kinematics as “knowns” into which microlensing studies must be fitted. There may well be surprises; and microlensing data may make an important contribution to uncovering them.

###### Acknowledgements.

This paper is based on observations made with the Hipparcos satellite and at the South African Astronomical Observatory (SAAO).


text:
# Polyelectrolyte gels in poor solvent: Equilibrium and non equilibrium elasticity

## I Introduction

Polyelectrolyte gels can absorb a considerable amount of water and swell up to 1000 times. This remarkable property can be monitored by the salt concentration and makes polyelectrolyte gels good candidates for superadsorbent materials. From a fundamental point of view, polyelectrolyte gels are less understood than neutral gels. In some respects, a neutral gel at swelling equilibrium is similar to a semi-dilute polymer solution at the overlap concentration, the mesh size in the gel grows up to the point where the meshes marginally overlap, this is the so-called $`c^{}`$ description. Neutral gels thus benefit to a large extent from the powerful scaling description developed for single neutral chains and solutions. The situation is, however, complicated by trapped entanglements that can, depending on the preparation process, severely restrict swelling and by heterogeneities (fluctuations in the crosslink density) that can be revealed by swelling .

There is no such strong analogy between polyelectrolyte gels and solutions, mainly because of the long ranged electrostatic interaction and the presence of counterions. Counterions are free in the solution whilst they are trapped in the equilibrium gel (immersed in a large excess of salt-free water) to ensure macroscopic electroneutrality. The osmotic pressure of the trapped counterions, with no counterpart in solution, turns out to be the main reason for swelling. Though it originates from the electrostatic interaction between polymer and counterions this pressure is not electrostatic in nature, it does not depend on the strength of the electrostatic interaction, say on the Bjerrum length. The forces responsible for the swelling may be identified as the electrostatic forces acting on the gel boundary where the counterion concentration vanishes. As often in these matters, the thermodynamic point of view, where the contribution from the narrow non neutral boundary is negligible, is more convenient.

Neutral gels also appear simple in the sense that the chain interactions do not play a major role in the elastic properties. The elastic modulus for unentangled networks for example can be simply determined from the density of crosslinks. Even very refined theories  show clearly that the modulus for untentangled networks is given by $`G=ck_\mathrm{B}T/N`$, where $`c`$ is the monomer density and $`N`$ is the number of monomers of the chain strand between two crosslinks. This simple result is quite independent of solvent quality; the excluded volume interactions influence mainly the bulk properties such as the compression modulus, rather than the shear elastic effects; the shear experiments probe, to first order, the strand elasticity and the shear modulus is proportional to the elastic energy density stored in the gel. The elasticity of polyelectrolyte gels is less simple because the gel structure is not simply that of the polyelectrolyte solution. Due to the stretching by the counterion pressure the shear modulus is actually higher than what would be simply given by the crosslink density; it is further expected to explicitly depend on the charge density, the quality of the solvent and the salt concentration.

In the present paper we are concerned with the behavior of polyelectrolyte gels in poor solvent. Khokhlov and coworkers have studied in details the behavior of such ”responsive” gels on the basis of the Flory-Rehner theory, i.e., the additivity of the different parts of the free energy. The starting point of these models is the Gaussian elastic free energy of the neutral networks to which all the other contributions such as, the entropy of the counterions, the electrostatic contributions or the solvent-monomer interactions are added .

The aim of the present paper is different. We study the elastic modulus of the polyelectrolyte gel in a poor solvent using the recent work on the conformation of polyelectrolyte chains. Polyelectrolytes in poor solvents are predicted to form pearl necklace structures . These structures result from a Rayleigh instability of the collapsed globules if their radius exceeds a certain size. In a previous paper we have studied the stretching of these pearl necklace structures . We calculated stress-strain relationships in various situations. For chains with many pearls, we found a continuous stress-strain relation. If the necklace has only a few pearls, it stretches discontinuously and the pearls are dissolved one by one.

The model used in the paper below consists of weakly charged crosslinked chains, with $`N`$ monomers of size $`b`$, a fraction $`f`$ of which are charged. The interaction strength of the Coulomb potential $`V(𝐫)=k_\mathrm{B}Tl_\mathrm{B}/r`$ is expressed by the Bjerrum length $`l_\mathrm{B}=e^2/(4\pi ϵk_\mathrm{B}T)`$; $`ϵ`$ is the dielectric constant of the solvent.

The paper is organized as follows: We first summarize the known results on polyelectrolyte gels in $`\theta `$-solvent and on necklaces . We then discuss equilibrium gels and nonequilibrium gels in section III and finally present our conclusions.

## II Polyelectrolyte gels in $`\theta `$-solvent and necklace conformations

In view of the following considerations it is useful to review first the main points of the present state of the art, via simplest scaling arguments. These include the theory of polyelectrolyte gels in $`\theta `$-solvent, the illustration of the role of the counterions, and the necklace description of polyelectrolytes in poor solvents. For the latter we are going to discuss single chain elasticity and semi-dilute solutions.

### A Polyelectrolyte gels in $`\theta `$-solvent

The size of an isolated polyelectrolyte chain in $`\theta `$-solvent can be estimated from the balance between entropy and electrostatic energy. A weakly charged chain in a $`\theta `$-solvent stretches provided that the fraction of charged monomers $`f`$ is large enough, $`f>N^{3/4}`$. The conformation remains locally isotropic at scales smaller than the electrostatic blob size

$$R_\mathrm{c}=b(l_\mathrm{B}f^2/a)^{1/3}$$
(1)
and the chain is extended at larger length scales. This leads to the usual blob model  with an overall radius : $`R=bN\left(\frac{l_\mathrm{B}f^2}{b}\right)^{1/3}`$ In solution, the chains overlap at concentrations $`c`$ larger than the overlap concentration $`c^{}=N/R^3`$ . In a semi-dilute solution ($`c>c^{}`$) the chains make a random walk of step length $`\xi _\mathrm{c}=\left(l_\mathrm{B}f^2/a\right)^{1/6}(ca^3)^{1/2}b`$ with an overall size

$$R_{1/2}=N^{1/2}\left(l_\mathrm{B}f^2/b\right)^{1/12}(cb^3)^{1/4}b$$
(2)
that becomes equal to the gaussian radius at high concentration when the electrostatic blobs begin to overlap at concentrations $`cc^{}f^{2/3}`$.

The swelling equilibrium of a charged gel results from a balance between the counterion contribution to the osmotic pressure $`\mathrm{\Pi }_{\mathrm{counter}}=k_\mathrm{B}Tfc`$ and the elastic contribution $`\mathrm{\Pi }_{\mathrm{elastic}}=k_\mathrm{B}TcR^2/N^2b^2`$, the equilibrium meshsize and concentration are:

$$R=Nf^{1/2}bc=N^2f^{3/2}b^3$$
(3)
where the so-called $`c^{}`$ theorem $`c=N/R^3`$ which assumes disentangled meshes is used. The meshsize $`R`$ is larger than the size of an isolated chain comprising $`N`$ monomers and does not depend on the Bjerrum length. The contribution of the direct electrostatic interaction to the energy density $`F_\mathrm{e}=k_\mathrm{B}Tl_\mathrm{B}f^2N^2/R^4`$ can be checked to be negligible.

The addition of salt in water at low concentrations $`n_\mathrm{s}<cf`$ does not affect the swelling. At higher salt concentrations, however, the gel is expected to deswell, as the small ion contribution to the osmotic pressure in the gel approaches the osmotic pressure in the surrounding bulk solution. The ion partitioning is fixed by a Donnan equilibrium , that fixes the inner salt concentration $`n_{si}`$ and to the reduced small ion pressure difference $`\mathrm{\Delta }\mathrm{\Pi }_{\mathrm{ion}}=c^2f^2/(4n_\mathrm{s})`$. This pressure difference is compensated by the elastic contribution of the chains in an equilibrium gel and the concentration is

$$c=N^{4/5}f^{6/5}(n_\mathrm{s}b^3)^{3/5}b^3$$
(4)
It is easily checked that the direct electrostatic energy density remains negligible as long as the mesh is weakly screened $`\kappa R<1`$. In the strong screening limit $`\kappa R>1`$ the interactions become short-range. If we impose a scaling form $`F_\mathrm{e}=k_\mathrm{B}T(l_\mathrm{B}f^2N^2/R^4)\mathrm{f}(\kappa R)`$ and require that the free energy density is quadratic in the polymer charge density we obtain $`\mathrm{f}(x)x^2`$ for the scaling function and the electrostatic free energy $`F_\mathrm{e}`$ is of the same order as $`\mathrm{\Delta }\mathrm{\Pi }_{\mathrm{ion}}`$. Equation(4) is thus expected to hold at the scaling level.

### B Necklace structure and elasticity

In a poor solvent a long neutral polymer chain collapses into a dense globule. The poor solvent conditions are characterized by a negative excluded volume $`v=\tau b^3`$, where $`\tau =(\theta T)/\theta `$ is the relative distance to the Flory compensation temperature $`\theta `$. In the absence of electrostatic interactions, the chain collapses into a dense globule that may be viewed as a small region of dense polymer phase at co-existence; the finite size of the globule is associated with an extra energy penalty due to the polymer-water surface tension. The balance of the osmotic pressure between the dense and dilute phases (the latter at almost vanishing concentration) yields the concentration inside the globule $`c=\tau /b^3`$. Alternatively, the dense phase may be described by a close packing of thermal blobs of size $`\xi _t=b/\tau `$, containing $`g=1/\tau ^2`$ monomers. A globule of $`N`$ monomers has then a radius $`R=b\left(N/\tau \right)^{1/3}`$. The corresponding surface tension is $`\gamma =k_\mathrm{B}T/\xi _t^2\tau ^2/b^2`$.

In a poor solvent, a long polyelectrolyte chain adopts an elongated shape determined by a balance between surface tension and electrostatic self-energy. The first model was proposed by Kokhlov who optimized the free energy of a cylindrical globule to obtain the transverse radius given by the electrostatic blob radius $`R_\mathrm{c}`$ eq.(1), and the length $`L`$ imposed by the globule volume $`V=Nb^3/\tau LR_c^2`$. At least at the scaling level, this description holds as long as there is no additional length scale. Recently it was suggested that the charged globule is subject to a Rayleigh instability and splits into connected droplets, the pearls, with a size of order $`R_\mathrm{c}`$. The polymer strands separating the pearls are stretched due to the electrostatic repulsion. However equilibrium between pearls and strands imposes the strand tension. In a first approximation, the strand tension is the equilibrium tension $`k_\mathrm{B}T\tau /b`$ when a polymer strand is pulled out of the dense equilibrium phase; the strands can then be viewed as linear arrays of Gaussian blobs of size $`\xi _t=b/\tau `$. Noting that, up to a logarithmic factor, the repulsion between two half-necklaces is the same as the interaction between two adjacent pearls, we find the strand length $`l`$ and the total chain length $`L`$.

$$l=b(b\tau /l_\mathrm{B}f^2)^{1/2}L=bN(l_\mathrm{B}f^2/b)^{1/2}\tau ^{1/2},$$
(5)
where we use the pearl mass $`m=\tau b/(l_\mathrm{B}f^2)`$ associated with the radius $`R_\mathrm{c}`$.

At low charge fraction or/and very poor solvent, the pearl surface potential calculated in the previous picture is larger than the thermal energy. The counterions are regulated by the pearls, this happens for $`\tau >\alpha ^3(bf/l_\mathrm{B})^{1/3}`$ with $`\alpha k_\mathrm{B}T`$ the counterion chemical potential. The condensation of counterions inside the globules actually suppresses the Rayleigh instability and a large globule is stable against the necklace.

Added salt at concentration $`n_\mathrm{s}`$, screens the electrostatic interactions, the screening length $`\kappa ^1`$ is related to the salt concentration through $`\kappa ^2=8\pi ł_\mathrm{B}n_\mathrm{s}`$. The necklace structure itself is only marginally affected as long as the interaction between adjacent pearls is not screened $`\kappa l<1`$. For higher salt concentrations however the distance between nearest pearls along the chain is almost fixed at the screening length.

In a salt-free solution, necklaces overlap above the concentration $`c^{}=N/R^3`$ and a single chain can be viewed as a random walk with a step length $`\xi _\mathrm{c}=(l_\mathrm{B}f^2/b)^{1/4}\tau ^{1/4}(cb^3)^{1/2}b`$, independent of the overall chain size, and a radius:

$$R=N^{1/2}(l_\mathrm{B}f^2/b)^{1/8}(cb^3)^{1/4}\tau ^{1/8}b$$
(6)
This description holds at moderate concentration as long as there are several pearls per correlation length $`\xi _\mathrm{c}`$, at higher concentration there remains one pearl per correlation length; we not consider this regime any longer.

When the single necklace is submitted to an external pulling force $`\phi `$, pearls are progressively dissolved and converted into strand when the force is increased. When there are only a few pearls, each pearl dissolution translates into a length jump under imposed force or a force drop under imposed length. For a higher number of pearls however the pearl number fluctuates by more than one unit, single pearl features are washed out and the force curve is continuous. This is the regime of interest for the gel system where elastically active paths can be quite long. The chain length can be related to the pulling force by a force balance on the half-necklace similar to the one leading to eq.(5) including the external force $`\phi `$:

$$L(\phi )=L(0)\left(1\frac{\phi b}{k_\mathrm{B}T\tau }\right)^{1/2}$$
(7)
(where all finite size corrections are neglected). The elastic energy stored in the chain is dominated by the stretched strands (provided $`L>L(0)`$) and is about $`k_\mathrm{B}T`$ per tensile-blob of radius $`\xi _t`$. This simple argument fails when most of the pearls are dissolved, i.e. for $`L>Nb\tau `$.The maximum length (when all pearls are used up) which represents an upper bound for the validity of the continuous model.

## III Polyelectrolyte gels in poor solvent

### A Equilibrium gels

We want here to include the details the pearl-necklace structure for the elastic chains of the gel; this cannot be done using complete calculations as for the $`\theta `$-gel. However, useful predictions can be made at the scaling level. We propose two different models : the $`c^{}`$-gel where the meshes are disentangled at swelling equilibrium and for which the preparation state is irrelevant, and the affine deformation model where the equilibrium swelling depends on the preparation state. The affine model considers that the deformation of the gel with respect to the preparation state is affine at all length scales; this imposes the constraint that the polymer content of a mesh volume remains constant upon swelling $`cL_m^3/N=c_\mathrm{p}L_\mathrm{p}^3/N=\alpha `$ where index $`p`$ refers to the preparation state. This constraint can be understood as imposing that initially entangled meshes do not disentangle. The special (limiting) case $`\alpha =1`$ corresponds to the $`c^{}`$-gel, obtained here by reacting all temporary contact points between different chains of a semi-dilute solution during the crosslinking process. Higher values of $`\alpha `$ are obtained if the gel is prepared in a semidilute or a dense solution when only a small fraction of the temporary contact points are crosslinked.

#### 1 $`c^{}`$-model

The $`c^{}`$-model is based on the Katchalsky picture . The equilibrium gel is immersed in a solvent reservoir that imposes its osmotic pressure. In the absence of salt (actually if the salt concentration $`n_s`$ is such that $`n_\mathrm{s}<cf`$), the osmotic pressure in the gel vanishes (it is much lower than the ideal gas counterion pressure $`k_\mathrm{B}Tcf`$). The contribution arising from the counterions and from the strand elasticity thus (almost) cancel.

$$\mathrm{\Pi }=\mathrm{\Pi }_{\mathrm{elastic}}+\mathrm{\Pi }_{\mathrm{count}}=0$$
(8)
The counterion contribution to the osmotic pressure, similar to that for the $`\theta `$-gel is dominated by the entropic term $`\mathrm{\Pi }_{\mathrm{counter}}=k_\mathrm{B}Tcf`$. The contribution arising from the strand elasticity can be derived from the elastic energy density that is proportional to the density of tensile-blobs :

$$F_{elastic}/k_\mathrm{B}T=(c/N)(N/m)(l/\xi )=\tau /(L_\mathrm{m}^2b)$$
(9)
it varies as the $`2/3`$ power of concentration. As for the $`\theta `$-gel, we assume that disentangled meshes are at equilibrium, $`c=N/L_\mathrm{m}^3`$. The contribution of the elastic energy to the osmotic pressure $`cdF_{\mathrm{elastic}}/dcF_{\mathrm{elastic}}`$ is negative, and only differs from $`F_{\mathrm{elastic}}`$ by a factor $`1/3`$. At the scaling level we thus write

$$\mathrm{\Pi }_{\mathrm{elastic}}/k_\mathrm{B}T=(c/N)(N/m)(L_\mathrm{m}/\xi )=\tau /(L_\mathrm{m}^2b)$$
(10)
This is compensated by the counterion pressure $`\mathrm{\Pi }_{\mathrm{count}}=k_\mathrm{B}Tcf`$ at the gel equilibrium concentration :

$$c=\frac{\tau ^3}{N^2f^3b^3}$$
(11)
where the meshsize is

$$L_\mathrm{m}=bfN/\tau $$
(12)
It is easily checked that the mesh is stretched with respect to a corresponding free necklace comprising $`N`$ monomers, and as for a $`\theta `$-gel, the structure of the gel in a poor solvent is not that of a $`c^{}`$-solution. On the other hand our description of the stretched necklace requires the necklace to be shorter than a string of thermal blobs as the pearl-necklace description remains valid only if most monomers belong to pearls: $`L<N\tau b`$. Eq.(12) thus holds for $`\tau >f^{1/2}`$; closer to the $`\theta `$-point, the solvent quality becomes irrelevant and the structure is that of a $`\theta `$-gel. It can actually be checked that the concentration smoothly crosses over with that of the $`\theta `$-gel when $`\tau =f^{1/2}`$.

For very poor solvents we know from the single necklace study that the counterions are no longer free and condense on the pearls. The gel then macroscopically collapses, this happens for

$$\frac{l_\mathrm{B}\tau ^3}{bf}>\left(\mathrm{log}\frac{N^2f^2b^3}{\tau ^3}\right)^3$$
(13)
where the logarithm accounts for the counterion entropy at the onset of condensation. Qualitatively, the gel is collapsed for $`\tau >f^{1/3}`$. In summary, the poor solvent conditions are irrelevant close to the $`\theta `$-point as long as $`\tau <f^{1/2}`$, the gel shrinks continuously for intermediate solvents $`f^{1/2}<\tau <f^{1/3}`$ and undergoes collapse at $`\tau f^{1/3}`$ due to counterion condensation. This behavior is represented in FIG (2).

The compression modulus and the shear modulus of the equilibrium gel, are proportional to the counterion free energy density and to the elastic energy density respectively, and are of the same order at equilibrium:

$$GK\mathrm{\Pi }_{\mathrm{count}}=\frac{k_\mathrm{B}T}{N^2f^2}\left(\frac{\tau }{b}\right)^3$$
(14)
When monovalent salt is added to the bulk solution, a Donnan equilibrium fixes the inner salt concentration and the small-ion excess pressure. For salt concentrations $`n_\mathrm{s}`$ smaller than the bound charge concentration $`cf`$, salt is irrelevant, and the small ion excess pressure remains equal to $`k_\mathrm{B}Tcf`$. For $`n_\mathrm{s}>cf`$ in contrast, the small ions contribute a smaller inner pressure excess $`\mathrm{\Delta }\mathrm{\Pi }_{\mathrm{ion}}=k_\mathrm{B}Tc^2f^2/(4n_\mathrm{s})`$.

This osmotic pressure difference is compensated by the elastic contribution given by eq.(10) at the equilibrium meshsize and concentration:

$$L_\mathrm{m}=\frac{N^{1/2}f^{1/2}}{\tau ^{1/4}(n_\mathrm{s}b^3)^{1/4}}bc=\frac{\tau ^{3/4}(n_\mathrm{s}b^3)^{3/4}}{N^{1/2}f^{3/2}b^3}$$
(15)
At low salt concentration $`n_\mathrm{s}b^3<N^2f^{1/2}(\tau /f^{1/2})^3`$ these results crossover the salt-free results eqs.(11,12), line (1) in the diagram FIG (2). Close to the $`\theta `$-point, $`n_\mathrm{s}b^3<N^2f^{1/2}(\tau /f^{1/2})^5`$ the negative excluded volume becomes irrelevant and the results crossover to those obtained for a $`\theta `$-gel, line (2) in the diagram. The weak screening regime described by eq.(15) is also bound at high salt concentrations : the screening length due to the added salt becomes shorter than the meshsize for $`n_\mathrm{s}b^3=\tau b^2l_\mathrm{B}^2N^2f^2`$; at higher salt concentration the salt imposes the screening length, line (3) in the diagram. On this line, the electrostatic interaction between polymer strands in the mesh contributes as much as the small-ions to the excess pressure. The crosslinks are thus only marginally relevant. One easily checks that the mesh size is indeed equal to the size of a free strand of mass $`N`$. This means that, at the scaling level, there is no difference between the structure of the gel and that of a semi-dilute solution above line (3). Lines (2) and (3) meet at point $`B(\tau =f^{2/3}(l_\mathrm{B}/b)^{1/3},n_\mathrm{s}b^3=f^{4/3}N^2(b/l_\mathrm{B})^{5/3})`$. For poorer solvents $`\tau >f^{1/3}`$, the small ions condense on the pearls and the gel collapses, line (4) in the diagram; this line has not been drawn precisely, in particular the condensation threshold weakly depends upon the salt concentration $`n_\mathrm{s}`$.

#### 2 The affine deformation model

The condition that the meshes only marginally overlap at swelling equilibrium is now released and replaced by the more general condition of affine deformation $`cL_m^3/N=\alpha `$ where the constant $`\alpha `$ is determined in the preparation state. As explained earlier, the physical content of this assumption is quite opposite to that of the $`c^{}`$-model. The $`c^{}`$-model is nonetheless recovered in the special case $`\alpha =1`$ corresponding to a semi-dilute solution where all temporary contact points between different chains have been made permanent by the crosslinking reaction. In practice $`\alpha >1`$ reaches its highest value $`N^{1/2}`$ for crosslinking in the neutral melt. The elastic energy density now depends on the preparation state through the parameter $`\alpha `$.

$$F_{el}=(c/N)(L_\mathrm{m}/\xi _{th})=\frac{\tau \alpha }{bL^2}=\frac{\tau }{b}\left(\frac{c}{N}\right)^{2/3}\alpha ^{1/3}$$
(16)
The equilibrium meshsize and concentrations have to be modified accordingly and typically differ from those in the $`c^{}`$-model by a power of $`\alpha `$.

For weakly entangled preparation states $`1<\alpha <f^{1/3}`$ the topology of the diagram of states presented in FIG.(2) is preserved. The salt concentration on dividing lines (1), (1$`\theta `$) and (2) is increased by a factor $`\alpha `$ whilst it is decreased by a factor $`1/\alpha `$ on line (3) and by a factor $`\alpha ^{2/3}`$ on line (3$`\theta `$). The points where these lines cross are shifted as shown in FIG.(3).

As a consequence the parameter range of the intermediate regimes $`\theta +S`$ and $`\tau +S`$ is reduced. For somewhat more entangled preparation states $`f^{1/3}<\alpha <f^{1/2}`$, lines (1) and (3) cross: For poor enough solvent $`\tau >1/\alpha `$, strong screening occurs before the Donnan effect is relevant as shown on FIG(4).

In this poor solvent regime a new dividing line (3$`\tau `$) directly separates the strongly screened region from the $`\tau `$ region where salt is irrelevant. On this line $`n_\mathrm{s}b^3=b\tau ^2/(l_\mathrm{B}N^2f^2)`$. There is direct contact between the strong screening and the salt irrelevant region over the whole $`\tau `$-range for $`\alpha =f^{1/2}`$ (see figure 3). For higher values of $`\alpha `$ the remaining lines (3) and (3$`\theta `$) no longer shift with $`\alpha `$.

The expression of the gel equilibrium concentration and elastic modulus in the various regimes are summarized in tableI. Let us stress that for $`\alpha >1`$, the gel structure on lines (3) is not that of the solution at the same concentration.

### B Shear modulus of non equilibrium gels

When a small amount of water is added to the gel, the added amount imposes the volume and the gel is not at swelling equilibrium. In this case the preparation state is of importance. Most of the (elastic) properties depend on the preparation state . Rather sophisticated theories have been developed for gels that have not reached swelling equilibrium, but here we argue on a much simpler level. The counterion contribution to the osmotic pressure is given as above by the ideal gas law. We can therefore conclude that the compression modulus is of the order of $`k_\mathrm{B}Tcf`$ with $`c`$ the actual gel concentration. The shear modulus on the other hand, depends on both the preparation state (crosslink configuration) and the swelling state. For a sake of simplicity it is assumed that the gel deforms affinely, and the preparation state is assumed to be an undeformed isotropic state. Thus the measure of the affine deformation with respect to the preparation size is determined by the ratio $`\lambda =(c_\mathrm{p}/c)^{1/3}`$. where $`c_\mathrm{p}`$ and $`c`$ are the initial and final gel concentrations, respectively. To estimate the elastic energy density stored in the gel we take the semidilute solution at the same concentration $`c`$ where the chains are stress free as a reference state. At high enough concentration, the elastic chains of the gel have gaussian statistics at large length scale and behave as stretched Gaussian chains of blobs of radius $`\lambda R_p`$. The modulus can be written as

$$G=\frac{k_\mathrm{B}Tc}{N}\left(\frac{\lambda R_\mathrm{p}}{R}\right)^2$$
(17)
where $`R`$ is the equilibrium radius of a chain in a semidilute solution at the same concentration. For moderate concentrations, $`c<(l_\mathrm{B}/b)^{1/2}f\tau ^{1/2}b^3`$, several pearls are found in one correlation length $`\xi _\mathrm{c}=(l_\mathrm{B}f^2/b)^{1/4}\tau ^{1/4}(cb^3)^{1/2}`$ and the chain radius is given by equation 6. We find then for the shear modulus

$$Gk_\mathrm{B}T\left(\frac{\tau }{f^2}\right)^{1/4}c^{5/6}$$
(18)
where the dependence upon the preparation state has been dropped. This analysis is valid in the weak stretching limit where the gel strand is not stretched. For more swollen gels, $`c<c_1`$, the radius of an elastic chain $`R_c=R_\mathrm{p}(c/c_\mathrm{p})^{1/3}`$, becomes larger than the radius of an isolated necklace given in eq.(5). The crossover concentration $`c_1`$ depends explicitly on the preparation state. For $`c<c_1`$, the shear modulus is dominated by the overstretched elastic chains and one should use the elastic free energy of equation 16

$$G=k_\mathrm{B}Tc^{2/3}c_\mathrm{p}^{1/3}\tau R_\mathrm{p}/(Nb)$$
(19)
Formally the modulus exhibits a sharp jump at the concentration $`c_1`$ within the scaling picture. Around $`c_1`$ a more careful consideration of the single chain elasticity is needed. This could be done using a non gaussian elasticity model such as the Langevin model for the chain of blobs, we are reluctant to include it in our rough scaling analysis.

Within the affine deformation model the equilibrium swelling concentration given by $`G=cf`$ depends explicitly on the preparation state, i.e.,

$$c_{eq}=c_\mathrm{p}R_\mathrm{p}^3\left(\frac{\tau }{Nf}\right)^3$$
(20)
in contrast to the $`c^{}`$-model. The two approaches only agree if the initial semi-dilute solution is cross-linked at all the contact points, where the $`c^{}`$-theorem is indeed expected to apply. In that case, $`c_\mathrm{p}R_\mathrm{p}^3=N`$ and eqs.(20,11) coincide. With this preparation state, the shear modulus of a gel that has not reached equilibrium reads:

$$G=k_\mathrm{B}T\left(\frac{c}{N}\right)^{2/3}\frac{\tau }{b}$$
(21)
Our results can be tested experimentally by checking the final charge fraction (or temperature) dependence of the shear modulus for identical preparation states. Such experiments, where the shear modulus was indeed found to decrease with the charge fraction at fixed final density, were performed in the group of J.Candau in Strasbourg .One should note however that the experiments were done with annealed polyelectrolytes by monitoring the final charge with the pH of the solution. Annealed polyelectrolytes, as shown recently , if the solvent quality is too low can undergo a sharp transition from a collapsed state to a stretched state without any intermediate regime of stable necklace. Our results would thus based on the pearl-necklace conformation would only be useful if the solvent is not too poor which might be the case for polyacrylic acid gels.

## IV Conclusion

In this paper we have studied the elastic behavior of polyelectrolyte gels in poor solvent which is markedly different from that of polyelectrolyte gels in good solvent. The study of the single chain behavior already illustrates the different physical behaviors. The stretching of polyelectrolyte chains in a good solvent can be carried out via a simple blob analysis, as the undeformed state of the chains can be described as an extended chain of electrostatic blobs. In a poor solvent, the undeformed chains are composed of pearls and strings. The pearls act as a ”reservoir” of monomers which can be pulled out at relatively small forces. In the continuum limit the force increases gently with the deformation. In the case of a low number of pearls, we expect length jumps that correspond to single pearl unwinding at an imposed force. These jumps are rounded by thermal fluctuations in many instances.

For a polyelectrolyte network, we do not expect any force jumps as the elastic chains carrying the stress are much larger than a single mesh; the continuous description is thus appropriate. A gel at swelling equilibrium in a large excess of $`\theta `$-solvent is known to swell under the effect of the trapped counterion pressure. Most polyelectrolytes however have hydrophobic backbones and the interplay between short range attractions and electrostatics, that leads to the pearl-necklace conformation, is expected to be important as well. Close to the $`\theta `$-point, for $`\tau <f^{1/2}`$, the solvent quality is irrelevant and $`cf^{3/2}`$. Note that the crossover between the $`\theta `$-solvent and the poor solvent behavior is different for the equilibrium gel and for isolated chains ($`\tau >f^{2/3}`$). For intermediate solvent quality $`f^{1/2}<\tau <f^{1/3}`$ the gel deswells and $`c(\tau /f)^3`$. When the solvent is very poor $`\tau f^{1/3}`$, the counterions condense on the pearls and the gel collapses.

When salt is added to the bulk solution, the gel shrinks. The various swelling regimes are represented on the diagram of Fig.(2). If the salt concentration is smaller than counterion concentration in the gel, the effect of salt is negligible. At higher salt concentration, the concentration in the gel increases as $`cn_\mathrm{s}^{3/4}`$. At salt concentrations higher than $`\tau /(N^2f^2)`$, the gel has shrunk to the extent that junction points no longer play any role, the properties of the gel are similar to those of a semi-dilute polyelectrolyte solution (the details of this description remain controversial). In all the equilibrium swelling regimes the compression and shear moduli are roughly proportional to each other and to the small-ion excess pressure. The modulus is independent of the strength of the electrostatic interactions, say of the Bjerrum length.

In the more realistic affine deformation model, the equilibrium swelling depends upon the preparation state through the constant $`\alpha =cL_\mathrm{m}^3/N`$ measuring the degree of mesh entanglement. The affine deformation constraint imposes that $`\alpha `$ remains constant during the swelling i.e., that the polymer content of a mesh volume remains constant and that meshes do not disentangle. The diagram of states is qualitatively similar to that obtained form the $`c^{}`$-model, meshsizes, the scaling laws for the equilibrium concentration and the elastic modulus are shifted by a power of $`\alpha `$ with respect to that obtained for the $`c^{}`$-gel that corresponds to $`\alpha =1`$. The results are summarized in tableI. For higher degrees of entanglement, $`f^{1/3}<\alpha <f^{1/2}`$, the diagram qualitatively changes: For poor enough solvent $`\tau >1/\alpha `$ there is a direct separation line between the salt irrelevant regime and the strong screening regime. For $`\alpha >f^{1/3}`$ this happens for any solvent quality.

A small amount of added solvent controls the gel volume, and the swelling remains below the equilibrium swelling. The compression and shear moduli are no longer proportional to each other, even can show opposite variations with the charge faction. Similar results where obtained previously for the $`\theta `$-gels and observed in experiments. We find a shear modulus $`G(\tau /f^2)^{1/4}c^{5/6}`$, that is shear softening upon charge increase when all other parameters are kept constant.

Our picture relies on a scaling description that allows us to account for the essential physics. Some additional subtle effects may arise from the long range character of the electrostatic interaction and its variation with the salt concentration. For gels in a good solvent  it has been suggested recently  that the elastic modulus could depend on the strength of the interaction and the range of the Debye - Hückel potential at least in some parameter range. Our rough approach would not capture those subtle effects.

### Acknowledgments

TAV acknowledges the financial support of the LEA. The warm hospitality of the Institut Charles Sadron is gratefully appreciated.


text:
# Comment on “Local dimer-adatom stacking fault structures from 3×3 to 13×13 along Si(111)-7×7 domain boundaries”

## Abstract

Zhao et al. \[Phys. Rev. B58, 13824 (1998)\] depicted several atomic structures of domain boundaries on a Si(111) surface and criticized the article by the present author and the co-workers. I will point out that their criticism is incorrect and their structure models have no consistency.

Zhao and co-workers referred to the work by the present author and the co-workers  in the former part of their article  and stated as if we merely tried to make a scheme to classify the domain boundaries (DBs) appearing on Si(111) surfaces with the $`7\times 7`$ dimer-adatom-stacking-fault (DAS) reconstruction . However, what we actually showed in Ref.  was that any structures of a reconstructed Si(111) surface including those in the DBs could be consistently accounted for by the atomic bondings which are already contained in the DAS structure. And this is why it was necessary to classify all the possible patterns of the DBs, followed by the structural analyses of the 14 patterns among the 27 possible ones. Then, we assumed the dimer-adatom interaction  and the reduction in the number of dangling bonds as the guiding principles to show that the atomic structures can be consistently determined with the use of the constituents of the DAS structure only.

In Ref. , Zhao et al. criticized Fig. 6(b) of Ref.  by stating that we had given too many dangling bonds and, accordingly, they drew Fig. 1 in Ref.  by adding one dimer and one rest atom in the DB in Fig. 6(b) of Ref. . Yet, they did not add such dimers to the DB in Fig. 2, which is the structure model of their own scanning tunneling microscopy (STM) image. In contrast to one of our assumptions, the dimer they added in Fig. 1 is not accompanied by any adatoms. It is clear from this that their structural analyses have no consistency by themselves because they added a dimer only to Fig. 1 and not to Fig. 2. Also, they added the rest atom to the position where its backbonds must suffer from the strong distortion to increase a stress energy. Although it is necessary to carry out STM observations with negative sample bias voltage $`V_s`$ to identify the locations of rest atoms , all of their STM images were obtained at $`V_s=2.0`$ V only . Based only on such STM images, nobody can discuss whether rest atoms are really present at surface sites which look dark in the empty-state images. For this reason, Fig. 1 of Ref.  is dubious.

In addition, although they stated that Fig. 2 of Ref.  was beyond our classification scheme, one can find the corresponding shift in the location of the $`(5,5)`$ pattern in Fig. 3 of Ref.  and, hence, their criticism is simply incorrect. Actually, based on our assumptions and with the help of a mirror symmetry as we discussed in Ref. , it is easy to deduce the atomic structures of the $`(n,n)`$ patterns for any integer values of $`n`$ ($`0n6`$), including the structure model they depicted in Fig. 2; it merely confirmed the validity of our scheme and the assumptions we adopted in Ref. .

According to their explanation, Fig. 4 was depicted based on the STM image of Fig. 6(a) in Ref. , in which the two $`3\times 3`$ structures are found apart from each other with their adjacent sites being either covered by the contaminations or made of the vacancies. Only with such a STM image, and without any filled-state images, again it is impossible to deduce the atomic structure correctly.

Finally, I will point out that a surface structure appearing after annealing depends on the annealing schedule due to an entropic effect , so that their argument in deducing the energy difference between the $`9\times 9`$ structure and the $`5\times 5`$ structure is insufficient to derive any definite conclusions.

In summary, firstly since they did not state any criteria to fulfill systematic analyses, their structure models have no consistency among them. Secondly since they did not observe any filled-state STM images, their results have no reliability on the presence or the location of rest atoms. Thirdly, they criticized the failure of our classification scheme based on the observation of the structure which is easily deduced by our classification scheme and by our assumptions and, hence, their criticism is simply incorrect. And finally they overlooked the importance of the entropic effect which is involved in the determination of atomic structures in an annealing process.