text:
# Electronic structure of the strongly hybridized ferromagnet CeFe2

## I Introduction

The 4$`f`$ states of rare-earth elements in solids usually retain free-ionic properties with a well-defined integer occupation number. However, there are also rare-earth compounds where the hybridization of the 4$`f`$ states with extended band states is important, in which case they may exhibit properties usually only found in actinide systems. In such systems, many unusual phenomena are typically observed, like for instance anomalously low saturation magnetization and Curie temperature $`T_C`$ (e.g., CeFe2), intermediate valence (e.g, SmS), heavy fermion behavior (e.g., YbBiPt), or non-Fermi-liquid behavior (e.g., CeCu6-xAux). Even more surprisingly, simultaneous magnetic ordering and superconductivity has been observed (e.g., CeCu2Si2). The superconductivity is unconventional, i.e., the order parameter suggests a d-wave superconducting state, as opposed to the conventional s-wave state. Finding a proper theoretical description of the 4$`f`$ states in these compounds remains one of the major problems in condensed matter physics.

As for the photoemission spectroscopy (PES) studies of Ce compounds, it is widely believed that the spectra are well described by the single-impurity Anderson model (SIAM).  Recently, however, it has been argued that systems, in which the Ce 4$`f`$ states hybridize strongly with the other valence electrons, calculations based on density functional theory (DFT) may give an equally good, or even better description of the photo-emission spectra than the SIAM analysis, provided that surface effects are properly taken into account in the analysis.  However, one should bear in mind that calculations based on DFT are not strictly applicable for excited state properties, instead the ground state properties, such as the magnetic moments of the ground state, which are typically the focus of these calculations. Nevertheless, the electronic structure given from such calculations are often compared with photoemission data and good agreement between experiment and calculations is frequently observed. In the limit of complete screening of the excited state, one would expect ground-state density functional calculations to be able to describe the spectra well.

CeFe2 is thought to belong to a class of strongly hybridized systems. This compound shows ferromagnetism below $`T_C`$ = 230 K with a saturation magnetization of 2.30 $`\mu _B`$/f.u. Above $`T_C`$, the magnetic susceptibility follows the Curie-Weiss law with an effective moment of 7.4 $`\mu _B`$/f.u.  If one compares CeFe2 with the other $`R`$Fe2 compounds ($`R`$: rare earth elements), a number of anomalies in its physical properties can be observed. The lattice constant is much smaller than an interpolation using the lattice constants of the other $`R`$Fe2 systems would suggest. Its Curie temperature is anomalously low: the other $`R`$Fe2 compounds have Curie temperatures ranging from 596 K to 796 K.  The saturation magnetization is unusually low compared to the other $`R`$Fe2 compounds (2.93 and 2.90$`\mu _B`$/f.u. for LuFe2 and YFe2, respectively ). Moreover, even if only a small fraction of the Fe atoms are substituted for Al, the ferromagnetic ordering is destroyed, and the system becomes anti-ferromagnetic.  In fact, even in pure CeFe2, recent neutron scattering experiments have revealed strong competition between the ferromagnetic ground state and an antiferromagnetic ground state. Together, these facts suggest that the Ce 4$`f`$ states in CeFe2 hybridize strongly with the other valence electrons, notably the Fe $`3d`$ valence states. This hypothesis is further supported by the X-ray absorption (XAS) experiments by Croft et al.

In this paper, we present high-energy spectroscopic results on CeFe2 including core-level x-ray photoemission (XPS), XAS, Ce 3$`d`$-4$`f`$ and 4$`d`$-4$`f`$ resonant PES, bremsstrahlung isochromat (BIS) and high-resolution ultraviolet photoemission spectroscopy (UPS) in order to elucidate the electronic structure of this system. In the case of Ce compounds with strongly hybridized 4$`f`$ states, it has been pointed out that surface effects are extremely important in the interpretation of the spectra.  Therefore, we have attempted to differentiate the electronic structure of bulk and that of surface for CeFe2 in the analysis of the spectra. As will be further elaborated on in section III C of this paper, electronic-structure calculations with the Ce 4$`f`$ states treated as valence states give a good description of the magnetism in CeFe2 It is of course highly interesting to assess the applicability of the same theory in describing also the photoemission spectra of CeFe2, even though as noted these calculations are not strictly applicable for excited state properties. Very recently, Sekiyama et al.  reported a high-resolution 3$`d`$-4$`f`$ resonant photoemission study of the strongly hybridized system CeRu2 and found that the Ce 4$`f`$ spectra can be explained by band theory. In the following, we first attempt to describe the spectra in the framework of the SIAM and obtain a set of SIAM parameters. In the analysis, surface effects on the spectra are explicitly taken into account. Next, the bulk component of the valence-band spectra is compared with the density of states calculated using density functional theory. All DFT results presented here have been calculated within the local spin-density approximation (LSDA). The use of more recently developed generalized gradient functionals would not, however, alter any of our conclusions.

## II Methods

### A Experiment

Polycrystalline samples of CeFe2 were prepared by arc-melting the pure constituent materials. Subsequently, the samples were annealed at 750 C for a week to obtain single phase samples. Magnetization measurements yielded the same $`T_C`$ as in the literature. The XPS spectra were taken with Mg $`K\alpha `$ radiation ($`h\nu `$ = 1253.6 eV) using a double-pass cylindrical-mirror analyzer, and the BIS spectra were obtained using a Pierce-type electron-gun and a quartz crystal monochromator which was set at $`h\nu `$ = 1486.6 eV. The Ce 4$`d`$-4$`f`$ resonant PES measurements were done at beam-line BL-2 of SOR-RING, Institute for Solid State Physics, University of Tokyo. The Ce 3$`d`$-4$`f`$ resonant PES and Ce 3$`d`$ XAS data were taken at beam-line BL-2B of Photon Factory, High Energy Accelerator Research Organization. Photoelectrons were collected using a double-pass cylindrical-mirror analyzer in the resonant PES measurements. The XAS spectra were obtained by measuring the total electron yield using an electron multiplier placed near the sample. All measurements were done in the range 50-80 K, i.e., below the Curie temperature. In the case of the Ce 4$`d`$-4$`f`$ resonant PES, additional measurements at room temperature, i.e, above the Curie temperature, were performed. The total energy resolution was $``$1.0 eV for XPS and BIS, $``$0.5 eV for Ce 4$`d`$-4$`f`$ resonant PES, $``$0.5 eV for XAS and $``$1.0 eV for Ce 3$`d`$-4$`f`$ resonant PES. The high-resolution UPS measurements were done around 17 K using a hemispherical analyzer and the He I ($`h\nu `$ = 21.2 eV) and He II ($`h\nu `$ = 40.8 eV) resonance lines. The energy resolution was $``$25 meV for both photon energies. The binding energies were calibrated using Au evaporated on the samples. For XAS and Ce 3$`d`$-4$`f`$ resonant PES, the photon energies were calibrated using the Cu 2$`p`$ edge of Cu metal and the Co 2$`p`$ peak of LaCoO3. Clean surfaces were obtained by scraping the sample repeatedly, while maintaining the sample under ultra-high vacuum, with a diamond file prior to each measurement. Cleanliness of the surfaces was checked by the absence of O 1$`s`$ and C 1$`s`$ XPS signals from contaminants in the case of the XPS, XAS, BIS and Ce 3$`d`$-4$`f`$ resonant PES measurements. In the case of the high-resolution UPS and Ce 4$`d`$-4$`f`$ resonant PES measurements, cleanliness was checked by the absence of a O 2$`p`$ feature which appears around 6 eV below the Fermi level ($`E_F`$).

### B Single-impurity Anderson Model

The SIAM calculations were made based on the variational 1/$`N_f`$-expansion method developed by Gunnarsson and Schönhammer.  Here, we performed the calculations to the lowest order in 1/$`N_f`$, where $`N_f`$ is the degeneracy of the Ce 4$`f`$ level and was taken to be 14. The $`f^2`$ configuration was also included in the calculation. The energy dependence of the hybridization matrix elements was taken from the off-resonant spectra, which approximately represent the Fe 3$`d`$ partial density of states. The configuration dependence of the hybridization strength was also taken into account, and was chosen to be the same as that obtained for $`\alpha `$-Ce by Gunnarsson and Jepsen.  In the calculations, we divided the band continuum into discrete levels following Kotani et al.  We further assumed that each spectrum was a superposition of two components which represent bulk and surface spectra. The weight of each component was treated as fitting parameters within a range consistent with the universal curve for the mean free path of photoelectrons.

### C Band-Structure Calculation

In the band-structure calculations presented here, we have used the full-potential linear muffin-tin orbital method (FP-LMTO).  In this method, the Kohn-Sham equations  are solved for a general potential without any shape approximation. The local (spin) density approximation (LSDA) in the Hedin-Lundqvist parameterization  was used for the density functional.

In the FP-LMTO method, space is divided into non-overlapping spheres, so-called called muffin-tin spheres,  surrounding each atomic site, and an interstitial region. The basis functions used are energy independent Bloch functions, whose construction is different in the spheres and in the interstitial.

A basis function in the interstitial is defined by the Bloch function of solutions to the spherical Helmholtz equation with nonzero kinetic energy $`\kappa ^2`$, or a linear combination of such solutions for different kinetic energies. The Fourier representation of this basis function is taken from the Fourier series of a function matching the basis in the interstitial region but not inside the spheres, a so-called pseudo-wave function, whose exact shape inside the muffin-tin sphere is of no importance for the final solution as long as it is continuous and differentiable at the sphere boundary and matches the true basis function in the interstitial.

Inside the spheres, where the charge density varies rapidly, the basis functions are Bloch functions of numerical radial functions times spherical harmonics. The radial part of a basis function is constructed from the numerical solutions $`\varphi _L(E_\nu ,r)`$ of the radial Schrödinger equation in a spherical potential at the fixed energy $`E_\nu `$, and their energy derivatives $`\dot{\varphi }_L(E_\nu ,r)`$. Here, the index $`L`$ stands for a collection of quantum numbers: the principal quantum number $`n`$, the orbital quantum number $`l`$, the magnetic quantum number $`m`$, and the kinetic energy $`\kappa ^2`$.

The treatment of the entire basis set within one single energy panel allows all states, including the semi-core states, to hybridize fully with each other. Our method is linear, i.e., the basis functions are constructed by expanding around fixed energies $`E_\nu `$. The expressions for the crystal wave functions in the muffin-tin spheres are matched to the interstitial crystal wave function at the sphere boundaries so that the total crystal wave function becomes continuous and differentiable in all space. In the present calculation, the expansion in spherical harmonics was taken up to $`l=6`$. For Ce, the $`6s`$, $`5p`$, $`6p`$, $`5d`$, and $`4f`$ orbitals were included in the basis set, with $`5p`$ as semi core. For Fe, we included the $`4s`$, $`4p`$ and $`3d`$ orbitals, i.e., no semi-core state was used for Fe. Four $`\kappa ^2`$-values were used in the calculation: -0.6 Ry and -0.1 Ry for the valence states, and -1.5 Ry and -1.0 Ry for the semi-core Ce $`5p`$ states, all with respect to the muffin-tin zero.

Reciprocal space was sampled with what would correspond to 1331 k-points in the full Brillouin zone (BZ) using special k-point sampling methods.  The non-overlapping muffin-tin spheres were chosen as 21% and 17% of the unit-cell lattice constant for Ce and Fe, respectively. With this choice, 36% or the unit cell volume is in the interstitial region and the closest muffin-tin spheres are 3% from touching.

The experimental lattice constant was used in the calculations. Furthermore, the calculations were spin-polarized but the spin-orbit interaction was not included. This latter approximation will be commented on further in conjunction with presenting and discussing the results from the band-structure calculation.

## III results

### A Experiment

Fig. 1 shows the Ce 3$`d`$ core-level XPS and XAS spectra. The XPS line-shape is a typical one for a strongly hybridized Ce compound, consisting of three peaks which correspond to the 3$`d^9`$4$`f^0`$, 3$`d^9`$4$`f^1`$ and 3$`d^9`$4$`f^2`$ final states in each of the $`j`$=3/2 and 5/2 spin-orbit components . In the XAS spectrum, the main peaks are due to the 3$`d^9`$4$`f^2`$ final-state multiplet and the satellite structures $``$5 eV above the main peaks are due to the 3$`d^9`$4$`f^1`$ final states. The rather distinct 4$`f^0`$ peaks in the XPS spectrum and the 4$`f^1`$ structures in the XAS spectra, together with the obscured 3$`d^9`$4$`f^2`$ final-state multiplet structures of the main XAS peaks, indicate strong hybridization of the 4$`f`$ states with the valence band in this system. The XPS spectrum reflects the surface electronic structure because of the rather low kinetic energies of photoelectrons from the Ce 3$`d`$ core level. A detailed analysis of this is given below.

The results of valence-band PES and BIS are shown in Fig. 2. The on- and off-resonance occurs, respectively, at $`h\nu `$ = 121 and 114 eV in the Ce 4$`d`$-4$`f`$ resonant PES and at $`h\nu `$ = 881.1 and 875.4 eV in the Ce 3$`d`$-4$`f`$ resonant PES. Identical spectra in the present resolution have been obtained for Ce 4$`d`$-4$`f`$ resonant PES at room temperature, which is above $`T_C`$ (not shown). We have obtained the Ce 4$`f`$ spectra by subtracting the off-resonance spectra from the on-resonance spectra as shown by solid curves in Fig. 2. As seen from this figure, there is a large difference between the Ce 4$`f`$ spectra obtained from the Ce 4$`d`$-4$`f`$ and Ce 3$`d`$-4$`f`$ resonant PES. While the former has a double-peak structure in the vicinity of -2 eV and near $`E_F`$, the latter is dominated by a single peak near $`E_F`$, implying stronger hybridization in the latter. This can be attributed to the difference in the surface sensitivity of the two spectra due to the different kinetic energies of photoelectrons. This also indicates that the Ce valency at the surface is closer to trivalent than it is in the bulk.

In Fig. 2 we also show the XPS spectrum of the valence band taken with Mg $`K\alpha `$ radiation. Owing to the higher kinetic energies of photoelectrons, this spectrum is considered to be more bulk sensitive than the above PES spectra. Considering the photoionization cross-sections,  the valence-band XPS spectrum should mainly reflect the Fe 3$`d`$ partial density of states (DOS) with significant contributions from Ce 4$`f`$ and Ce 5$`d`$. As seen in the figure, the XPS spectrum shows a line-shape similar to the off-resonance spectra of Ce 4$`d`$-4$`f`$ and Ce 3$`d`$-4$`f`$ resonant PES.

It is expected, that the BIS spectrum should also reflect the bulk electronic structure rather well. There is a peak near $`E_F`$ and a broader feature at $``$6 eV. They originate mainly from the Ce 4$`f`$ states, although there are contributions from the Ce 5$`d`$ and Fe 3$`d`$ states, too. The peak near $`E_F`$ and the structure around $``$6 eV correspond to the 4$`f^1`$ and 4$`f^2`$ final states, respectively.  The broad line-shape of the structure at $``$6 eV is due to the 4$`f^2`$ final-state multiplet.  The strong intensity of the peak near $`E_F`$ again indicates strong hybridization of the Ce 4$`f`$ states with the valence states.

Fig. 3 shows high-resolution UPS data. In this photon energy range, the cross sections of the Fe 3$`d`$, Ce 4$`f`$ and Ce 5$`d`$ states varies rapidly with photon energy.  The relative cross-section of Ce 4$`f`$ to the other orbitals increases when going from $`h\nu `$ =21.2 to 40.8 eV while that of Ce 5$`d`$ rapidly decreases. Therefore the structure at -(2 $``$ 3) eV which appears only in the $`h\nu `$ = 40.8 eV spectrum, originates from the Ce 4$`f`$ states, and corresponds to one of the double peaks in the Ce 4$`f`$ spectrum obtained by the Ce 4$`d`$-4$`f`$ resonant PES. This observation is also consistent with the fact that the 40.8 eV spectrum is surface sensitive, according to the “universal curve” of the mean free path of photoelectrons.  In the near $`E_F`$ region, structures just below $`E_F`$ and at $``$ -0.3 eV are somewhat enhanced in the 40.8 eV spectrum. These structures originate from the Ce 4$`f`$ states and correspond to the tail of the Kondo resonance (possibly with unresolved fine structures due to crystal-field splitting) and the spin-orbit side band, respectively. These structures are also expected to be dominated by surface contributions.

### B Single-Impurity Anderson Model

The SIAM parameters obtained in our calculation are listed in Table I. Here, $`ϵ_f`$ is the position of the bare 4$`f`$ level (4$`f^1`$ $``$ 4$`f^0`$ ionization level) relative to $`E_F`$, $`U_{\mathrm{𝑓𝑓}}`$ is the 4$`f`$-4$`f`$ on-site Coulomb energy, $`U_{\mathrm{𝑓𝑐}}`$ is the 4$`f`$-core-level Coulomb energy and $`V`$ is the Ce 4$`f`$-valence-band hybridization strength, in accordance with the definitions in Ref. . Using those parameters, the $`4f`$ occupation number $`n_f`$ has been calculated and listed in the last column of Table I. Since we have fitted many different types of spectra using the SIAM shown in Fig. 4, many constraints have lead to a rather unique set of SIAM parameters.

In Fig. 4, comparison is made between the SIAM calculations and the experimental spectra. As seen, we obtain good overall agreement with all experimental spectra. The main discrepancy between the SIAM results and experimental spectra is found in the BIS spectrum (Fig. 4c), where the position of the calculated $`f^1`$ peak is about 0.5 eV lower than in the experimental spectrum. Noticeable from this figure is also the large difference between the bulk and surface spectra obtained through the SIAM analysis. For instance, in the Ce $`4f`$ spectrum obtained from Ce $`4d`$-$`4f`$ resonant PES (Fig. 4b), which is basically a double-peak structure, the relative strength of the two peaks is very different between the bulk spectrum and the surface spectrum. Similar differences in directly measured bulk- and surface spectra have previously been reported for Ce-metal, and can be explained as due to larger hybridization in the bulk.

Table I shows that, apparently, the Ce atoms belonging to the surface are as good as completely trivalent, with $`n_f1.0`$, whereas in the bulk, the 4$`f`$ states are strongly hybridized, having the significantly lower occupation of 0.78.

However, there are noticeable amplitudes of the $`f^0`$ and $`f^2`$ configurations also at the surface, indicating that also here, some hybridization between the $`4f`$ and valence states is taking place.

### C Band-Structure Calculation

Experimentally, the partial moments in CeFe2 have been studied using several different experimental methods: polarized neutrons,  Compton scattering and, very recently, x-ray magnetic circular dichroism (XMCD). In all experiments, an anti-parallel coupling of the Ce and Fe moments is found. This coupling is also reproduced in our calculation, as well as in earlier calculations. As is well known, this anti-parallel coupling of the moments is a strong indication that the Ce $`4f`$ states in CeFe2 are delocalized. This can easily be understood from the following argumentation. If the $`4f`$ electrons are localized, the $`4f`$ spin moment would be dictated by the polarization of the $`spd`$-electrons of the Ce atom, that via hybridization effects are known to be anti-parallel to the $`3d`$ moment of the Fe atom. Hence the spin moments of the Ce atom and the Fe atom are always anti-parallel, both in the localized and delocalized case. For localized $`4f`$ electrons, the $`4f`$ spin moment is accompanied by an orbital moment (larger than the spin moment) that (via Hund’s third rule) is anti-parallel to the Ce spin moment. Hence, for localized $`4f`$ electrons the total (spin+orbital) Ce-Fe coupling is ferromagnetic, whereas if the Ce $`4f`$ orbital moment is quenched, due to band formation, the coupling is anti-parallel.

In Fig. 5, the spin-resolved partial DOS for the Ce $`4f`$, Ce $`5d`$ and Fe $`3d`$ states are shown. The first and third panels show the majority spin channel for Ce and Fe, respectively, and the second and fourth panels show the minority spin channel. Comparing the DOS for Ce and Fe, we see that the Ce $`4f`$ and Fe $`3d`$ states have opposite spin polarization. Furthermore, the Ce $`5d`$ band width is seen to be much larger than that of the Ce $`4f`$ and Fe $`3d`$ states, with the magnitude of the Ce $`5d`$ DOS roughly an order of magnitude smaller than that of the Ce $`4f`$ and Fe $`3d`$ states.

Fig. 6 shows the band-structure of spin-polarized CeFe2 along high-symmetry directions in the Brillouin zone. The flat bands clustered just above the Fermi level are predominantly of $`4f`$-character. In the region from $``$2 eV to $``$10 eV, and around -5 eV with respect to the Fermi level, the spin splitting of the bands is clearly visible.

Our calculated total spin magnetic moment amounts to 2.48 $`\mu _\mathrm{B}`$ per formula unit (experimental saturation magnetization: 2.30 $`\mu _\mathrm{B}`$, as stated earlier in this paper), with the main contributions being the Ce $`4f`$ moment -0.54 $`\mu _\mathrm{B}`$, the Ce $`5d`$ moment -0.23 $`\mu _\mathrm{B}`$, and the Fe $`3d`$ moment 1.75 $`\mu _\mathrm{B}`$. The partial occupation numbers summed over spin, within the muffin-tin spheres, are 1.07 for Ce 4$`f`$, 1.28 for Ce 5$`d`$, and 6.18 for Fe 3$`d`$. Note that the partial spin magnetic moments are calculated using the partial occupation numbers inside the muffin-tin spheres, which is somewhat arbitrary. The total spin moment is on the contrary, of course, well defined. An obvious point, which seems to have been overlooked so far, is that not only in band-structure calculations, but also experimentally, the division of space between individual atomic species in a compound is in fact not unique nor even well defined. It is reasonable to believe that different experimental procedures differ in their “volume of sensitivity” around each atom, and thus effectively correspond to different ways of dividing up the total space in the compound between the atoms. This could be one reason why different experimental techniques find quite different values for the partial magnetic moments, and also why, in order to find the total Ce moment from experimental results, assumptions have to be made regarding the ratio of the number of $`5d`$ and $`4f`$ electrons contributing to the magnetization.

An analysis of experimental data along this direction of thought might help resolve controversies regarding the electronic structure of CeFe2.

A calculation of the moments including spin-orbit coupling and orbital polarization, using the FP-LMTO method, has been performed earlier by Trygg et al.. The difference between the presently reported spin moments and the ones reported by Trygg et al., which include spin-orbit coupling, is very small, around 1%. Thus, the effect of including spin-orbit coupling is shown to have only a very minor effect on the magnitude of the spin moment. Calculations by Eriksson et al. using the atomic sphere approximation (ASA) give somewhat different values for the magnetic moments than the present method, in which no such geometrical approximation regarding the form of the potential, wave functions or charge density is made. As demonstrated in Ref. , the spin density in CeFe2 is highly non-spherical, which may well be the reason for the differences in results from full-potential and ASA calculations. To summarize, the arguments presented above justify our present calculational approach, i.e, using a full-potential method, but neglecting spin-orbit coupling.

Concerning the absolute magnitudes of the individual moments, the discrepancies between different experimental approaches can be quite large, for instance, the Ce $`4f`$ spin moment is measured to be -0.37 $`\mu _\mathrm{B}`$ with XMCD, whereas polarized neutrons find the corresponding moment to be only about a fourth as large: -0.10 $`\mu _\mathrm{B}`$. As already touched upon above, one reason for these discrepancies between different experimental techniques may well be that they differ in the way the space in the compound is effectively divided up between the atoms. With this in mind, the magnitudes of our calculated moments must be said to be in satisfactory agreement with experimental findings, although the overall trend appears to be that the calculations overestimate the moment magnitudes.

In Fig. 7, we compare the valence-band Ce 4$`f`$ PES and BIS spectra with the DFT DOS. As mentioned earlier, the Ce 3$`d`$-4$`f`$ resonant PES spectrum and the BIS spectrum are rather bulk sensitive, and thus it is relevant to compare these spectra with bulk DFT calculations. In the PES part of the spectra, comparison is made between the experimentally obtained Ce 4$`f`$ spectra and the Ce 4$`f`$-projected DFT DOS. As for the BIS part, the Ce 4$`f`$, 5$`d`$ and Fe 3$`d`$ partial DOS have been added taking account of the atomic photoionization cross-sections.  Agreement between experiment and theory is satisfactory almost to the same extent as in the SIAM calculation. In the DFT DOS, the structure around 6 eV in the BIS spectrum is of course not reproduced, since this structure corresponds to the 4$`f^2`$ final state, and thus is a purely excited-state property of the system. Furthermore, in the DFT calculation, the energy of the near $`E_F`$ peak in the Ce 4$`f`$ PES spectrum and also in the BIS spectrum is slightly higher than in the experimental spectra. Also, the intensity on the higher binding energy side of the PES spectrum is underestimated in the calculation. However, one should note that the relative intensities depend also on the transition matrix elements, which are not included in the DOS curves.

Fig. 8 shows the valence-band XPS, Ce 3$`d`$-4$`f`$ and 4$`d`$-4$`f`$ off-resonance and UPS (He II) spectra. The spectral weight comes primarily from the Fe 3$`d`$ states, and thus we compare these spectra with the Fe 3$`d`$ partial DFT DOS. Although there are differences in the surface sensitivity and in the contributions from other orbitals, all the experimental spectra have similar band widths and line shapes. In comparison with the DFT DOS, although overall features are well reproduced, the experimental spectra have larger spectral weight near $`E_F`$ than the calculated DOS.

## IV Discussion and Conclusions

Generally, for any material, the surface electronic structure can differ substantially from the bulk one. For valence fluctuating systems like CeFe2, it is very likely that the electronic structure of the Ce atoms close to the surface is not the same as that of the bulk Ce atoms. Therefore, in order to study the bulk electronic structure of valence fluctuating Ce compounds by means of high-energy spectroscopic methods, it is essential to take into account the effects of the sample surface when interpreting the spectra. In the present work, by assuming that the spectra are superpositions of the surface and bulk components, we have shown that all measured spectra of CeFe2 are fairly well reproduced by the SIAM calculations using the same set of parameters. At the surface, Ce is found to be nearly trivalent. The bulk set of parameters places CeFe2 in the strongly intermediate-valent regime, giving a $`4f`$-occupation number $`n_f`$ as small as 0.78 (with considerable amplitude of the $`f^2`$ configuration both in the surface and bulk). This means that the Ce 4$`f`$ states are strongly hybridized with the Fe 3$`d`$ states in the bulk and that the states around $`E_F`$ have a large amount of $`f`$ character. Also, the differences between our SIAM-derived bulk- and surface spectra for CeFe2 are similar to the differences between directly measured bulk- and surface spectra of Ce-metal, a difference which can be explained as due to larger $`4f`$ hybridization in the bulk than at the surface.

Apart from the SIAM analysis, a number of features of the measured spectra force us to draw the same conclusion regarding the nature of the $`4f`$ states in the bulk and at the surface, notably the large difference between the Ce $`4f`$ spectra obtained from the Ce $`3d4f`$ and Ce $`4d4f`$ resonant PES, and the strong intensity of the peak near $`E_F`$ in the BIS spectrum,

As a result of the strong hybridization of the Ce 4$`f`$ states in the bulk, the “$`f^0`$” final state feature in the Ce 4$`f`$ spectrum deduced from the Ce 3$`d`$-4$`f`$ resonant PES has a very weak intensity in contrast to the Ce 4$`d`$-4$`f`$ resonant PES spectrum, where the $`f^0`$-final-state feature leads to the well-known double-peak structure. In such a case, i.e, where there is strong screening of excitations, it is expected that the bulk 4$`f`$ spectrum can be interpreted in terms of a one-electron picture, and thus the DFT DOS should compare well with the $`4f`$ spectrum. This is also seen to be the case. Comparison of the 4$`f`$ spectra with the DOS (Fig. 7) shows that the Ce 4$`f`$ partial DFT DOS describes the valence-band Ce 4$`f`$ spectra well, except for the $`f^2`$ structure in the BIS spectrum of course, since the $`f^2`$-peak is due to incomplete screening. This fact poses the question of how the “$`f^1`$” final state (which is commonly referred to as the “Kondo peak”) of the SIAM picture and “the 4$`f`$ band” in the band picture are related to each other, since according to DFT, this peak is a one-electron feature, and in the SIAM, this is due to a many-body effect.

We also draw the conclusion that due to the strong hybridization in the bulk, the spin-orbit side band seen in the high-resolution UPS spectra (Fig. 3) must have surface origin since it is known that the spectral weight of the spin-orbit side band is strongly reduced when the hybridization is strong.

We now turn to a more detailed comparison of how well the DFT calculations and the SIAM analysis perform for the different spectra. Regarding the position of the near $`E_F`$ peak in the BIS spectrum, the DFT calculation predicts a higher energy than experiment while the SIAM calculation predicts a lower energy than experiment. The intensity on the higher binding energy side of the near $`E_F`$ peak in the PES spectrum is underestimated in the DFT calculation (although strictly speaking, intensities cannot be expected to be reproduced with a DOS, since the transition matrix elements are neglected) while it is overestimated in the SIAM calculation. The position of the near $`E_F`$ peak in the PES spectra is calculated to be too close to $`E_F`$ both in the SIAM and DFT calculations. As for the Fe 3$`d`$ component, the DOS calculated using DFT does not give a good account of the intensity in the experimental spectra near $`E_F`$. This may be due to that the $`4f`$ transition matrix elements are large close to the Fermi level compared to the $`3d`$ transition matrix elements.

All in all, the above discussion amounts to that the Ce $`4f`$ electrons in the bulk hybridize strongly with the Fe $`3d`$ electrons. This conclusion agrees perfectly with the experimentally observed anti-parallel coupling of the Ce and Fe moments, which is also reproduced in the DFT calculation. In the SIAM, the $`4f`$ electron is assumed to be localized, which indirectly implies that a parallel coupling of the Ce and Fe moments is expected.

Finally, we wish to mention some sources of error in the present work. Our measurements were done on scraped surfaces, which might make the surface rough, and therefore ill-defined. Furthermore, the precise values of photoelectron mean free paths are difficult to estimate, which naturally also has the effect of making the border between “bulk” and “surface” somewhat ill-defined. Our SIAM is not the most elaborate one, for instance we assume a degeneracy of 14 of the $`4f`$ level, thereby neglecting spin-orbit coupling and anisotropic hybridization effects, which leads to a crystal-field splitting. Furthermore, as in all DFT calculations, the functional used treats electron correlation only to a limited extent, i.e., it is not meaningful to expect perfect agreement between the DFT results and experiment.

## Acknowledgment

The authors thank the staff of SOR-RING, Y. Azuma, T. Miyahara and the staff of Photon Factory for technical support. The authors also thank A. Sekiyama, J. Okamoto, T. Tsujioka and T. Saitoh for help in the experiment. A.D. acknowledges financial support from the Swedish Foundation for International Cooperation in Research and Higher Education (STINT). Part of this work has been done under the approval of the Photon Factory Program Advisory Committee (Proposal No. 94G361). A.D. and O.E are grateful to J. M. Wills for supplying the full potential code used in this study.


format: text

text:
# The Effect of the Mass Accretion Rate on the Burst Oscillations in 4U 1728–34

## 1 Introduction

Low mass X-ray binaries (LMXBs), systems in which a compact object accretes matter via Roche lobe overflow from a low mass ($`<`$ 1 $`M_{}`$) main sequence companion, provide a fascinating window into the physics and dynamics of accretion disks and compact objects. When the compact object is a neutron star, these systems also permit detailed studies of the burning of the accreted material on the neutron star surface. Thermonuclear flashes (called type-I X-ray bursts) occur in many LMXBs and are caused by explosive unstable burning of the accreted matter. The properties of these type-I X-ray bursts are affected by the rate at which the matter is deposited onto the neutron star (see, e. g., Lewin, van Paradijs & Taam, 1993, for an extended review of X-ray bursts). Thus, in order to understand these bursts, it is necessary to investigate the effect that variations in the mass accretion rate ($`\dot{M}`$) onto the neutron star have on them. The first attempts at characterizing $`\dot{M}`$ in LMXBs centered around the persistent flux of the sources, since intuitively the larger the amount of matter accreted, the greater the amount of X-ray flux emitted. However, it has become clear that the persistent flux in general may not always be a good measure of $`\dot{M}`$ (e. g., van der Klis, 1995) and hence cannot be directly correlated with the burst properties.

X-ray color-color diagrams (CDs) were introduced to study subtle spectral changes in LMXBs (see, e. g., van der Klis, 1995, and references therein). Hasinger & van der Klis (1989) studied the correlated rapid X-ray variability and spectral behavior of LMXBs and showed that they naturally fall into two distinct categories. These two groups are called the atoll and Z sources, after the tracks they trace out in the CD. The motion of the sources along these tracks is thought to be caused by changes in the mass accretion rate (Hasinger & van der Klis, 1989) and it has become clear that the position of the sources on these tracks is a better indicator of the $`\dot{M}`$ of the system than the persistent flux (see van der Klis 1995 for LMXBs in general but Méndez 2000 for 4U 1728–34 in particular).

With the launch of the Rossi X-Ray Timing Explorer (RXTE) satellite, previously unknown high frequency timing phenomena were discovered and now offer us new views of LMXBs (see, e. g., van der Klis, 2000, for an extended review of the recently discovered timing phenomena in LMXBs). In particular, nearly coherent brightness oscillations (NCOs) with frequencies between 330 and 589 Hz were discovered in Type I bursts from several sources (Strohmayer et al., 1996; Strohmayer et al., 1997a; Smith, Morgan, & Bradt, 1997; Zhang et al., 1997, 1998; Markwardt, Strohmayer, & Swank, 1999; Wijnands, Strohmayer & Franco, 2000). These NCOs have extremely high coherence values ($`\nu /\mathrm{\Delta }\nu >4000`$; Strohmayer & Markwardt, 1999; Muno et al., 2000), high amplitudes (up to $``$50% rms; Strohmayer, Zhang & Swank, 1997b; Strohmayer et al., 1998b), and they exhibit long term frequency stability of up to years (Strohmayer, Zhang & Swank, 1998a). These properties strongly argue for the rotation of the neutron star as the origin of the NCOs. In this interpretation, during a thermonuclear burst inhomogeneities in the fuel or in the fuel burning may arise causing “hot spots” on the surface of the neutron star. These hot spots rotate with the spin frequency of the star and periodically modulate the X-ray flux.

These NCOs can be present in the rise of the bursts, very near the burst onset, and/or in the decay phase (or tail) of the bursts. When the NCOs are observed in the tail, their frequency can increase by 1 – 3 Hz, reaching an asymptotic limit as the bursts progress (see e. g., Strohmayer et al., 1998b). It has been proposed that this frequency increase is caused by thermal expansion of the thermonuclear burning layer on the stellar surface. Heating from the burning expands the shell by several tens of meters at the onset of the burst. Due to angular momentum conservation this layer will slow down causing the observed frequency to decrease. During the decay of the burst, the burning layer settles back on the neutron star surface and spins up causing the frequencies to increase again. In this picture, the asymptotic limit represents the true neutron star spin frequency (Strohmayer, Zhang & Swank, 1997b). The atoll source 4U 1636–53 has shown evidence for a subharmonic to the NCO frequency, indicating the possible existence of two, nearly antipodal hot spots (perhaps at the magnetic poles; Miller, 1999). This detection of the subharmonic raises issues about how the burning front propagates on the surface of the neutron star and how the fuel can be concentrated in two antipodal regions.

The burst oscillations in the atoll source KS 1731–260 have been recently studied by Muno et al. (2000). They reported the presence of a correlation between the inferred $`\dot{M}`$ of the system and the presence of NCOs in the bursts. They find that the bursts from KS 1731–260 can be separated into two classes: the so called “fast” bursts, which occur at relatively high $`\dot{M}`$, exhibit NCOs, and show episodes of photospheric radius expansion (PRE, the expansion of the neutron star photosphere occurring when the burst luminosity is higher than the Eddington limit), and the “slow” bursts, which occur at lower $`\dot{M}`$, do not exhibit NCOs, and do not show PRE episodes. Studies like this require a large sample of bursts observed over a considerable range of accretion rates. Such data bases have so far only been obtained for the persistent sources because the outbursts of the X-ray transients are rare and when they do occurr only a limited number of X-ray bursts are seen. In this paper, we will concentrate on the persistent source 4U 1728–34.

The LMXB 4U 1728–34 is a prolific source of type I X-ray bursts (Basinska et al., 1984) and one of the best studied LMXBs. It has been classified as an atoll source (Hasinger & van der Klis, 1989; Méndez & van der Klis, 1999) and it exhibits most of the behaviors seen in those systems. It was the first neutron star LMXB to show twin kilohertz quasi-periodic oscillations (kHz QPOs) in the persistent emission and the first to show NCOs (Strohmayer et al., 1996). The kHz QPOs in 4U 1728–34 have been studied by several groups (Strohmayer et al., 1996; Ford & van der Klis, 1998; Méndez & van der Klis, 1999; Di Salvo et al., 2000b) and the NCOs by Strohmayer et al. (1996), Strohmayer, Zhang & Swank (1997b, 1998a), and Strohmayer & Markwardt (1999). In this paper we present a study of the behavior of NCOs in a set of 38 type-I X-ray bursts from 4U 1728–34, observed over a wide range of mass accretion rates. We found that 4U 1728–34 is similar to KS 1731-260 in that the NCOs are only present in bursts which occur at high $`\dot{M}`$ (Muno et al., 2000). However, contrary to the case of KS 1731–260, we find an anti-correlation between the presence of PRE episodes in the bursts and the inferred mass accretion rate of the system. We also find a more complicated behaviour with $`\dot{M}`$ for the NCOs in the rise and in the decay phases of bursts.

This paper is organized as follows: in §2 we describe in detail the observations being used for this study. In §3 we explain the methods used to characterize the persistent emission and the burst properties, and the techniques used to study the oscillations. In §4 we discuss the spectral states (i. e., the CD) of the source. In §5 we discuss the behavior of the NCOs and study them in correlation with other burst characteristics and the spectral state of the source. Finally in §6 we discuss the implications of these correlations for our current understanding of 4U 1728–34 in particular and for thermonuclear burning on the surface of neutron stars in general.

## 2 The observations

The source 4U 1728–34 was observed with the Proportional Counter Array (PCA) on board RXTE on many occasions between 1996 February 15 and March 1 (hereafter referred to as observation set AO1), 1997 September 19 and October 1 (AO2), and 1998 September 30 and 1999 January 18 (AO3). The actual observation dates are listed in Tables 2, 2 and 4, along with the main properties of all 65 observations in our sample. In total $``$545 ks of on-source data were obtained. The AO1 observations were obtained with epoch 1 gain settings for the PCA instrument, while AO2 and AO3 observations correspond to epoch 3 gain settings. Hereafter, when we refer to an observation we mean those data that can be identified by an unique observation identification number or ObsID (see Tables 2, 2, and 4). All 5 PCA detectors were on during all observations except for 3.6 ks during the observation with ObsID 20083-01-01-020 and 3.6 ks during 20083-01-03-020. No bursts were detected during those two intervals, so we have excluded those data from our study. The analysis of portions of our data have already been published in the literature. Strohmayer et al. (1996) reported the discovery of NCOs using burst 5 in our sample and Strohmayer, Zhang & Swank (1997b) presented a temporal and spectral analysis of bursts 3, 4 and 5, and reported the existence of NCOs in bursts 4 and 5. Méndez & van der Klis (1999) published a CD of AO1 and a portion of AO2 data (see also Di Salvo et al., 2000b).

A variety of data modes was available for each observation set in addition to the Standard 1 and Standard 2 modes. Because the same data modes were

not always available across observation sets, we describe here the extant modes and how they were used to extract the quantities of interest to this study. The AO1 data were obtained in one burst trigger mode (TLA\_1s\_10\_249\_1s\_5000\_F), two burst catcher modes (CB\_125us\_1M\_0\_249\_H with 122 $`\mu `$s time resolution in one energy channel; CB\_8ms\_64M\_0\_249\_H with 8 ms resolution in 64 channels) and one single bit mode (SB\_125us\_0\_249\_1s with 122 $`\mu `$s resolution in one channel). The burst catcher modes provide data typically for 1 s before the burst onset and for as long as the count rate remains above the trigger level. In our case the length of data ranges from 13 to 16 s for these modes. For AO2 we had the same burst trigger and catcher modes available as during AO1, however, instead of the single bit mode we now had one event mode (E\_16us\_64M\_0\_1s with 16 $`\mu `$s resolution and 64 channels) available. The data for AO3 were obtained in three single bit modes (SB\_125us\_0\_13\_1s, SB\_125us\_14\_17\_1s and SB\_125us\_18\_23\_1s all with 122 $`\mu `$s resolution and one channel) and one event mode (E\_125us\_64M\_24\_1s with 122 $`\mu `$s time resolution and 64 energy channels). When combined, all four data modes for AO3 data cover the complete RXTE energy range. To generate the CD we extracted light curves from the Standard 2 data in 4 energy ranges to calculate the colors. We also used the Standard 2 data to generate spectra of the persistent emission and extract average fluxes. The power density spectra of the persistent emission were calculated from the single bit mode data for AO1, the event mode data for AO2, and the data of the combined modes for AO3. These same data modes were also used to search for NCOs in the bursts. Finally, the time series of the burst spectra were obtained from the event mode for AO2, the combined data modes for AO3, and the 8 ms resolution burst catcher mode for AO1 (hence there is at most 16 s of time resolved spectra for each burst in AO1).

## 3 The Analysis

To characterize the persistent emission, we calculated average fluxes for each observation. The fluxes were obtained by fitting (using Xspec; Arnaud 1996) a cut-off power law spectrum to the background-subtracted Standard 2 data. A Gaussian line around 6 keV was added in order to obtain acceptable fits. Since we are only interested in the flux of the persistent emission, we did not investigate the significance of this feature further (but see, e. g., Di Salvo et al., 2000a). The resulting averaged fluxes in the 3–18 keV range are listed in Tables 2, 2, and 4. In order to

verify the reliability of our fluxes we fit other spectral models (3-component models: a comptonized spectrum plus either a multi-color disk blackbody or a plain black body plus a Gaussian line around 6 keV; see, e. g., Di Salvo et al. 2000a). We find that these alternative models give the same average flux of the persistent emission within the quoted errors.

To generate the CD, for AO1 we define the soft color as the logarithm of the ratio of the count rate in the energy range 2.9–4.5 keV to the count rate in the range 4.5–6.3 keV. Similarly we define the hard color as the logarithm of the ratio of the count rate in the range 6.3–9.0 keV to that in the range 9.0–18.2 keV. For AO2 and AO3 the energy ranges are slightly different: 3.0–4.4 keV to 4.4–6.2 keV for the soft color and 6.2–9.1 keV to 9.1–18.2 keV for the hard color. The colors were calculated from background-subtracted light curves from the Standard 2 data and hence have 16 s time resolution. However, to reduce the errors on the colors we averaged them in 256 s intervals.

We searched for bursts by generating light curves of the Standard 1 data with 0.125 s resolution. We found 38 bursts in our data set. We have numbered the bursts according to time and have listed their main properties in Table 4. The start time of the burst is determined from the first time bin in the background-subtracted light curve with count rate above 1000 counts s-1. We characterize the bursts by obtaining time resolved series of spectra. We extracted dead-time corrected spectra of 0.125 s duration for a time interval of 20 s during the bursts.999Except for AO1 observations where the only data mode available with high enough time and energy resolution was a burst catcher mode with only up to 16 s per burst. We

single trial significance of a little better than $`1\times 10^3`$. Since pulsations are known to exist in this source we do not need to search a wide frequency space and we therefore have few effective trials. Using this criteria we found that 16 out of our sample of 38 bursts exhibit NCOs. Typical examples of bursts which show NCOs are shown in Figure 2. We established a measure of the strength of the oscillations by integrating the RMS amplitude of the oscillations over their duration in the burst, ignoring the wings of the power features which fall below our detection threshold.

## 4 The Color-Color Diagram of 4U 1728–34

The CD of 4U 1728–34 in Figure 3a shows a track typical of atoll sources (Hasinger & van der Klis 1989, see also Méndez 2000). It is important to note that the shift in the CD of the track for AO1 with respect to the tracks of AO2 and AO3 is not an intrinsic behavior of the source but due to the difference in gain settings. Our observations of 4U 1728–34 sampled both the island state and the banana branch, though not during a single observation set. Both AO1 and AO3 exhibit the complete island state track, from the extreme island state at the upper right of the CD, believed to be the state with the lowest mass accretion rate, to just before the island state–banana branch vertex (hereafter simply referred to as ’the vertex’). Only during our AO2 observations was 4U 1728–34 found in the banana branch. This observation set covers the entire banana branch, from the lower banana branch on the left corner to the upper banana branch on the bottom rightmost portion of the track, where the mass accretion rate is thought to be the highest. Our CD looks very similar in shape to those published previously for 4U 1728–34 (Méndez & van der Klis, 1999; Piraino et al., 2000; Di Salvo et al., 2000b). Our island state appears somewhat longer, but this difference is due to the different choice of colors and disappears when similar colors are used.

To determine the position of the source in the CD prior to each burst, we determine the color of the source from an interval of 320 s before the burst, ending 16 s prior to the onset of the burst. The results can be seen in Figure 3a where the numbers indicate the spectral state of the source just prior to the occurrence of each burst. From the figure it is clear that we detected bursts throughout almost the entire track, from the extreme island state to the upper banana branch although not in the uppermost part of the upper banana branch. However, it is not yet clear whether this is due to a real lack of bursts at the highest $`\dot{M}`$ or due to the limited amount of data so far obtained when 4U 1728–34 was at such high accretion rates.

To verify the spectral state of the source prior to each burst, average power density spectra (PDS) were obtained by integrating over the entire length of the continuous data train prior to each burst. Figure 4a shows the PDS of the source prior to burst 22, which is typical

of the source when it is in the extreme island state. This PDS clearly shows band-limited noise at frequencies above 100 Hz, with no high frequency QPOs, a bump around 10 Hz, a flat spectrum at low frequencies, and a break around 1 Hz, typical for the island state (see, e. g., van der Klis, 1995). At the vertex, e. g., prior to burst 21, the source PDS (Fig. 4b) shows a kHz QPO and a break around 20 Hz typical of the lower banana branch in transition to the island state. Figure 4c shows the PDS for the source prior to burst 15, when 4U 1728–34 is in the upper banana branch. The PDS shows a typical featureless power law noise component. As mentioned above, the difference in gain settings between epoch 1 and 3 causes a shift in the island state tracks between AO1 and AO3. To see where the AO1 track might join the AO2 track if there were no gain differences, we show in Figure 4d the PDS of the source prior to burst 3. This PDS is very similar to that of burst 21 (Fig. 4b), indicating a transition between the island state and the lower banana branch.

## 5 The Behavior of Burst Oscillations

### 5.1 The NCOs and the Inferred Mass Accretion Rate

As shown in Figure 6, we find no correlation between the average flux of the persistent emission prior to the bursts and the existence or the strength of the NCOs. However, if the presence and possibly the strength of the NCOs is correlated with the mass accretion rate, as was the case in KS 1731–260 (Muno et al. 2000), then the lack of correlation is not surprising, because, as already mentioned in § 1, the flux is not a robust tracer of the mass accretion rate. Instead, a better measure of $`\dot{M}`$ is the position of the source on the atoll track in the CD.

In Figure 3b, we plot all bursts in the CD; with diamonds the bursts that do not show NCOs, and with circles those that do. For bursts with oscillations the size of the circles is proportional to the strength of their NCOs as defined in § 3. From this figure it is clear that bursts with NCOs occur only on the bottom part of the diagram where the accretion rate is thought to be highest, i. e., near the vertex and on the banana branch. Close inspection of this figure reveals that the strength of the NCOs generally increases with decreasing hard color, i. e., the NCOs become stronger as the inferred mass accretion rate of the system increases. This relation is especially evident for bursts within a single observation set, and hence separated in time by only a few days to a few weeks, but it is not so clear for bursts from different observation sets, separated in time by up to years (i. e., between the AO1 and the AO2 data sets).

### 5.2 The NCOs and the Burst Flux

In Figure 6 we show the peak burst fluxes for the bursts in our sample. Bursts without NCOs are marked with asterisks, while bursts with NCOs are shown with circles whose radius is proportional to the strength of the oscillations. Burst 3 which shows an unusual count rate profile is marked with a square. The first trend evident from Figure 6 is that the bursts in our AO2 observations (except for bursts 20 and 21) have considerably lower peak fluxes than those in our AO1 and AO3 observations. We also observe that the bursts without NCOs have on average higher peak fluxes. Furthermore, among the bursts with oscillations in a single observation set (AO1 or AO2), the dimmest bursts also have the strongest oscillations. This correlation does not hold when we consider bursts from AO1 and AO2 together.

### 5.3 The NCOs and Bursts with Photospheric Radius Expansion

In their study of KS 1731-260, Muno et al. (2000) found that all the bursts with PRE episodes occur at high mass accretion rates and they all show NCOs. They go further to distinguish between oscillations occurring during the rise of the burst and those occurring during the cooling tail of the burst. They note that oscillations in the rise of bursts are not correlated with PRE episodes since burst 7 in their sample exhibits such NCOs but not a PRE episode. However, their data suggests that all bursts with NCOs in the decay portion of the burst do exhibit PRE episodes.

A careful look at Figure 3c reveals a different picture for 4U 1728–34. Bursts with PRE episodes (shown with squares) occur throughout the entire island state and at the vertex, but not in the banana branch itself. This is precisely the opposite of the behavior in KS 1731-260 where all bursts with PRE episodes lie in the banana branch, at the highest $`\dot{M}`$. We also note that all 7 bursts from our sample that lie in the banana branch (but not bursts 20 and 21 which lie near the vertex) have very low peak measured fluxes (see Fig. 6) while all the PRE bursts in the KS 1731–260 study are “fast” bursts which show high peak fluxes.

The behavior of the bursts in 4U 1728–34 is richer still. Muno et al. (2000) hinted at a possible difference between the NCOs observed during the rise of bursts and those observed in the cooling tails of the bursts. If we make this distinction we find in our sample 2 bursts with NCOs only in the rise, 7 bursts with NCOs in both the rise and the tail, and 7 bursts with NCOs only in the tail. Representative bursts of each of these categories are shown in Figure 2. In Figure 3d we can see that these distinctions in the NCOs are also borne out on the CD. Bursts with NCOs only in the decay tail of the burst are seen (marked with squares) at intermediate $`\dot{M}`$, in the island state and the vertex. Those bursts with NCOs in both the rise and the tail of the bursts are found (marked with triangles) only at high $`\dot{M}`$, in the banana branch. The two bursts with NCOs only in the rise of the burst are shown with circles and they appear in the upper banana branch, at the highest inferred $`\dot{M}`$. However, with only two such bursts it is difficult to determine whether or not their location only in the UB branch is significant. The bursts without NCOs are shown with diamonds.

## 6 Discussion

We have searched three years of archival RXTE data of the neutron star LMXB and atoll source 4U 1728–34 for type-I X-ray bursts. We have found 38 bursts, 16 of which show nearly coherent brightness oscillations (NCOs). The presence or absence of the NCOs in the bursts is not well correlated with the persistent flux of the source, but it is very well correlated with the position of the source on the atoll track in the CD. The NCOs are only present in the bursts that occur when the source is near the island state-banana branch vertex or on the banana branch. It is thought (Hasinger & van der Klis, 1989) that the position of the source in the atoll track is determined by the mass accretion rate. If true, then the presence of the NCOs is very well correlated with the inferred mass accretion rate and the lack of any correlation with the persistent X-ray flux is further proof that the persistent flux is not a good measure of the accretion rate of the system.

A similar study for KS 1731–260 by Muno et al. (2000), although with less bursts (9 bursts, 5 with NCOs), also found that the NCOs are only present in bursts which occur at relatively high mass accretion rates. However, Muno et al. (2000) also report that the bursts of KS 1731–260 which have episodes of photospheric radius expansion (PRE) occur at these high accretion rates and that they all exhibit NCOs. Such a correlation can be excluded for 4U 1728–34. When we considered bursts that exhibit (PRE) episodes, we discovered a strong anti-correlation with inferred mass accretion rate: bursts with PRE episodes occur only when the source is at its lowest and intermediate mass accretion rates and not at its highest accretion rates. Also, the presence of NCOs during the bursts is not correlated with the presence of PRE episodes: both bursts with and without NCOs can exhibit periods of PRE. Yet another difference we found is that the bursts which occured during relatively high mass accretion rates were the dimmest in our sample, while the bursts in KS 1731–260 were bright when the accretion rate in that system was high (Muno et al., 2000). It is unclear why the positions of bursts in the CD are very similar with respect to the existence of NCOs in both 4U 1728–34 and KS 1731–260, but other burst properties differ so significantly with CD position (and thus with the inferred mass accretion rate).

One clue to these differences may be found in the theory of thermonuclear bursts. Atoll sources accrete mass at rates in the range 0.01 $`\dot{M}_{\mathrm{Edd}}`$$`<`$ $`\dot{M}`$$`<`$ 0.1 $`\dot{M}_{\mathrm{Edd}}`$ (Hanawa & Fujimoto, 1982). In this regime, the model of Fujimoto, Hanawa & Miyaji (1981) for thermonuclear bursts (see also Bildsten, 1997) predicts two main behaviors dominated by the accretion rate. At relatively low mass accretion rates, the hydrogen on the neutron star surface burns stably forming a layer of pure helium. As more matter accumulates this layer is compressed and heated until it ignites unstably, producing a pure helium burst. As the mass accretion rate increases, hydrogen is accreted faster than it can be burned and when the helium flash ignities it does so in a mixed hydrogen/helium environment. For atoll sources, bursts occurring in the island state at low inferred mass accretion rates would then be the result of pure helium flashes, with the fraction of hydrogen fuel increasing as the source progresses toward the banana branch.

Because of the relative timescales for energy release of the reactions involved (the CNO cycle for hydrogen burning is limited by the rate of $`\beta `$-decay and thus much slower than the triple alpha process for helium burning), bursts with high fractions of helium are expected to release their energy faster and reach higher peak fluxes. This behavior seems borne out in the bursts from 4U 1728–34. From Table 4 and Figure 3a we can see that in our AO1 observations, the bursts in the extreme island state (bursts 9–12) have higher peak fluxes than the rest of the bursts in that observation set. Similarly within the AO3 observations the bursts at the lowest $`\dot{M}`$ (22–27 and 38) also have the highest peak fluxes among their group. It is natural then that episodes of photospheric radius expansion are anti-correlated with mass accretion rate; the lower the $`\dot{M}`$, the higher the helium fraction, the faster the energy release and the brighter the bursts, making it more likely that they will reach the Eddington limit and expand the photosphere.

In the case of KS 1731–260, Muno et al. (2000) argue that the inverse behavior of this source with respect to the theory’s predictions may be explained by considering not the overall mass accretion rate but the rate of mass accretion local to the fuel surface. If the area over which the matter is being accreted increases with mass accretion rate, the local accretion rate may in fact decrease, causing the inverse behavior (Bildsten, 2000). The question then becomes, why does the area over which the accreted matter is distributed increase with increased $`\dot{M}`$ in KS 1731–260 but not in 4U 1728–34? The answer may lie in the intrinsic differences of the neutron stars in each of these systems. For example, a significant difference in the magnetic field (i. e., its strength and orientation) of the neutron star in 4U 1728–34 compared to that in KS 1731–260 may inhibit the motion of matter from the equator to higher latitudes. Alternatively, if the frequency of the NCOs is indeed the true spin frequency of the neutron star, then KS 1731–260 with an NCO asymptotic frequency of 524 Hz is spinning nearly 45% faster than 4U 1728–34 (with an asymptotic NCO frequency of 363 Hz). This may reduce the effective gravity at the equator of the neutron star in KS 1731–260 and allow more significant spread of the accreted matter.

If KS 1731–260 and 4U 1728–34 are so different with respect to their burst behavior, why then are burst oscillations not observed at low inferred mass accretion rates in either of these sources? Evidently the mechanism producing the oscillations is strongly influenced by the mass accretion rate of the system. Since the NCOs are likely closely tied to the properties of the nuclear burning it is very possible that we are seeing the direct influence of mass accretion rate on the burning behavior. Muno et al. (2000) point to recent work by Cumming & Bildsten (2000) which indicates that modulations of the burst emission are more likely observable during helium-rich bursts, which they argue occur in the banana branch for KS 1731–260 where the NCOs are observed in this source. This cannot be the case for 4U 1728–34 as we have argued that the bursts with NCOs in 4U 1728–34 have in fact high fractions of hydrogen mixed in which causes them to be dim. Further work on possible “hot spot” or other alternative mechanisms will greatly illuminate this issue.

Although it is likely that a low mass accretion rate directly inhibits the oscillation production mechanism, we explore the possibility that these oscillations are in fact present during the burst but that environmental changes (due to the differences in $`\dot{M}`$) cause the oscillations to become undetectable at low $`\dot{M}`$. Recently, the X-ray eclipsing source X 1658–298 was found to also exhibit NCOs during five bursts (Wijnands, Strohmayer & Franco, 2000). Several more bursts were found which did not show these oscillations. However, there is strong evidence that the absence of NCOs in those bursts was caused by environmental changes in the system. In X 1658–298, the oscillations were undetectable or only marginally detectable at times when the source exhibited clear episodes of dipping behavior, making it likely that the mechanism behind the dips causes the oscillations to become undetectable (Wijnands, Strohmayer & Franco, 2000). The dips can be understood as an obscuration of the inner system by matter periodically coming into the line of sight. This same matter then also attenuates the oscillations making them undetectable. Although no X-ray dips have been observed for both 4U 1728–34 and KS 1731–260, it is possible that the absence of NCOs at low mass accretion rates in these sources might also be due to changes in the environment. It is known that when $`\dot{M}`$ decreases in all atoll sources, the X-ray spectrum becomes significantly harder and follows a power law spectrum (see, e. g., Barret et al., 2000, for a recent study of the hard X-ray spectrum of several atoll sources). It has been postulated that this power law spectrum is due to a comptonizing medium (e. g., a corona) around the inner system. If such a medium is indeed present at low $`\dot{M}`$ but not at high $`\dot{M}`$, it might be responsible for smearing out the oscillations making them unobservable at low $`\dot{M}`$. Similar mechanisms have also been introduced to explain the absence of the expected coherent pulsations in the persistent X-ray emission of neutron star LMXBs.

However, this interpretation is not without problems. It has been shown that the only LMXB for which coherent millisecond oscillations have so far been found during the persistent emission (SAX J1808.4–3658; Wijnands & van der Klis 1998) was in the island state (as judged both from the aperiodic rapid variability and the power law spectrum; Wijnands & van der Klis 1999; Gilfanov et al. 1998) during the observations which led to its discovery as a millisecond X-ray pulsar. So, if a comptonizing medium is present in 4U 1728–34 and KS 1731–260 which blurs the NCOs during the type-I X-ray bursts, why are the pulsations in SAX J1808.4–3658 not attenuated? This suggests that the mass accretion rate has a more direct influence on the production mechanism of the NCOs.

Such a direct influence of $`\dot{M}`$ on the NCOs is also suggested by the fact that for seven bursts in 4U 1728–34 the oscillations were visible during both the rise and the decay phases of the burst, while seven bursts show oscillations only in the decay phase and two bursts show oscillations during the rise but not during the tail of the burst. We find that in 4U 1728–34 bursts with NCOs only during the tail occur at intermediate mass accretion rates, while bursts with NCOs in both the rise and the decay occur at high mass accretion rates. Our data further suggest that bursts with NCOs only in the rise occur at the highest mass accretion rates. Clearly, the exact moment when the NCOs will be present during the bursts is determined by the mass accretion rate, which strongly suggests that the NCOs are indeed heavily influenced by the mass accretion rate directly and not by an $`\dot{M}`$-induced environmental change.

The behavior of the bursts, in particular the burst oscillations, in 4U 1728–34 is very rich and many questions remain for which we have no adequate answers. Although the production mechanism behind the NCOs is not well understood, their coherence and their frequency stability over several years strongly suggest that indeed the oscillations are due to the spin frequency of the neutron star. But it is still difficult to imagine how the NCOs in the tail of bursts are formed. The entire surface of the neutron star is expected to be involved in the thermonuclear flash after the first few seconds. In such an environment a hot spot would probably not survive. A possible solution to this quandary is to invoke a different mechanism to produce the NCOs in the decay phase of bursts. However, burst 17 in our sample (see Fig. 2) shows what appears to be a continuous evolution of the NCO frequency. This strongly suggest that only one mechanism is behind the NCOs.

When we observe the oscillations in the tail of the burst, we usually also see an increase in the frequency by a few Hertz as the burst progresses (see also Strohmayer et al., 1996; Strohmayer, Zhang & Swank, 1998a; Strohmayer & Markwardt, 1999). The observed evolution of the oscillation frequency can be convincingly attributed to the expansion of the neutron star atmosphere and conservation of angular momentum (see § 1). However, for at least one of our bursts (burst 17, see Fig. 2) we observe an increase in frequency of about 3.5 Hz. This would mean that the burning layer must have expanded by almost 50 meters. It is not clear if such a large expansion can be accounted for by the expanding layer interpretation, although the expansion can be at least as large as 40 meters (Cumming & Bildsten, 2000). Such a large increase in frequency is not unique to 4U 1728–34. An increase of 5 Hz was recently found in the new NCOs source X 1658–298 (Wijnands, Strohmayer & Franco, 2000), which is even larger than what we have found for 4U 1728–34. However, because of the higher frequency of the NCOs in X 1658–298 ($``$567 Hz) compared to that in 4U 1728–34, the burning layer has to have expanded more in 4U 1728–34 than in X 1658–298 (50 meters versus 40–45 meters).

When we investigated the strength of the oscillations in detail, we found that it is well correlated with inferred mass accretion rate: the oscillations increase in strength when the accretion rate of the system increases. This correlation is more evident within bursts observed only a few weeks apart, and becomes less clear for bursts separated in time by several months to years. It is unclear why the correlation is strong on the short term (weeks) but fails to hold on longer timescales (months to years). Whatever the modulation mechanism, it will likely depend on the state of the fuel layers on the neutron star, which may change considerably with long term trends in the mass accretion rate. If we think of the state of the neutron star surface as determining the moment for the onset of the modulation mechanism, then one can imagine that the overall changes in $`\dot{M}`$ over months and years are resurfacing the neutron star and hence the “critical” moment will be reached at slightly different accretion rates on subsequent trips down the atoll track. In other words, if the time averaged accretion rate has been very low for a long interval of time, the structure of the neutron star surface might be different than when the time averaged accretion rate was considerably larger.

The influences of the mass accretion rate on the burst oscillation properties strongly suggest that they are indeed manifestations of physical properties of the atoll systems. We would expect that this behavior of the NCOs will be also present in other sources. However, no detailed studies have been performed for the other sources similar to our study of 4U 1728–34 and that of KS 1731–260 by Muno et al. (2000). It is important to determine whether the same correlation will hold for all the sources which have NCOs. However, we can definitely exclude the possibility that all X-ray burst sources will show NCOs when they are at high mass accretion rates. For at least two sources (Serpens X-1 and 4U 1705–44), X-ray bursts have been observed when they were on the banana branch, thus at relatively high mass accretion rates. No NCOs were discovered during those bursts (R. Wijnands 2000, private communication; Ford, van der Klis & Kaaret, 1998). A more intrinic reason (e. g., binary inclination, the strength and/or orientation of the magnetic field) must be found to explain the many burst sources that do not show NCOs. On the other hand, it is still possible that a considerable fraction of those sources have indeed only been observed in their island states and might still show NCOs when they can be observed to burst in their banana branches.

To conclude it is clear that from the study of 4U 1728–34 and KS 1731–260 we can already glean some common behaviors and many differences. More atoll systems must be searched for NCOs in their bursts to begin differentiating between those trends that are universal accross atoll sources and those that are determined by intrinsic differences of the systems.

Note added in manuscript. When we were about to submit our paper we became aware of the paper in preparation by van Straaten et al. (2000) presenting a similar analysis of 4U 1728–34. Although their study is limited to the AO1 and AO2 observations, they confirm our main conclusions: the NCOs are only present during bursts at the highest inferred mass accretion rates and the strength of the NCOs increases with accretion rate. They distinguish between two types of PRE episodes, those marked by a small increase in radius after the main PRE episode, which they termed unusual, and the more standard continuous radius decrease after the PRE episode. They find that standard PRE bursts only show NCOs after the peak of the burst and unusual PRE episodes occur only in bursts at the lowest $`\dot{M}`$. We note that these two correlations also hold when including our AO3 data set.

I thank Craig Markwardt for extremely helpful and well documented software, Angela Olinto for unwavering support, Rudy Wijnands for several critical readings of the manuscript and Tod Strohmayer for unselfishly donating his data, time and knowledge and for teaching me all about RXTE and LMXBs. I thank Steve van Straaten and collaborators for their willingness to exchange manuscripts before publication. This work made use of public RXTE data obtained through the GSFC HEASARC archive, and was supported in part by NASA through grant NRA-99-01-ADP-195, by NSF through grant AST-0071235 by DOE through grant DE-FG0291.


text:
# Proof of the thermodynamical stability of the E′ center in SiO2

## Abstract

The E center is a paradigmatic radiation-induced defect in SiO2 whose peculiar EPR and hyperfine activity has been known since over 40 years. This center has been traditionally identified with a distorted, positively-charged oxygen vacancy V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$. However, no direct proof of the stability of this defect has ever been provided, so that its identification is still strongly incomplete. Here we prove directly that distorted V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$ is metastable and that it satisfies the key requirements for its identification as E, such as thermal and optical response, and activation-deactivation mechanisms.

Understanding defects in solids is a key factor in improving device performance and materials quality. Defect identification combines experimental observation and theoretical predictions, and a major ingredient in this process is the theoretical justification of the thermodynamical stability of the defect. If this information is missing, the identification is arguably incomplete or uncertain. Surprisingly, this is the case for the E center , a radiation-induced defect first observed experimentally as early as 45 years ago  in SiO2, a material of prime current importance in microelectronics and fiber optics . E is traditionally identified with a positively-charged distorted oxygen vacancy V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$ , with important support from calculations of hyperfine couplings  and optical activity ; its thermodynamical stability, however, was never theoretically proven, and the mechanisms involved in its activation and deactivation are still debated. Here, using the ab initio theory of defect formation in solids , we show that the conditions for E stability are naturally realized in stoichiometric or neutron-irradiated SiO2, and conclusively put on firm ground the identification of the E center.

Our model of the stability of E is based on two native defects: the oxygen vacancy VO and the oxygen interstitial Oi. The motivation is that E is observed chiefly (though not only) in neutron-irradiated material , where VO and Oi are both abundant. Indeed, we find that it is their concurrent presence that produces the conditions for the existence of E, in neutron-irradiated as well as non-neutron-irradiated material. In the former, vacancies VO and interstitials Oi are essentially produced in pairs by knock-on, kick-out events; in the latter, they form in thermal equilibrium and, as it turns out, in similar concentrations. As will become apparent below, our argument on E applies to both cases without modifications, except for the fact that neutron-produced defects appear in concentrations determined by the irradiation dose, whereas the concentration of thermally formed defects depends on formation energies, which can be directly predicted.

At a growth temperature Tg and with Ns available sites, the equilibrium concentration of a defect is D$`=\mathrm{N}_s\mathrm{exp}(\mathrm{F}_{\mathrm{form}}/\mathrm{k}_\mathrm{B}\mathrm{T}_\mathrm{g})`$. The formation free energy Fform = Eform – T Sform depends  on the chemical potentials of atoms added or removed, on the defect charge state, i.e. the charge released to or captured from the thermodynamic reservor constituted by the surrounding crystal, and on the electronic chemical potential $`\mu _e`$ of the latter. Once the formation energies of the all relevant defects (vacancy and interstitial in our case) are known, the defect concentrations and the chemical potential $`\mu _e`$ are determined self-consistently, subject to charge neutrality, as detailed in . A specific defect configuration or charge state is predicted to exist if its formation energy is lower than that of all other defect states for some value of $`\mu _e`$. Also, the defect is metastable if a non-zero energetic barrier prevents its deactivation or disappearance into other lower-energy configurations of the same defect, or recombination with other defects. The formation energy for our defects in charge state $`Q`$ reads

$$\mathrm{E}_{\mathrm{form}}(Q)=\mathrm{E}_{\mathrm{tot}}^{\mathrm{def}}(Q)\mathrm{E}_{\mathrm{tot}}^{\mathrm{undef}}+Q\mu _e+\mathrm{M}(Q)+\mathrm{P}$$
(1)
where E$`{}_{}{}^{\mathrm{def}}{}_{\mathrm{tot}}{}^{}`$ and E$`{}_{}{}^{\mathrm{undef}}{}_{\mathrm{tot}}{}^{}`$ are the total energies of the defected and undefected system, respectively, $`\mu _e`$ is the electron chemical potential (equaling the Fermi level EF in our T=0 calculations), M$`(Q)`$ is the defect-dependent multipole correction for charge state $`Q`$ of Ref. , P = $`\mu _\mathrm{O}`$ for VO and P = –$`\mu _\mathrm{O}`$ for Oi, and $`\mu _\mathrm{O}`$ is the oxygen chemical potential. The latter is fixed to stoichiometric conditions, i.e. at the center of its variation range $`\mu _{\mathrm{mol}}/2+\mathrm{\Delta }\mathrm{H}/2<\mu _\mathrm{O}<\mu _{\mathrm{mol}}/2`$ determined by the total energy $`\mu _{\mathrm{mol}}`$ of the O2 molecule, and the calculated formation enthalpy $`\mathrm{\Delta }\mathrm{H}`$ of SiO2. Ionization energies, i.e. the energy needed to promote (say) an electron from the valence band to an empty acceptor level, are defined via total energy differences of different charge states. Formation entropies are beyond the scope of the method used here; plausible estimates are used when needed.

Energies and forces are accurately calculated from first-principles within density-functional theory in the local approximation, using the ultrasoft pseudopotential plane-wave method as implemented in VASP . An isolated defect is simulated in periodic boundary conditions via the repeated supercell approach. We use crystalline $`\alpha `$-quartz SiO2 supercells of tetragonal symmetry, comprising 71 to 73 atoms and of linear dimensions 18.49, 16.02, and 20.44 atomic units (theoretical lattice parameters , matching experiment to about 0.5%). Atomic geometries of the defects are optimized for each $`Q`$ (obtained by removing or adding electrons as appropriate) until all residual force components in the system are below 0.02 eV/Å. No symmetry restriction is imposed on geometry optimization. Improving slightly on the setting of Ref., a (222) mesh is used  for k-space summation (4 special points in the supercell Brillouin zone). Our use of a crystalline SiO2-based model of the E defect, which is observed both in amorphous and crystalline phases, is justified by the closely similar behavior of several E variants in crystalline and amorphous SiO2 in experiment  as well as theory . In addition, the simulated structure of amorphous silica  deviates moderately from that of crystalline a-quartz SiO2.

Fig.1 shows the ground state formation energies of the interstitial Oi and the vacancy VO. The vacancy acts as a double acceptor (double donor) in extreme $`n`$ ($`p`$) doping conditions, but it is in fact neutral for most of the Fermi level range. As to the interstitial, we find that an oxygen atom initially placed near the center of the hexagonal channels of quartz, relaxes sideways towards the helical chains and, after overcoming a small ($``$ 0.2 eV) barrier, it stabilizes into a split-interstitial (s-Oi) configuration with the nearest bridging oxygen in the helical tetrahedra chain, with $``$ 1 eV of energy gain with respect to the starting site. Details will be discussed elsewhere: here we note the main consequences of this result: a) the stabilization in the s-Oi configuration (which is found to be quite close to that suggested in Ref.) prevents Oi from easily recombining with vacancies, because the detachment from the split-interstitial configuration costs about 1.2 eV; b) as seen in Fig.1, s-Oi is a (negative-$`U`$) deep double acceptor with second ionization energy at 3.3 eV above the valence band top Ev. The Fermi level, calculated as in , is pinned at EF = Ev \+ 3.3 eV (vertical solid line in Fig.1). Thus a consequence of vacancy-interstitial pair formation, is that moderately $`p`$-type conditions are achieved; this is indeed often observed in irradiated samples . In the absence of s-Oi, the Fermi level would be at midgap, Ev \+ 4.4 eV. From the formation energies one can estimate the chemical concentrations of s-Oi and VO: for a typical Tg of 1500 K, and assuming a reasonable formation entropy Sform= 5 kB, the concentration of both defects is $``$1014 cm-3. This figure matches well E concentrations measured  in non-neutron irradiated samples after UV, $`\gamma `$, or X illumination fall, and therefore supports the hypothesis that the vacancy is the parent defect of E. Neutron irradiation of course produces dose-dependent , typically much higher concentrations ($`>10^{19}`$ cm-3). The key point is that, since the concentrations of vacancies and interstitial are essentially the same in both cases, and because only the interstitial has electrical activity (via its double acceptor level), the Fermi level is pinned at the same value in both cases, so that the theory of E stability and activation discussed below applies identically to both cases.

It appears from Fig.1 that the +1 charge state of the vacancy, V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$, is not among the thermodynamically stable ground states of the defect. Therefore, if this state of VO is to be identified with E, it must at least be proven metastable; if it is metastable, a mechanism for its creation starting from the ground state (the neutral vacancy) must be identified. As to the first point, since the +1 vacancy has a formation energy that increases linearly with the Fermi level while the neutral one remains constant, V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$ or a distorted variant thereof, may only be metastable in a limited range of that variable: we show below that the Fermi level pinning due to split-Oi produces naturally the conditions for the metastable existence of E. As to the second point, experiments indicate that E is activated by ionizing radiation  such as $`\gamma `$, X, or UV photons shone onto vacancy-containing samples, or concurrently with neutron irradiation (causing knock-on vacancy creation). Indeed, since the vacancy ground state is the neutral undistorted configuration, the distorted +1 state (alias E) can only be accessed by ionization of V$`{}_{}{}^{0}{}_{\mathrm{O}}{}^{}`$: our picture provides consistently such activation mechanism.

We proceed to study the behavior of the vacancy when subjected to the undistorted-to-puckered transition as proposed in earlier studies . In accord with the results discussed above, we fix the Fermi level at Ev \+ 3.3 eV. The creation of a vacancy starting from the perfect lattice results in moderate local distortions in both the neutral and +1 charge states. The puckered configuration is obtained by moving one of the two vacancy-adjacent Si1 and Si2 atoms (specifically the “long-bond  Si2) away from the vacant site, and pushing it across the basal plane of the incomplete tetrahedron centered on Si2 itself. When Si2 pokes through this triangular constriction, it gets strongly and suddenly bound to a backbonding oxygen, Ob. Upon completion of the distortion, Si2 regains 4-fold coordination, and backbonding Ob becomes 3-fold coordinated (see also ), while it was originally 2-fold coordinated as all tetrahedron-bridging oxygens in SiO2. Si1 remains instead 3-fold coordinated: in the +1 charge state, its dangling bond is half-filled, and causes the observed  and predicted  EPR signature which identifies E .

It may appear at first sight that the puckering distortion should be symmetric in Si1 and Si2. This is not the case, however, because of the intrinsic asymmetry of the quartz structure. As already noted earlier on , only Si2 finds the backbonding oxygen Ob in the correct position. This applies largely also to amorphous silica, whose structure is moderately different at the local level from that of quartz . The backbonding oxygen is therefore the main stabilizing agent of the E defect.

The total energy of the system in charge state $`Q`$ is calculated as a function of the separation between Si1 and Si2. Only the modulus dSi-Si of the Si1-Si2 connecting vector is constrained, and all other degrees of freedom are fully relaxed: the minimum energy path is thus mapped out for the undistorted-puckered transformation in the constrained-dSi-Si configurational subspace. In Fig.2 we display the full level diagram for the neutral, +1, and +2 charge states of the vacancy as a function of the puckering distortion, quantified by dSi-Si. All energy curves depend on the Fermi level through Eq.1; they can be directly compared on the same energy scale because EF is fixed at the natural value of Ev \+ 3.3 eV determined above.

The outstanding feature of Fig.2 is that for the natural Fermi level of stoichiometric or neutron-irradiated silica, the candidate E, i.e. distorted V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$, is indeed the stable defect state for the distorted geometry. We stress that the Fermi level position is essential here: if EF were at midgap, the +1 curve would be 1.1 eV above its position in Fig.2. Then E would be unstable towards magnetically-inactive V$`{}_{}{}^{0}{}_{\mathrm{O}}{}^{}`$. Globally, E is metastable with a confining barrier of 0.8 eV. The barrier to enter the metastable state is 1.1 eV, and the undistorted state is lower than the distorted by 0.3 eV \[this difference to Ref. , where the distorted state was lower by the same amount, is possibly due to our improved k-sampling\]. Clearly, in the absence of excitation, V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$ will remain trapped in the metastable E state and will show EPR activity. When thermally activated to overcome the barrier, the puckered center will transform into undistorted, and by electron capture it will become neutral. Therewith, E disappears permanently, and so does its magnetic activity, because Si1 and Si2 combine their dangling bonds to bind into a dimer . (The same deactivation route is not readily available for the level ordering of Ref. , which implies that i) E remains activated at equilibrium, and that ii) concurrent barrier jump and electron capture are required to quench it.)

The calculated puckered-to-undistorted barrier is now compared with an estimate extracted from the measured relative drop in E population upon isocronal thermal annealing  in irradiated silica. In the simplest model, the distorted-state population N0 diminishes by a factor of N/N0 = exp (–R(T) $`\tau `$) upon annealing over a time $`\tau `$ at temperature T, with escape rate R = $`\nu `$ exp (-$`\mathrm{\Delta }`$F/kBT), with $`\mathrm{\Delta }`$F the free energy barrier for escape from the puckered state, and $`\nu `$ an average vibrational frequency in that state. Using the data  for E$`{}_{\gamma }{}^{}{}_{}{}^{}`$ and assuming $`\nu 50`$ THz, we obtain an experimental $`\mathrm{\Delta }`$F$``$1.1 eV, in reasonable agreement with our calculated $`\mathrm{\Delta }`$F$`\mathrm{\Delta }`$E=0.8 eV; account for the transition entropy, which we neglect, should further improve agreement since the transition state (through the tetrahedron basal plane) is severely constrained geometrically, and has a higher average vibrational frequency.

Let us now come to E activation. In the present picture, E is created transforming undistorted V$`{}_{}{}^{0}{}_{\mathrm{O}}{}^{}`$ into distorted V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$ via two routes. The first proceeds on path A in Fig.2 with two successive one-photon ionizations of V$`{}_{}{}^{0}{}_{\mathrm{O}}{}^{}`$ into V$`{}_{}{}^{+2}{}_{\mathrm{O}}{}^{}`$, followed by non-radiative decay into E. This path is efficient since the +1 undistorted state is kept populated by sustained illumination (a much less efficient double-photon excitation of V$`{}_{}{}^{0}{}_{\mathrm{O}}{}^{}`$ into V$`{}_{}{}^{+2}{}_{\mathrm{O}}{}^{}`$ may also occur). The excitation energies for path A are both near 4–4.5 eV if the ionized electron is transfered to the Fermi level, i.e. to the EF-pinning impurity; if it is promoted to the conduction band, the energies are instead about 7 eV. Both processes are possible with X or $`\gamma `$ radiation, whose energy vastly exceeds that needed in the transition. In UV irradiation, the center is often activated by pumping at 5 eV, and clearly only Fermi-level capture matches this figure. (There are, however, qualitative indications that the 5 eV excitation may activate E via alternate routes involving other pre-existing defects.)

The second excitation route, path B, involves an optical excitation of V$`{}_{}{}^{0}{}_{\mathrm{O}}{}^{}`$ into undistorted V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$, and a thermal excitation of the latter into the puckered state. The energy difference (0.3 eV) between the two V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$ states implies that only a fraction of 10-5 of the vacancies gets promoted into the distorted state in equilibrium at room temperature, on sustained illumination. Therefore, though admissible, this path is preempted by path A. With the level ordering of Ref. , also path B competes with path A. (Our ordering, however, matches better the thermal behavior, as discussed above.)

Other calculated observables of our E model are consistent with previous studies . For instance, the optical absorption of the neutral undistorted state at 6.9 eV correlates well with the 7.6 eV absorption band usually attributed to the neutral vacancy . For E (metastable puckered V$`{}_{}{}^{+}{}_{\mathrm{O}}{}^{}`$) we find an absorption at 4.7 eV (defect-to-conduction promotion) followed by slow non-radiative decay back into E. Since it is not followed by any emission or E deactivation, this absorption correlates with the broad 5.8-eV band typical of E , which exhibits the same behavior in experiment.

In summary, we conclusively put on firm ground the identification of the singly positive O vacancy in SiO2 with the E center proving its thermodynamical stability via first principles calculations. Our picture provides naturally activation and deactivation mechanisms, and other optical signatures, in agreement with known experimental observations. In addition, our picture naturally explains the moderate $`p`$ conditions produced by irradiation in SiO2.

We acknowledge support from Istituto Nazionale per la Fisica della Materia under “Iniziativa Trasversale Calcolo Parallelo”.