text:
# The DEEP2 Galaxy Redshift Survey: Spectral classification of galaxies at 𝑧∼1

## 1. Introduction

The classification of galaxies is of fundamental importance for understanding galaxy populations, and for this reason is a very important aspect of any galaxy redshift survey. Having a data set of many thousands of galaxy spectra allows one to test the validity of galaxy formation and evolution scenarios with unprecedented accuracy. However, the sheer size of the full spectral data set presents its own unique problems. In order to make such a galaxy data set more ‘digestible’ some form of data compression is necessary, whether this be through the adoption of morphological segregation, colors or some other compression/classification scheme. If these quantities (and their associated distributions) can be determined consistently over a wide range of redshifts, they can be compared with theoretical predictions and simulations, and hence set constraints on scenarios for galaxy evolution. This will be especially true if consistent classification regimes can be used for both the high redshift ($`z1`$) surveys currently underway (the DEEP2 Galaxy Redshift Survey, Davis et al. 2002 and the VIRMOS-VLT Deep Survey, Le Fevre et al. 1999) – and the large $`z0`$ surveys now approaching full completion; the Sloan Digital Sky Survey (SDSS, Strauss et al. 2002) and the 2dF Galaxy Redshift Survey (2dFGRS, Colless et al. 2001)

A number of different approaches to the classification of galaxy spectra have been adopted for local galaxy surveys. These include the calculation of rest-frame colors (e.g. Strateva et al. 2001); principal component analysis (PCA) based spectral classifications (e.g. Connolly et al. 1995; Bromley et al. 1998; Folkes et al. 1999; Madgwick et al. 2002; de Lapparent et al. 2003), and other more sophisticated discriminations (e.g. Heavens, Jimenez & Lahav 2000; Slonim et al. 2001), based upon information theory. The underlying theme of all these alternative methods is that they characterize the galaxy population exclusively in terms of their observed spectra. The work presented here is the first attempt to apply one of these methods (PCA) to the classification of such a distant sample of galaxies.

### 1.1. The role of spectral classification

There are three methods which have generally proved to be popular for the classification of galaxies: morphological segregation, rest-frame colors and direct spectrum based classifications. Each of these methods has its own unique drawbacks and advantages.

To understand galaxy evolution out to redshifts of $`z1`$, it is essential to have a consistent implementation of these classifications over a wide range of look-back times. For this reason it can be argued that morphological segregation – although perhaps the most natural form of classification – may not be the optimal solution to adopt over such a large range of redshifts. This is due to both the degradation of morphology with redshift and the absence of a robust and repeatable methodology to perform this classification (see e.g. Conselice 2003 for further discussion and possible solutions to this situation). For this reason we focus in this paper on alternative classification methods to morphology, which we hope will complement earlier studies based on this method, whilst at the same time providing a new perspective by more directly reflecting the physical properties of each galaxy.

The remaining two options – rest-frame colors and spectral ‘types’ – are linked, in that they both provide some compressed representation of the observed galaxy’s spectral energy distribution (SED), hence providing a relatively direct insight into physical processes such as star-formation currently occurring in each galaxy. However, in terms of how each is calculated for high-z galaxies there are significant differences which are important to understand.

The main complication for the calculation of rest-frame colors is the accurate determination of $`k`$-corrections to account for the different pass-bands sampled by each photometric filter at different redshifts. These generally cannot be estimated from the observed spectrum, but rather must be determined by matching each observed galaxy to a set of template SEDs with full rest-frame wavelength coverage.

In the case of spectral classification, one must consider that the rest-frame wavelength coverage varies with the redshift of the galaxy, hence to adopt a uniform classification over a large range of redshifts one needs to significantly restrict the rest-frame wavelength range considered (e.g. by focusing in on particular line features through equivalent width measurements) or to determine some way of filling in the gaps that are not observed in each spectrum.

These two problems are very closely related in that they express the need to determine the form of a galaxy’s SED over a wavelength interval that is not necessarily observed. In the case of $`k`$-corrections only $`10`$ template SEDs are available for this task (Kinney et al. 1996; Coleman, Wu & Weedman 1980); reddening is manually incorporated and we have little control over the wavelength intervals involved. However as will be shown in this paper, with principal component analysis (PCA), we can achieve this task for spectral classification, using the observed spectra themselves as templates, giving us $`N_{\mathrm{gal}}`$ galaxy templates, with which to interpolate or extrapolate a given observed spectrum. In addition, because we have much more control over the wavelength interval adopted for the classification we can modify the analysis to use only the parts of the observed spectra which have the most uniform rest-frame wavelength coverage.

It is primarily for this latter reason that adopting PCA-based spectral types will provide a classification scheme that is particularly uniform over the large range of redshifts encountered in the DEEP2 redshift survey, hence providing a robust probe of evolution in the galaxy population. In addition we note that such methods of classification are timely, in that there is now a considerable body of low redshift spectroscopic data from e.g. 2dFGRS and SDSS, which enables us to make detailed comparisons. We note that PCA-based classifications do suffer from one particular drawback, which is that they are not as straightforward to interpret as other classifications, and this is an issue we will attempt to address in this paper.

The outline of this paper is as follows: In Section 2 we briefly discuss the DEEP2 Redshift Survey data that we will use in this analysis. Section 3 describes the implementation of PCA we have adopted for this paper, and gives a discussion as to why this particular formalism has been used. In particular, issues regarding the restframe wavelength coverage are discussed in some detail. In Section 4 we discuss the method of spectral classification we will adopt, based upon the results of the PCA we have implemented. This spectral classification is contrasted with that of the 2dF Galaxy Redshift Survey in Section 5, in which the selection effects of the two surveys are also discussed. We then conclude this paper in Section 6 with a brief discussion as to future applications of this work.

## 2. DEEP2 galaxy spectra

In its first season of observations (August – October, 2002), the DEEP2 Redshift Survey has already accurately measured the redshifts of $`5600`$ galaxies, out of a proposed total of 60,000. These galaxies have been pre-selected to have $`z0.7`$ and an $`R_{\mathrm{AB}}`$ limiting magnitude of 24.1, from a set of $`B`$,$`R`$ and $`I`$ CFHT 12k x 8k images covering $`4`$ deg2 on the sky. Foreground galaxies have been excluded using a simple photometric cut, based upon the observed $`RI`$ and $`BR`$ color of each galaxy (Davis et al. 2002). In this paper we make use of all this data from the first observing season.

A sophisticated automated pipeline has been developed to efficiently extract and reduce the spectra observed in the survey, details of which will be presented by Davis et al. (2003). The observed spectra themselves are taken at high resolution ($`R5000`$) using the DEIMOS spectrograph mounted on the 10-m Keck-II telescope (Davis & Faber, 1998), and generally span the wavelength range of $`6400<\lambda <9200`$Å. For galaxies with $`z>0.7`$ this allows us to very accurately measure the redshift of each galaxy, particularly when the resolved \[Oii\] doublet is present (out to a redshift of $`z=1.5`$). Absorption based redshifts are also readily determined, primarily using the Ca H+K features which are visible to redshifts of $`z1.3`$.

All spectra used in this analysis have been corrected for the DEIMOS throughput efficiency of the 1200 line/mm grating used in the survey111see http://www.ucolick.org/ripisc/Go1200.html, and are presented in units of counts/pixel/hour. Note that the flux calibration is only approximate at present ($``$10%), however we have found the components derived from our PCA to be robust to the exact value of the throughput efficiency. Additionally, the spectra have been normalized to have a mean count of unity, after being minimally smoothed (3 pixels $`0.9`$Å observed-frame, $`0.5`$Å restframe) according to the inverse of their estimated variance (in order to remove sky residuals and other artifacts remaining from the spectral reduction).

## 3. Principal Component Analysis

The spectral classification presented here is based upon a Principal Component Analysis (PCA) of the DEEP2 galaxy spectra. PCA is a powerful technique, allowing us to easily visualize and quantify a multi-dimensional population in terms of just a handful of significant components. It does this by identifying the components of a data set (in this case the galaxy spectra) which are the most discriminatory between each galaxy (where the significance is determined in terms of its contribution to the variance over the entire sample). This allows us to identify just the most significant components for future use. It is clear from such a formalism that any clustering in the PCA space is indicative of distinct sub-populations within the sample.

PCA has been used with considerable success by a variety of authors to deal with large multi-dimensional data sets. Several similar mathematical formulations of PCA for galaxy spectra have appeared in the literature, in particular Connolly et al. (1995); Bromley et al. (1998); Galaz & de Lapparent (1998); Folkes et al. (1999). The most significant difference between the various techniques is that some methods utilize mean-subtracted covariance matrices for the PCA (by subtracting the mean of the normalized spectrum of each galaxy), while others do not. This has little effect on the subsequent analysis since the latter methods simply yield a mean spectrum as the first component.

In this paper we present a new variation on previous formalisms for carrying out the PCA of galaxy spectra, designed to accommodate the features of our data set that make the standard method very difficult to apply. In particular there are two complications that must be addressed to analyze the DEEP2 galaxy spectra: the first is the very high resolution of the observed spectra which require extremely time-consuming matrix computations. The second is that the observed spectra cover very different rest-frame wavelength ranges once de-redshifted, and also have a number of effective ’gaps’ present due to the presence of strong night-sky lines. For this reason a method of performing interpolations and (small) extrapolations is required to ensure the classification is uniform.

### 3.1. Implementing PCA for high-resolution data

A significant drawback to implementing PCA on large or very high-dimensional data sets is the required computation time, particularly to determine the covariance matrix. For $`n`$ galaxies - each with $`p`$ spectral channels - this requires $`O(np^2)`$ operations. Given that each DEEP2 spectrum contains $`O(10^4)`$ channels this would be very time consuming indeed if the full spectral resolution is used.

Fortunately, it is possible to solve for the PCA eigenspectra without calculating the entire covariance matrix. The key is developing a probabilistic formalism for the PCA, compatible with an expectation-maximization (E-M) algorithm (see Roweis 1997 and Tipping & Bishop 1999). Adopting such a formalism allows one to solve for the first $`k`$ eigenspectra in only $`O(npk)`$ operations. It has been shown that this method for performing PCA is robust, in that it has only one stationary point that is not a saddle point, which guarantees there will be no false convergences (Tipping & Bishop 1999).

Computationally, the E-M algorithm proceeds in iterations of two steps. The $`k`$ eigenvectors to be calculated are assumed to be spanned by the columns of a $`p\times k`$ matrix $`𝐂`$. We start by making an initial (random) guess for the columns of this matrix and use this to determine the $`k\times n`$ matrix, $`𝐗`$, of $`k`$ ‘states’ for each galaxy. These states correspond to the principal component projections of each galaxy in the non-orthogonal space defined by the columns of $`𝐂`$. These states are then used in conjunction with the original $`p\times n`$ data matrix, $`𝐘`$, to make a better estimate for $`𝐂`$. This proceeds until convergence. These steps can be summarized as,

1. Step I:

$$𝐗=(𝐂^T𝐂)^1𝐂^T𝐘$$
2. Step II:

$$𝐂^{\mathrm{new}}=\mathrm{𝐘𝐗}^T(\mathrm{𝐗𝐗}^T)^1$$

Once the algorithm converges the columns of $`𝐂`$ will span the space of the first $`k`$ eigenvectors. This can therefore be used to construct the orthogonal basis that defines the principal components and their projections. We note that this method for PCA gives identical principal components to the other (simpler) formalisms discussed previously. The sole reason we adopt an E-M based approach is to make the PCA efficient for the large and very high resolution data set we are using. Throughout the remainder of this paper we will denote the eigenvectors determined by the PCA (herein eigenspectra) as $`𝐏_\mathrm{𝟏}`$,$`𝐏_\mathrm{𝟐}`$ etc. and the projections onto these axes by $`p_1`$, $`p_2`$ etc.

### 3.2. Dealing with incomplete data in the PCA

The treatment of incomplete data in PCA has already been discussed in the literature, e.g. Everson & Sirovich (1995) and Connolly & Szalay (1999). However, the implementation is complicated and deserves further discussion. The issues that must be resolved are two-fold: First, how to determine the eigenspectra (or principal components) when very few galaxies cover the full wavelength range, and second, how to project a galaxy onto these eigenspectra when its spectrum does not cover the entire restframe wavelength range considered. It can be shown that the latter issue has a relatively straightforward, clean solution involving de-correlating the eigenspectra (which are not orthogonal when we do not use their entire wavelength coverage). This has been presented in Connolly & Szalay (1999). However, estimating the eigenspectra in the first place is much more difficult to address.

Consider the situation when the eigenspectra, $`𝐏_i`$, have already been estimated using the PCA. Each of the observed, mean-subtracted spectra, $`𝐟^{}=𝐟\overline{𝐟}`$, can then be expressed as a linear combination of these eigenspectra,

$$𝐟^{}=\underset{i}{}p_i𝐏_i,$$
(1)
where $`\overline{𝐟}`$ is the mean spectrum of all the galaxies and $`p_i=𝐟^{}𝐏_i`$ is the projection of each galaxy onto the $`i`$th principal component. The index $`i`$ ranges from 1 to the number of resolution elements in each spectrum. However, because most of the variance in each spectrum is contained in only the first few eigenspectra (as will be demonstrated later), the sum in Eqn. 1 can be truncated to include only the first $`k`$ most significant elements.

We wish to determine how to derive the projections $`p_i`$ when the spectrum, $`𝐟^{}`$, does not cover the full restframe wavelength range spanned by the eigenspectra (e.g. because of poor sky subtraction or the effects of redshift). The eigenspectra are defined to be orthonormal,

$$𝐏_i𝐏_j=\delta _{ij}.$$
(2)
However, when we only consider the wavelength range over which the observed spectrum, $`𝐟^{}`$, is defined this is no longer the case. We denote the incomplete spectrum by $`\mathrm{𝐰𝐟}^{}`$, where $`𝐰`$ is a vector with non-zero elements only where $`𝐟^{}`$ is defined. As shown in Connolly & Szalay (1999), the correlations between the portions of the eigenspectra relevant for $`f^{}`$ can be quantified in terms of a correlation matrix,

$$𝐌_{ij}=\mathrm{𝐰𝐏}_i𝐏_j.$$
(3)
Using this matrix it is possible to show that, when a spectrum does not cover the entire wavelength range of the eigenspectra, unbiased estimates of the true projections can be calculated simply using the inverse of this correlation matrix, $`M`$,

$$\widehat{p}_i=\underset{j}{}𝐌_{ij}^1𝐏_j,$$
(4)
where again this sum need only be carried out over the most significant principal components (in our case we will use the first $`k=20`$). In what follows we use this procedure for two separate purposes: to estimate the corrected projections, $`\widehat{p}_i`$, when the wavelength range of the galaxy does not fully span the eigenspectra, and (using these corrected projections) to interpolate over any gaps in a spectrum resulting from sky-subtraction etc., using Eqn. 1. This is discussed further in the next section.

Before this formalism can be adopted one first requires an estimate of the eigenvectors of the PCA, $`𝐏_i`$. Several approaches to estimating these eigenspectra in cases where few or no galaxies cover the full wavelength range of interest have been suggested. Everson & Sirovich (1995) discuss approaches involving extrapolation of the observed spectra followed by iteration. Another approach is to make a least-squares generalization to the E-M PCA (Roweis 1997). Each of these methods assumes that the ‘gaps’ in each spectrum are randomly positioned, which may be true of low-z galaxy spectra but is clearly unsatisfactory for high-z surveys since almost all of the observed spectra will miss a substantial amount of either the extreme blue or extreme red end of the full rest-frame wavelength range. For this reason a simple approach to estimating the eigenspectra will normally fail, since it is difficult to correlate (and hence extrapolate) one end of the observed spectrum with what would be expected at the other end.

A more practical solution, which we follow in this paper, is to restrict the rest-frame wavelength range of interest such that there is a more significant subsample of galaxies which span the entire wavelength range of interest. Two such divisions are outlined in Fig. 1; the first involves considering only the rest-frame wavelength range of 3700-4200Å (covering \[Oii\] and beyond the 4000Å break) while the second is extended to 3700-5050 Å (essentially covering \[Oii\] to \[Oiii\]). There is an obvious draw-back to this approach in that we are not making use of all the information contained in our observed spectra; however, as will be shown below this is a necessary compromise.

### 3.3. Limited $`\lambda `$ PCA

Using the first wavelength cut (3700-4200Å), $`50`$% of our galaxy sample span the entire rest-frame wavelength range (with the exception of small gaps masked out due to bad sky subtraction, gaps between the CCDs and bad pixels) and can be used in a PCA analysis. For the second wavelength cut (3700-5050Å), we are using a greater portion of the observed galaxy spectra and so the analysis is potentially more informative; however, only $`10`$% of galaxies span this entire wavelength range.

In each case the procedure we adopt is as follows:

1. The subset of our galaxies which span the entire restframe wavelength range of interest are passed through the E-M PCA to determine a first set of $`k=20`$ eigenspectra, corresponding to the orthogonalized columns of the matrix $`𝐂`$ (c.f. Section 3.1).
2. Gaps in the galaxy spectra used in the previous step are interpolated over by re-calculating the ‘de-correlated’ projections, $`\widehat{p}_i`$ (as described in the previous section) and are then used to reconstruct the missing parts of the spectrum. The new galaxy spectra are then re-processed through the PCA to produce improved estimates of the PCA eigenspectra.
3. The new set of eigenspectra is then used to compute projections for the entire galaxy sample regardless of rest-frame wavelength coverage, by means of the de-correlation procedure outlined in Eqn. 4.

We do not iterate further after projections have been determined for the entire galaxy sample, since this tends to propagate errors through the analysis and can lead to unphysical eigenspectra.

The values of the first two projections derived in this manner are plotted in Fig. 2, for both the rest-frame wavelength ranges considered. This figure allows us to determine whether the PCA is robust beyond the narrow redshift range that corresponds to each of the rest-frame wavelength limits adopted. To judge this we have compared in each case the distribution of the first two components ($`p_1`$ and $`p_2`$) for only those galaxies used in the initial PCA (the restricted $`z`$ sample), and those for the entire galaxy sample (the unrestricted $`z`$ sample). Although some small degree of systematic offset is expected, the two should be broadly similar, which is clearly not the case for the second sample used (3700-5050Å). This lack of agreement effectively highlights the limitations of performing this PCA formalism for spectra that do not have a sufficient degree of overlap.

Clearly the extended sample is more interesting, in so far as it covers a larger rest-frame wavelength range and therefore includes more spectral features. However, Fig. 2 demonstrates that this classification cannot in fact be performed robustly since too few galaxies cover this full rest-frame wavelength range. On the other hand, the more restricted wavelength choice (3700-4200Å) is much more robust. We therefore choose to adopt this latter wavelength coverage to define our spectral classification.

## 4. Spectral Classification

PCA is not a classification algorithm, rather, it is simply a method of data compression. For this reason another (usually manual) step must be performed, in which the insight gained from this compression is used to divide a galaxy sample. This is not necessarily straightforward and will usually involve some degree of arbitrariness – especially since galaxy spectra appear to span a single continuum of types rather than falling into distinct classes.

By far the most significant output of the PCA is the first principal component ($`p_1`$ contains 10 times more variance than any other component in our analysis, see Fig. 3), which appears to be dominated by the nebular emission line strengths in a spectrum. In particular the strength of the \[Oii\] emission feature figures prominently, as shown in Fig. 4. Note however that, as demonstrated in Fig. 3, this component does in fact also encapsulate a great deal more information about each spectrum and can be used to very accurately reconstruct a wide variety of galaxy spectra. In particular, because the Ca H&K absorption features are inverted in this spectrum, it can act as a classifier for both absorption and emission dominated galaxy spectra.

Other eigenspectra, examples of which are shown in Fig. 4, also contain useful information about each spectrum. For example, the second eigenspectrum, $`𝐏_\mathrm{𝟐}`$, appears to primarily quantify the width of the various emission lines that is not incorporated into the first principal component. The component $`𝐏_\mathrm{𝟑}`$ reflects the stellar continuum, and in particular the height of the 4000Å break relative to \[Oii\] that is not already encapsulated in the first principal component. Note that each of these eigenspectra represent variations in the galaxy spectra relative to the mean spectrum and so, for example, the \[Oii\] line in $`𝐏_\mathrm{𝟏}`$ represents the additional emission line strength present in each spectrum, whether this be less than ($`p_1<0`$) or greater than ($`p_1>0`$) the value for the mean spectrum. Beyond the fourth principal component the eigenspectra become dominated by unphysical features in the galaxy spectra.

Because the first component is so much more discriminatory between galaxies (in terms of variance) this appears to be the most logical (and simplest) classification to adopt, particularly since this component alone appears to do well at reconstructing a large variety of observed spectra. It is additionally reassuring to note that $`𝐏_\mathrm{𝟏}`$ also encapsulates the greatest variety of features in each galaxy’s spectrum (Fig. 4), and so is also potentially the most astrophysically interesting component available from the PCA. Note that excluding the second component means that we are not making use of the detailed line-width information, some of which is contained in $`p_2`$. As the DEEP2 reduction pipeline improves more precise linewidth information, and even detailed rotation curves, will be extracted from the full 2-dimensional spectrum of each galaxy (rather than the compressed 1-dimensional spectra used here). In so doing we will be able to provide a much more accurate characterization of the kinematic properties of each galaxy that will be complimentary to the relative emission/absorption line-strengths measured by $`p_1`$.

A similar dominance of the first principal component, $`p_1`$, has been noted in previous analyses of low-z galaxy surveys, suggesting that the distribution of ‘normal’ galaxy spectra can be well approximated by a one dimensional sequence. However, we must be careful to ensure that this projection is in fact a robust measure of a galaxy’s spectrum if we are to use it alone to determine spectral type. For example, it has been noted for low-z galaxy samples observed through small fiber optic apertures that the first principal component was not stable between multiple observations of the same object (Bromley et al. 1998; Madgwick et al. 2002). We find that this is not the case for our most significant projections, which proved to have relatively small measurement uncertainties. In particular, we measured the value of $`p_1`$ obtained from the instances where a given galaxy spectrum was multiply observed and found that the distribution of errors was gaussian with a standard deviation of only $``$5. This is most likely owing to the use of slit-apertures instead of fiber optics. We therefore choose to adopt $`p_1`$ as our measure of spectral type for galaxies observed in the DEEP2 Redshift Survey and denote it by,

$$\eta _{\mathrm{DEEP}}p_1.$$
(5)
This notation has been specifically chosen to analogous to low redshift samples (e.g. the 2dFGRS) for which a similar classification scheme has been adopted.222Note that in the 2dFGRS, $`\eta `$ was defined as a linear combination of the first two components; $`p_1`$ and $`p_2`$. The two components were needed because of the additional calibration uncertainties introduced by the small fiber apertures used in that survey. Despite this the classifications are in fact comparable as that combination was also the most statistically significant output of the PCA which was robust to the known instrumental uncertainties.

A histogram of $`\eta `$ is shown in Fig. 5, in which a significant bimodality appears to be present. Since this is the only discontinuity in the distribution of this projection, it is natural to divide the sample there ($`\eta 13`$). The presence of this bimodality is interesting in that similar effects have been observed with galaxy colors (Strateva et al. 2001; Bell et al. 2003; Weiner et al. 2003) and also in spectral classifications of other samples such as in the 2dFGRS (Madgwick et al. 2002).

Having divided the galaxies into two ‘types’, we proceed to calculate the mean spectrum for each type, shown in Fig. 6. Clearly, the difference between these two classes is quite pronounced, and is differentiating primarily between absorption spectra (associated with older stellar populations) and emission dominated spectra (associated with recent star-formation activity). Taking this interpretation further, it should be possible to directly relate this classification to the underlying star-formation activity in each galaxy as was done for the 2dFGRS (see Madgwick et al. 2003). We return to this point in Sec. 5. Galaxies with values of $`\eta `$ less than $`13`$ will be referred to as either ‘early type’ or passively evolving galaxies, whilst those with higher values will be referred to as ‘late type’ or actively star-forming galaxies.

Figure 7 compares the the restframe $`(BR)_0`$ colors of each galaxy with its spectral type, as defined by our $`\eta `$ parameter. It can be seen from this figure that there is a correlation between the restframe color and the strength of emission/absorption features in each observed spectrum, particularly for those galaxies with the largest $`(BR)_0`$.

## 5. Comparison with $`z=0`$ classification

In this Section we briefly compare the classification derived for the DEEP2 Survey with that adopted in the 2dF Galaxy Redshift Survey (Colless et al. 2001), which comprises over 200,000 low redshift galaxy spectra. The reason we make this particular comparison is that the two classification regimes have been similarly defined, in that they are both derived from the most significant component of the PCA that is robust to the known instrumental uncertainties (Madgwick et al. 2002). Such a classification for the Sloan Digital Sky Survey is also in progress but is not yet available.

Both classification methods provide a continuous parameterization of the spectral type of a galaxy, denoted by $`\eta `$. At a practical level, this spectral classification is simply a dot-product between each observed galaxy spectrum and the chosen classifying eigenspectra, which are shown in Fig. 8 for the two surveys. It can be seen from this figure that despite the very different selection criteria in the two surveys, the chosen classifications both sample the same spectral features in the same relative proportions over their common wavelength interval. Hence the two classifications should be very comparable for broad studies of galaxy evolution.

The distributions of the $`\eta `$ spectral types are also shown in Fig. 8. Note however that because the DEEP2 survey spectra have much higher resolution, the locus of projections for this survey spans a much larger range (since each spectrum has been normalized to have mean flux of unity regardless of resolution). For this reason we have scaled the $`\eta _{2\mathrm{d}\mathrm{F}}`$ spectral parameter of the 2dFGRS survey to match $`\eta _{\mathrm{DEEP}}`$ in this comparison. We find that multiplying by a factor of 8 suffices to match the two classifications, a similar factor would be derived by noting that early type galaxies are separated from late types at $`\eta _{2\mathrm{d}\mathrm{F}}=1.4`$ (Madgwick et al. 2002), as opposed to $`\eta _{\mathrm{DEEP}}=13`$.

Figure 9 shows a comparison of the average spectra for early and late type galaxies as calculated from the 2dFGRS and DEEP2 spectra (with the latter smoothed to the 9Å resolution of the former). This figure again highlights that the correspondence between the two classifications is quite striking, despite being derived over different wavelength ranges and at different resolutions, they essentially encapsulate the same physical information.

It has already been shown (Madgwick et al. 2003) that the spectral classification $`\eta `$ adopted for the 2dFGRS corresponds most naturally to the relative amount of star formation activity currently occurring in each galaxy as compared with its past average (the Scalo birthrate, $`b`$, parameter, Scalo 1986). Objects with high $`\eta `$ values are galaxies with particularly strong recent star formation, whereas the lowest-$`\eta `$ sample of galaxies has $`b0.1`$ (i.e. their current star formation rate is only 10% of their past averaged value). Because there is such a simple one-to-one correspondence between this spectral classification and the star formation activity of each galaxy, we can confirm from Fig. 8 our intuition that selecting galaxies using restframe $`U`$ magnitudes (as is the case for the DEEP2 Survey, see e.g. Weiner et al. 2003) is very much more biased towards galaxies with recent episodes of star-formation than the $`B`$ selection adopted for the 2dFGRS. This result will have important repercussions for the interpretation of future analyses of the DEEP2 Survey data.

A more detailed assessment of the exact relative frequency of star forming galaxies in each survey will be forthcoming in a future paper, in which the selection function of each survey will be fully incorporated. Once this is available it will also be particularly interesting to experiment with how consistent the evolution between the populations is with that expected from spectral synthesis models (e.g. Bruzual & Charlot 1993; Fioc & Rocca-Volmerange 1997).

## 6. Conclusions

In this paper we have presented a new PCA-based spectral classification, $`\eta _{\mathrm{DEEP}}`$, for the galaxies observed to date in the DEEP2 Redshift Survey. The main goal in developing this classification was to provide a consistent and robust measure of the type of a galaxy over the large range of redshifts encountered in this survey. To do this special handling of incomplete and ‘gappy’ observed spectra is required.

Although this classification, $`\eta _{\mathrm{DEEP}}`$, appears to primarily identify galaxies with differing strengths of nebular emission, Fig. 3 demonstrates that this component alone can reconstruct the spectra of a wide variety of galaxies extremely well.

We have not directly related this classification to the role of star formation activity in a given galaxy However, there is strong evidence from a similar analysis carried out for the 2dFGRS that $`\eta _{\mathrm{DEEP}}`$ should also correlate well with the relative amount of star formation. A more detailed study of this correlation, together with the role of higher order PCA projections, will be forthcoming when higher resolution spectral synthesis models become publicly available in the near future (the DEEP2 galaxy spectra are at a resolution that is $`20`$ times higher than any synthesis model currently available).

This particular classification will be particularly useful in subsequent analyses e.g. of the galaxy luminosity function or correlation functions, as it is easily comparable with other classifications at $`z=0`$, for which large spectroscopic samples are now publicly available. In addition, previous work (e.g. Madgwick et al. 2003) has shown that it is straightforward to make direct comparisons between such a spectral classification regime and the output of semi-analytic galaxy models (e.g. Kauffmann, White & Guiderdoni 1993; Cole et al. 1994; Somerville & Primack 1999), allowing us to directly constrain the assumed models of galaxy formation and evolution between $`z=1`$ and $`z=0`$ using this form of classification.

When used in conjunction with spectral synthesis models, we expect that parameters similar to the one presented here will provide a wealth of information on both the amount and ‘type’ of evolution that has occurred in the galaxy population. This is particularly relevant now that such large samples of galaxy spectra are becoming available over such a wide range of redshifts. For this reason we expect studies of the evolution in the galaxy population to become a particularly rich field of research in the near future, and that the DEEP2 Redshift Survey will play an especially prominent role in this field.

This work was supported in part by NSF grants AST00-71048 and KDI-9872979. The DEIMOS spectrograph was funded by a grant from CARA (Keck Observatory), by an NSF Facilities and Infrastructure grant (AST92-2540), by the Center for Particle Astrophysics, and by gifts from Sun Microsystems and the Quantum Corporation. DPF is supported by a Hubble Fellowship. The DEEP2 Redshift Survey has been made possible through the dedicated efforts of the DEIMOS staff at UC Santa Cruz who built the instrument and the Keck Observatory staff who have supported it on the telescope.


format: text

text:
# Fluctuation properties of strength functions associated with giant resonances

## I Introduction

Strengths of exciting high lying states in nuclei, for instance, a giant resonance, are spread over a certain energy interval due to the coupling to a background which is complicated and has huge degrees of freedom. As a result, the strength function, at a large energy scale, exhibits a global shape profile like a Lorentzian. The values of the center of energy or the total width depend on the properties of a probe such as multiporality, isospin, and so on, and change smoothly as the mass number of a nucleus. Most of these values are well understood up to now . On the other hand, at small energy scale limit comparable to the level spacing of the background states, the strength function may rapidly fluctuate from state to state on top of the global profile, if we neglect the escaping width due to the coupling to the continuum. It is believed that the fluctuation properties at this energy scale can be well simulated by the random matrix theory of the Gaussian orthogonal ensemble (GOE) , which has been verified by several experimental data: The nearest-neighbor level spacing distribution (NND) of neutron resonances near the threshold approximately follows the Wigner distribution . The distribution of the reduced width (proportional to an absolute square of a component of the wave function) of these resonance also shows the Porter-Thomas distribution . Furthermore, the NND of the Nuclear Data Ensemble shows the Wigner distribution and the $`\mathrm{\Delta }_3`$ statistics of them shows $`\mathrm{\Delta }_3(L)\mathrm{ln}L`$ . All of these are typical signatures of the GOE.

Thus, we know well the behavior of the strength function at both extreme small and large energy scales. From this knowledge, one may construct the following model to describe the strength function of highly excited states in nuclei: A collective state couples democratically to each of unperturbed background states, and the background Hamiltonian itself is the GOE random matrix. This GOE background model is closely connected with the pandemonium picture . If we adjust the energy of the collective state and the coupling strength between the collective state and background states, the resultant strength function should be consistent with the giant resonance strength function at least at both extreme energy scale limits. Therefore, in order to find the difference between a realistic strength function and that from this model, we must investigate fluctuation properties of the strength function at intermediate energy scale .

Recently, we proposed a new method to perform fluctuation analysis of the strength function . This method is devised to quantitatively characterize a fluctuation property as a function of energy scale by a new measure called the local scaling dimension (LSD), and is suitable to investigate fluctuation at intermediate energy scales. We applied this method to a simple doorway damping model. The doorway damping model is different from the above GOE background model with respect to two points: First, unperturbed background states are divided into two different classes so that unperturbed states in one class, which are called the doorway states, couple directly to a collective state, while states in another class do not. The second point is that doorway states and other unperturbed background states have a finite spreading width due to the mixing among them. (Note that the GOE background model corresponds to the large spreading width limit.) Fluctuation properties of the strength functions of both models are same at small energy scale limit. The result of the analysis clearly showed that fluctuation properties of both strength functions deviate from each other at a certain intermediate energy scale, and this energy scale is closely related to the spreading width of the doorway states in the doorway damping model.

In this paper, we apply the same method to a more realistic strength function of the giant resonances in a nucleus. Strength functions of the isoscalar (IS) and the isovector (IV) giant quadrupole resonance (GQR) in 40Ca are obtained by the shell model calculation within up to the 2p2h configurations. The strength functions reproduce the experimental center of energy and the total width, and also show the GOE fluctuation at the small energy scale limit, as will be shown later. Furthermore, in this model, the doorway structure and the deviation from GOE in the background may be introduced in a natural way, while they are input by hand in the doorway damping model. Therefore these strength functions may provide the better test for the analysis.

Note that a similar method of analysis was proposed and applied to the ($`e,e^{}`$) or ($`p,p^{}`$) experimental data on 208Pb .

The paper is organized as follows: In Sec. II, we briefly explain the LSD. In Sec. III, we apply the method to the strength functions of the isoscalar (IS) and the isovector (IV) GQR. We also discuss the physical origin of the deviation from the GOE background model in detail. Finally, Sec. IV is devoted to conclusions.

## II Local Scaling Dimension

We briefly explain the local scaling dimension in this section. See Ref.  for detail.

The strength function is expressed as

$$S(E)=\underset{i}{}S_i\delta (EE_i+E_0),$$
(1)
for exciting the nucleus with excitation energy $`E`$. Here $`E_i`$ and $`E_0`$ are the energy of discrete levels and the ground state energy, respectively, and $`S_i`$ denotes the strength of exciting the $`i`$th energy level. Let us assume strengths are normalized as $`_iS_i=1`$.

In order to quantitatively characterize how the strength function fluctuate at various energy scale, we consider binned distribution of the strength function $`S(E)`$ by dividing whole energy interval under consideration into $`L`$ bins with length $`ϵ`$. Strength contained in $`n`$th bin is denoted by $`p_n`$,

$$p_n\underset{in\mathrm{th}\mathrm{bin}}{}S_i.$$
(2)
To characterize the distribution of the binned strengths, we introduce the moments of it, which are called in literature the partition function $`\chi _m(ϵ)`$ defined by

$$\chi _m(ϵ)\underset{n=1}{\overset{L}{}}p_n^m=Lp_n^m.$$
(3)
Clearly, the partition functions contains, at small $`ϵ`$ limit, the information of the strength distribution , as well as the energy-strength correlation as a function of a bin width or an energy scale $`ϵ`$.

It is then helpful to employ the idea of the generalized fractal dimension , which is derived from the partition function and has a definite physical meaning. The concept of the generalized fractal dimension is useful only when a system has a multifractal nature, namely, has a self-similar structure against the change of the scale. Strength functions at a nucleus generally do not have such a property. Thus, we have to extend the concept of the generalized fractal dimension and finally reach the local scaling dimension (LSD) defined by

$$D_m(ϵ)\frac{1}{m1}\frac{\mathrm{log}\chi _m(ϵ)}{\mathrm{log}ϵ}.$$
(4)
Different from the generalized fractal dimension, the local scaling dimension characterizes “how the partition function scales for different energy scales” and accordingly is a function of a bin width or an energy scale $`ϵ`$. The local scaling dimension, however, has an easily understandable physical meaning similar to those of the generalized fractal dimension: For instance, the value of $`D_m(ϵ)`$ close to unity means that the strength function looks like a one dimensional object distributed smoothly over an energy interval when we look at the strength function at an energy scale $`ϵ`$. On the other hand, it looks like a zero dimensional object, namely a dot if the value of $`D_m(ϵ)`$ is close to zero. Note that the local scaling dimension reduces to the generalized fractal dimension when the value of $`D_m(ϵ)`$ is constant over a whole interval of $`ϵ`$.

In actual calculation of the local scaling dimension, we define it by means of the finite difference under the change of a factor 2,

$$D_m(\sqrt{2}ϵ)=\frac{1}{m1}\frac{\mathrm{log}\chi _m(2ϵ)\mathrm{log}\chi _m(ϵ)}{\mathrm{log}2},$$
(5)
rather than the derivative in Eq. (4). Using the finite difference, the calculation is very simple for all the moments.

One of the important strength functions may be the one obtained from the GOE random matrix. The local scaling dimension of this GOE strength function can be used as a reference to study other strength functions which are supposed to have a similar fluctuation properties to the GOE, as in the case of nuclear strength functions in the highly excited region. We show in Fig. 1 an example of the GOE strength function and the local scaling dimension $`D_m(l)`$, where $`lϵ/d`$ (the bin width measured in the unit of the level spacing $`d`$) is used as the scaling parameter.

## III Analysis of the giant quadrupole resonances in 40Ca

We apply the local scaling dimension to the strength functions of the giant quadrupole resonances in 40Ca.

### A Numerical Calculation of Strength Functions

We calculated the strength functions of the isoscalar and the isovector giant quadrupole resonances in 40Ca within the second Tamm-Dancoff approximation including the 1p1h and 2p2h excitations. Single-particle wave-functions and energies were obtained in terms of a Woods-Saxon potential including the Coulomb interaction. To simulate the bare (Hartree-Fock) single-particle energies $`\epsilon _{HF}`$, the above single-particle energies $`\epsilon _{WS}`$ were scaled by the effective mass $`m^{}/m`$ as $`\epsilon _{HF}=\epsilon _{WS}/(m^{}/m)`$. As the residual interaction, the Landau-Migdal-type interaction  including the density-dependence was adopted. The model space was constructed in terms of single-particle states within the four major shells, two below and two above the Fermi surface, and included all 1p1h states and 2p2h states whose unperturbed energies were less than 50MeV. Resultant number of 1p1h states and 2p2h states are 34 and 4144, respectively. We diagonalized the Hamiltonian within this model space and obtained the strength functions for the isoscalar and the isovector quadrupole operators.

Figure 2 (a) and (b) show the strength functions for ISGQR and IVGQR, respectively. The mean energy and the total width (the standard deviation) are 23MeV(31MeV) and 4.0MeV(7.0MeV) for the ISGQR(IVGQR). The effective mass $`m^{}/m=0.9`$ was used to reproduce the empirical mean energy for IVGQR in 40Ca . These strength functions were already used in Refs.  and  to study the generic properties of nuclear giant-resonance decay.

### B Several Measures

Before going to the detailed discussion for the local scaling dimension, we briefly show other fluctuation measures.

Figure 3 shows the nearest-neighbor level spacing distribution (NND), where levels between 20MeV and 30MeV are considered (the number of levels is 804), and are unfolded by means of the Strutinsky method with a smoothing width 5MeV . The NND follows rather well the Wigner distribution, which indicates at the very small energy scale the level fluctuation almost obeys the GOE. We also plot in Fig. 4 the $`\mathrm{\Delta }_3`$ statistics. We again find that at small energy range the $`\mathrm{\Delta }_3`$ follows the GOE line. At $`L^{\mathrm{max}}15`$, however, $`\mathrm{\Delta }_3`$ starts to deviate from the GOE line to upward. This value of $`L^{\mathrm{max}}`$ roughly corresponds to the energy scale $`E0.18`$MeV.

Figure 5 shows the strength distribution where the distribution of the square-root of normalized strengths $`\overline{S}_i`$ is plotted as a histogram; strengths are normalized also by Strutinsky method, detailed procedure of which will be explained in the next subsection. We can see that the strength distribution follows considerably well the Porter-Thomas distribution.

### C Results of the Local Scaling Dimension

Since we are not interested in the global shape of the strength function, it is convenient to introduce the normalized strength function where the global smooth energy dependence is removed from the original strength function as follows:

$$\overline{S}_i=𝒩\frac{S_i\stackrel{~}{\rho }(E_i)}{\stackrel{~}{S}(E_i)}.$$
(6)
Here, $`\stackrel{~}{S}(E)`$ and $`\stackrel{~}{\rho }(E)`$ are obtained by averaging the calculated strength function, Eq.(1), and the level density $`\rho (E)=_i\delta (EE_i+E_0)`$, respectively, by the Strutinsky method. We adopted 5MeV as a smoothing width. $`𝒩`$ is an overall normalization factor to guarantee $`_i\overline{S}_i=1`$ for the considered energy interval. Moreover, we adopt as energy levels the equidistant ones , namely, $`E_i=id`$, where $`d`$ denotes average level spacing. By adopting these levels, we can neglect global energy dependence of the level density and the local level fluctuation. (We verified that the inclusion of the local level fluctuation affects only the behavior of the LSD at small energy scales, and this effect can be described by the GOE.) Thus, we obtain the normalized strength function $`\overline{S}(E)`$;

$$\overline{S}(E)=\mathrm{\Sigma }_i\overline{S}_i\delta (EE_i+E_0).$$
(7)
The local scaling dimension is derived from this normalized strength function. Finally, since the visible strengths of ISGQR are almost concentrated in the energy range between 20MeV and 30MeV, we adopt this energy interval for analysis.

The normalized strength functions for ISGQR and IVGQR are plotted in Figs. 6 (a) and (b), respectively.

First, let us discuss the ISGQR case in detail. The case of the IVGQR will be discussed in Sec. III D 4. Figures 7(a) and (b) show the partition functions $`\chi _m(ϵ)`$ defined by Eq. (3) and the local scaling dimensions $`D_m(ϵ)`$ defined by Eq. (5), respectively. The horizontal axes in both figures represent the bin width of energy axis in unit of $`d`$, where $`d`$ represents the average level spacing over the energy range 20MeV and 30MeV ($`d12`$keV). From the partition function $`\chi _m(ϵ)`$, one sees that the fluctuation does not follow the power scaling law (or linear relation in the log-log plot). Correspondingly, the local scaling dimension $`D_m(ϵ)`$ varies as a function of $`ϵ`$. At small energy scale, the values of $`D_m(ϵ)`$ gradually increase as the energy scale $`ϵ`$ increases, and almost follow the GOE line. As the energy scale $`ϵ`$ increases further, the values of $`D_m(ϵ)`$ bend downwards and then again turn to increase. The behavior of $`D_m(ϵ)`$ around this energy scale clearly deviates from the GOE limit. The value of the energy scale where the deviation is maximum is roughly $`ϵ_\mathrm{f}140d1.7`$MeV. From these observations, we again find that at small energy scale, fluctuation of the strength function is essentially governed by the GOE, as one knows from other measures such as the NND or the strength distribution. The scaling analysis by means of the local scaling dimension, however, can reveal the new feature which may not be realized by other measures. Namely, the GOE fluctuation is limited to the small energy scale, while, at intermediate energy scales, a different fluctuation from the GOE arises.

One may wonder that the statistical error could produce a fictitious dip, since the number of levels considered 804 is not so large. Of course, there may be a statistical error. Accordingly, we should take the energy scale $`ϵ_\mathrm{f}140d1.7`$MeV mentioned above as a rough estimate. The existence of the dip itself, however, is not due to a statistical error. This would be verified from a systematic analysis in the next subsection. Note that the simple estimate of the standard deviation, of the $`m`$th local scaling dimension is $`\sigma (D_m)\sqrt{3}m/\sqrt{N_{\mathrm{tot}}}`$ , and in the present case for the second local scaling dimension, $`\sigma (D_2)0.12`$.

A type of the smoothing procedure to obtain the normalized strength function does not much affect the results. Indeed, we verified that the Gaussian smoothing instead of the Strutinsky smoothing leads to the essentially same results.

### D Origin of the Deviation from the GOE

We shall investigate the origin of the observed energy scale $`ϵ_\mathrm{f}=1.7`$MeV characteristic to the fine structures in the isoscalar quadrupole strength function.

#### 1 Doorway damping mechanism

For this purpose we consider a damping mechanism by separating the shell model space to 1p1h and 2p2h subspaces, and we focus on the residual interactions within and between the subspaces. The giant resonance is a collective vibrational state whose dominant component is a coherent superposition of 1p1h configurations. The collective state can be described by the Tamm-Dancoff (TD) approximation, which corresponds to a truncated shell model where only 1p1h configurations are taken into account. The result of the TD approximation (TDA) is shown in Fig. 8(a), and is compared with the full shell model calculation including up to 2p2h’s (Fig. 2(a)). It is seen that the peak position is well reproduced by the TDA. Note also that the isoscalar quadrupole strength in the TDA is spread over many eigenstates, which are correlated 1p1h states. This spreading of the quadrupole strength among 1p1h states is often called the Landau damping. Comparing Figs. 8(a) and 2(a), the Landau damping is found to give dominant contribution to the total FWHM of the strength function.

If we take into account the 2p2h configurations, the interaction acting between 1p1h and 2p2h states (abbreviated as $`V_{12}`$ ) takes part in. Since the 2p2h states have much more degrees of freedom, $`V_{12}`$ causes spreading of the TD 1p1h states. Because of this spreading of 1p1h states, the quadrupole strength is also fragmented by $`V_{12}`$ on top of the Landau damping at the TDA level. To see the effect of $`V_{12}`$, we performed a calculation where unperturbed 2p2h states and the interaction $`V_{12}`$ is introduced. The resultant strength distribution, shown in Fig. 8(b), appears similar to the final result (Fig. 2(a)) as far as the gross profile of the strength distribution is concerned.

Note also that there exists the residual interaction acting among 2p2h states, which we call $`V_{22}`$ in the following. $`V_{22}`$ causes configuration mixing among the 2p2h states. Comparison between Figs. 8(b) and 2(a) indicates effects of $`V_{22}`$ on fine details of strength distribution, although the gross profile of strength distribution is not influenced very much.

The above observations lead us to the following doorway damping picture. The quadrupole strength distribution is fragmented first in the 1p1h subspace. The 1p1h TD states then spread further over 2p2h states through $`V_{12}`$. Here the 1p1h TD states may be considered as doorway states in the whole damping processes. The 2p2h states on the other hand may be regarded as states that play roles of background in the main damping processes since they do not dominate the total damping width, but they influence the fine structures through the coupling to the doorway states.

Assuming the above picture, we can consider energy scales that may be relevant to the fine structures of the quadrupole strength distribution. Concerning the 1p1h TD states which are considered as the doorway states, we have as energy scales 1) the spreading width (which we denote $`\gamma _{12}`$) of the TD 1p1h states caused by the interaction $`V_{12}`$, in addition to 2) the average level spacing $`d_{1p1h}`$ between the 1p1h states. Concerning the 2p2h states, 3) the level spacing $`d_{2p2h}`$ between 2p2h states is to be noted. ($`d_{2p2h}`$ is almost identical to the level spacing $`d`$ for the whole set of energy spectra since the number of 2p2h states is much larger than that of 1p1h states. ) Note also that mixing among 2p2h states caused by the interaction $`V_{22}`$ is characterized by 3) the spreading width ($`\gamma _{22}`$ named in the following) of 2p2h states. We have found previously that, in the case of a schematic model which incorporate doorway states coupled to background states, the spreading width of doorway states (corresponding to $`\gamma _{12}`$ in the present context) determines the energy scale where the local scaling dimension deviates from the generic GOE limit . Note also that the deviation from the GOE is seen in the energy level fluctuation by the analysis of the $`\mathrm{\Delta }_3`$ statistics. The energy level fluctuation is sensitive to the residual interaction that causes mixing among unperturbed configurations. We found that the energy level fluctuation is governed by $`V_{22}`$, but not by $`V_{12}`$, as shown in Fig. 9. Thus we can consider another possibility that the energy scale of the deviation from the GOE limit may be related to the spreading width $`\gamma _{22}`$ of the 2p2h states. In the following, we shall focus on the effects of $`V_{12}`$ and $`V_{22}`$ which are connected to the two spreading widths $`\gamma _{12}`$ and $`\gamma _{22}`$.

#### 2 Spreading widths of 1p1h and 2p2h states

Let us quantify $`\gamma _{12}`$ and $`\gamma _{22}`$. If we assume the Fermi golden rule, we obtain a simple estimate of $`\gamma _{12}`$ by $`\gamma _{12}^{\mathrm{FG}}=2\pi \overline{1p1h|V_{12}|2p2h}^2/d_{2p2h}`$ where $`\overline{1p1h|V_{12}|2p2h}`$ is a root mean square value of the matrix elements between 1p1h and 2p2h states. Similarly, we can make an estimation of the spreading width $`\gamma _{22}`$ as $`\gamma _{22}^{\mathrm{FG}}=2\pi \overline{2p2h|V_{22}|2p2h}^2/d_{2p2h}.`$ The above estimate gives $`\gamma _{12}^{\mathrm{FG}}5.3`$MeV, and $`\gamma _{22}^{\mathrm{FG}}6.0`$MeV. Here we have inserted $`\overline{1p1h|V_{12}|2p2h}=98`$ keV, $`\overline{2p2h|V_{22}|2p2h}=104`$ keV, $`d_{2p2h}=11.3`$keV, which are evaluated for all the unperturbed 1p1h and 2p2h states in the energy interval of $`E=2030`$MeV. Note also $`d_{1p1h}0.5`$MeV for the 1p1h states in the same interval. Very similar values are obtained for $`\overline{1p1h|V_{12}|2p2h}`$ when we use the correlated 1p1h states (TD states) instead of unperturbed 1p1h’s.

It is also possible to evaluate $`\gamma _{12}`$ and $`\gamma _{22}`$ in a more direct way by using the strength functions of 1p1h and 2p2h states. Strength function of a 1p1h state $`|i,1p1h`$ is given by

$$S_{i,1p1h}(e)=\underset{a}{}i,1p1h|a^2\delta (e(E_aE_{i,1p1h})),$$
(8)
where $`|a`$ and $`E_a`$ are the eigenstates of the full calculation and their energies, respectively. $`e`$ denotes the relative energy from the energy centroid $`E_{i,1p1h}`$ of the 1p1h state considered. Evaluating an average of $`S_{i,1p1h}(e)`$ over all the 1p1h states, the average strength function of 1p1h states is obtained, as shown in Fig. 10(a). Its FWHM gives an evaluation of the average spreading width $`\gamma _{12}`$ of 1p1h states. Similar method can be used to evaluate the spreading width $`\gamma _{22}`$ of the 2p2h states (cf. Fig. 10(b)). This direct method gives $`\gamma _{12}^{\mathrm{dm}}=1.5`$MeV and $`\gamma _{22}^{\mathrm{dm}}=5.2`$MeV.

It is noted that the values of the spreading width $`\gamma _{12}`$ are very different between the golden rule and the direct method. To investigate the difference, we evaluated root mean square values of interaction matrix elements of $`V_{12}`$ for individual 1p1h states, and found that about one third of the 1p1h states have significantly large value $`120200`$ keV as compared with the other 1p1h states which have small value $`3070`$keV. Corresponding golden rule estimates of the spreading width are $`\gamma _{12}^{\mathrm{FG}}=8.322`$MeV and $`\gamma _{12}^{\mathrm{FG}}=0.52.8`$MeV, respectively. The FWHM of the average strength function of 1p1h states used in the direct method is dominated by the 1p1h states having small FWHM. Accordingly the direct method gives much smaller value ($`\gamma _{12}^{\mathrm{dm}}=1.5`$ MeV) than that of the golden rule ($`\gamma _{12}^{\mathrm{FG}}=5.3`$ MeV).

#### 3 Dependence of LSD on $`V_{12}`$ and $`V_{22}`$

In order to identify the origin of the observed energy scale $`ϵ_\mathrm{f}=1.7`$MeV, we shall analyze the dependence of the fine structures of the strength functions on $`V_{12}`$ and $`V_{22}`$, which are intimately related to the spreading widths $`\gamma _{12}`$ and $`\gamma _{22}`$. For this analysis we rescale the interaction strength $`V_{12}`$ and $`V_{22}`$ by factors $`\lambda _{12}`$ and $`\lambda _{22}`$ respectively, and look into the local scaling dimension to see how the fine structures change.

It is seen from Fig. 11 that the fine structures probed by the local scaling dimension $`D_m(ϵ)`$ strongly depend on $`\lambda _{12}`$. When the value of $`\lambda _{12}`$ is small, $`\lambda _{12}=0.1`$, the LSD exhibits very strong deviation from the GOE limit in almost all energy scales except at very small energy scale $`ϵ/d1`$. As the value of $`\lambda _{12}`$ increases, the deviation becomes smaller and the LSD at small energy scales coincides with the GOE limit. The energy scale $`ϵ^{}`$ where the LSD start to deviate from the GOE limit moves to larger values as $`\lambda _{12}`$ increases. If $`\lambda _{12}`$ is increased more than the original value (e.g. $`\lambda _{12}=1.7`$), the LSD almost reaches the GOE limit over almost whole energy scales. The above dependence on $`\lambda _{12}`$ agrees very well with a similar behavior found in the simple doorway damping model . (In that schematic model one can account for the spreading width $`\gamma _{\mathrm{dw}}`$ of the doorway states analytically.) We found in the schematic model that the energy scale $`ϵ^{}`$ is proportional to $`\gamma _{\mathrm{dw}}`$ as $`ϵ^{}\gamma _{\mathrm{dw}}/5`$ and the energy scale $`ϵ^{}`$ where the LSD decreases most steeply corresponds to $`\gamma _{\mathrm{dw}}`$. $`ϵ^{}`$ and $`ϵ^{}`$ moves with the interaction $`v`$ between the doorway states and the background states since $`\gamma _{\mathrm{dw}}v^2`$. This behavior is consistent with the $`\lambda _{12}`$ dependence in the present model. For quantitative evaluation, we compare the spreading width $`\gamma _{12}`$ of the 1p1h TD states, which is considered here as the doorway states, and the energy scale $`ϵ_\mathrm{f}`$ where the deviation of the LSD from the GOE limit is maximum. Here we adopt the energy scale $`ϵ_\mathrm{f}`$ in place of $`ϵ^{}`$, whose extraction is not easy in the present model due to fluctuation in the LSD. The energy scale $`ϵ_\mathrm{f}=1.7`$MeV found for the original interaction (i.e. $`\lambda _{12}=1`$) agrees with the spreading width $`\gamma _{12}`$ of the 1p1h states if we adopt the value $`\gamma _{12}^{\mathrm{dm}}=1.5`$MeV obtained by the direct evaluation. The agreement is a little worse if we refer to the value $`\gamma _{12}^{\mathrm{FG}}=5.3`$MeV derived from the golden rule estimate. We expect that the fine structures are dominated by the group of 1p1h states having smaller spreading width $`\gamma _{12}1`$MeV, while the remaining 1p1h states with $`\gamma _{12}>10`$MeV, which is larger than the total width of the giant resonance, do not contribute to the fine structure. The comparison with the spreading width $`\gamma _{12}^{\mathrm{dm}}`$ by the direct evaluation would be more appropriate since $`\gamma _{12}^{\mathrm{dm}}`$ reflects mainly the first group. (On the other hand, $`\gamma _{12}^{\mathrm{FG}}`$ based on the root mean square of interaction matrix elements puts emphasis on the other group. See the previous subsection.) The $`\lambda _{12}`$-dependence and the comparison with $`\gamma _{12}`$ indicate an approximate relation $`ϵ_\mathrm{f}\gamma _{12}`$, for which we estimate an ambiguity by a factor of two.

Fig. 12 shows $`\lambda _{22}`$-dependence of the LSD. We immediately see that dependence on $`\lambda _{22}`$ is much weaker than that on $`\lambda _{12}`$. For $`\lambda _{22}<0.2`$, the LSD exhibits deviation from the GOE limit for rather wide interval of energy scale, but the deviation itself is not very large. As $`\lambda _{22}`$ increases, the LSD approaches the GOE limit at small energy scales. (At $`\lambda _{22}=0.5`$, for instance, the local scaling dimensions almost follow the GOE when $`ϵ10d`$.) For $`\lambda _{22}>0.5`$ the LSD is not sensitive to the change in $`\lambda _{22}`$, and there remains the same deviation from the GOE limit even up to the maximal choice of $`\lambda _{22}=2.0`$.

The insensitivity to $`\lambda _{22}`$ for $`\lambda _{22}>0.5`$ may be understood as follows: We first note that the energy scale of the mixing among 2p2h states is given by the spreading width $`\gamma _{22}`$ of the 2p2h states, and that the mixing properties are expected to show up in the GOE behaviors for small energy scales satisfying $`ϵ<\gamma _{22}`$. In the case of the original strength $`\lambda _{12}=\lambda _{22}=1`$, the spreading width $`\gamma _{22}^{\mathrm{dm}}`$ =5.2MeV of the 2p2h states is larger than the spreading width $`\gamma _{12}^{\mathrm{dm}}=1.5`$MeV of the 1p1h TD states. If $`\lambda _{22}<0.5`$, opposite relation $`\gamma _{22}<\gamma _{12}`$ is realized (n.b. $`\gamma _{22}\lambda _{22}^2`$ and $`\gamma _{12}\lambda _{12}^2`$), and hence the energy scale where the deviation from the GOE limit starts is determined by $`\gamma _{22}`$. In the case of $`\lambda _{22}>0.5`$ (including the case of the original strength), on the other hand, we have $`\gamma _{12}<\gamma _{22}`$, and the deviation from the GOE limit is governed by $`\gamma _{12}`$. In principle, we can consider a possibility that the LSD exhibits the two energy scales $`\gamma _{12}`$ and $`\gamma _{22}`$ separately. However the deviation related to $`\gamma _{22}`$ is small even when $`\gamma _{22}<\gamma _{12}`$ ($`\lambda _{22}<0.5`$), making it difficult to extract the energy scale associated with $`\gamma _{22}`$. In the case of the original interaction strength, the spreading widths $`\gamma _{22}^{\mathrm{dm}}=5.2`$MeV is close to the total width $`\gamma _{GR}`$ ($``$ 5MeV) of the giant resonance, and thus it is difficult to be probed. (Notice that we adopted a smoothing width 5MeV.)

Summarizing the above two analyses, the characteristic energy scale of the fine structures is related to the spreading width $`\gamma _{12}`$ of the 1p1h TD states, which is considered as the doorway states in the giant resonance decay.

#### 4 Isovector strength function

We have also analyzed the strength distribution for the isovector quadrupole operator. Figures 13(a) and (b) show the partition function and the local scaling dimension, respectively, for the isovector strength function in the energy interval $`E=`$20-30MeV. The behavior of the LSD is similar to that of isoscalar strength function. Namely, the local scaling dimension almost follows the GOE limit at smaller energy scales whereas the deviation from the GOE limit is recognized at intermediate energy scales. We also verified that the change of the value of $`ϵ_\mathrm{f}`$ is quite sensitive to the change of $`\lambda _{12}`$ compared to the change of $`\lambda _{22}`$ as in the case of the ISGQR. Characteristic energy scale of the deviation reads $`ϵ_\mathrm{f}70d0.9`$MeV. This value is smaller than that for the isoscalar strengths ($`ϵ_\mathrm{f}1.7`$MeV) by about a factor of two. Both are consistent with the spreading width $`\gamma _{12}^{\mathrm{dm}}=1.5`$MeV.

## IV Conclusions

We have analyzed fine structures of the giant quadrupole resonances described by means of the microscopic shell model including up to 2p2h configurations. By applying the scaling analysis of strength fluctuations that utilizes the local scaling dimension, we extracted the energy scale $`ϵ_\mathrm{f}=1.7`$ MeV for the ISGQR and $`ϵ_\mathrm{f}=0.9`$ MeV for the IVGQR that characterize the fine structure. We discussed the origin of the characteristic energy scale in terms of the spreading widths of the 1p1h and 2p2h states. We found a clear correlation between the characteristic energy scale and the spreading width $`\gamma _{12}`$ of the 1p1h Tamm-Dancoff states, which play a role of doorway states in the present shell model description of the giant quadrupole resonances in 40Ca.

In the present paper, we intended to illustrate how the scaling analysis can be utilized to investigate the fine structures and their origin. Although the adopted model has some realistic features, we do not claim that the present shell model predicts all the quantitative features of the quadrupole strength distributions since the particle escaping width for instance is neglected in the present model. Also, other mechanisms such as the surface phonon coupling are not explicitly taken into account in the present model. Apparently, one may need more realistic model if one intends to make a quantitative prediction or comparison with the experiment. Instead we would like to emphasize that the local scaling analysis provides a general tool for quantitative investigation of the fine structures, and the scaling analysis itself can be applied to any kinds of models which exhibit fine structures in the strength functions.


text:
# Split-off dimer defects on the Si(001)2 × 1 surface

## I Introduction

There are currently several exciting proposals to use the (001) surface of silicon for the construction of atomic-scale electronic devices, including single electron transistors Tucker and Shen (2000), ultra-dense memories Quaade et al. (2001) and quantum computers Kane (1998); O’Brien et al. (2001). However, since any random charge or spin defects in the vicinity of these devices could potentially destroy their operation, a thorough understanding of the nature of crystalline defects on this surface is essential. The Si(001) surface was first observed in real space at atomic resolution using scanning tunneling microscopy (STM) by Tromp et. al.Tromp et al. (1985) in 1985. In this study they observed the surface consisted of rows of “bean-shaped” protrusions which were interpreted as tunneling from the $`\pi `$-bonds of surface Si dimers, thereby establishing the dimer model as the correct model for this surface. Since then, STM has been instrumental in further elucidating the characteristics of this surface, and in particular atomic-scale defects present on the surfaceHamers and Köhler (1989); Zhang et al. (1996); Hata et al. (2000a).

The simplest defect of the Si(001) surface is the single dimer vacancy defect (1-DV), shown schematically in Figs. 1(a) and 1(b). This defect consists of the absence of a single dimer from the surface and can either expose four second-layer atoms (Fig. 1(a)) or form a more stable structure where rebonding of the second-layer atoms occurs Wang et al. (1993) as shown in Fig. 1(b). While the rebonded 1-DV strains the bonds of its neighboring dimers it also results in a lowering of the number of surface dangling bonds and has been found to be more stable than the nonbonded structure. Owen et al. (1995); Wang et al. (1993) Single dimer vacancy defects can also cluster to form larger defects such as the double dimer vacancy defect (2-DV) and the triple dimer vacancy defect (3-DV). More complex clusters also form, the most commonly observedKoo et al. (1995); Wang et al. (1993) example is the 1+2-DV consisting of a 1-DV and a 2-DV separated by a single surface dimer, the so-called “split-off dimer”. The accepted structure of the 1+2-DV, as proposed by Wang et. al. based on total energy calculations,Wang et al. (1993) is shown in Fig. 1(c) and consists of a rebonded 1-DV (left), a split-off dimer, and a 2-DV with a rebonding atom (right). Recently we have observed another DV complex that contains a split-off dimer, called the 1+1-DV, which consists of a rebonded 1-DV and a nonbonded 1-DV separated by a split-off dimer, as shown in Fig. 1(d).

Here we present a detailed investigation of DV defect complexes that contain split-off dimers. Using high-resolution, low-bias STM we observe that split-off dimers appear as well-resolved pairs of protrusions under imaging conditions where normal Si dimers appear as single “bean-shaped” protrusions. We show that this difference arises from an absence of the expected $`\pi `$-bonding between the two atoms of the split-off dimer but instead the formation of $`\pi `$-bonds between the split-off dimer atoms and second layer atoms. Electron charge density plots obtained using first principles calculations support this interpretation. We observe an intensity enhancement surrounding some split-off dimer defect complexes in our STM images and thereby discuss the local strain induced in the formation of these defects. Finally, we present a model for a previously unreported triangular-shaped split-off dimer defect complex that exists at SB-type step edges.

## II High-resolution variable-bias STM imaging of defect complexes

Experiments were performed in two separate but identical variable temperature STM systems (Omicron VT-STM). The base pressure of the ultra-high vacuum (UHV) chamber was $`<5\times 10^{11}`$ mbar. Phosphorus doped $`10^{15}`$ and $`10^{19}`$ cm-3 wafers, orientated towards the  direction were used. These wafers were cleaved into $`2\times 10`$ mm2 sized samples, mounted in sample holders, and then transferred into the UHV chamber. Wafers and samples were handled using ceramic tweezers and mounted in tantalum/molybdenum/ceramic sample holders to avoid contamination from metals such as Ni and W. Sample preparationSwartzentruber et al. (1989) was performed in vacuum without prior ex-situ treatment by outgassing overnight at 850 K using a resistive heater element, followed by flashing to 1400 K by passing a direct current through the sample. After flashing, the samples were cooled slowly ($`3`$ K/s) from 1150 K to room temperature.

### II.1 Split-off dimers

The sample preparation procedure outlined above routinely produced samples with very low surface defect densities. However, the density of defects, including split-off dimer defects, was found to increase over time with repeated sample preparation and STM imaging, as reported previously.Hata et al. (2000b) It is known that split-off dimer defects are induced on the Si(001) surface by the presence of metal contamination such as Ni, Zandvliet et al. (1995) and W Matsuura et al. (2000). The appearance of these defects in our samples therefore points to a build up of metal contamination, either Ni from in-vacuum stainless steel parts, or more likely W contamination from the STM tip. After using an old W STM tip to scratch a $``$ 1 mm line on a Si(001) sample in vacuum and then reflashing, the concentration of split-off dimer defects on the surface was found to have dramatically increased, confirming the STM tip as the source of the metal contamination.

Figure 2 shows an STM image of a Si(001) surface containing a $``$ 10% coverage of split-off dimer defects. The majority of the defects in this image can be identified as 1+2-DVs, however, two 1+1-DVs are also present, as indicated. The most striking feature of this image is the difference in appearance of the split-off dimers in contrast to the surrounding normal surface dimers. Each split-off dimer in this image appears as a double-lobed protrusion, while the surrounding normal Si dimers each appear as a single “bean-shaped” protrusion, as expected at this tunneling bias. Hata et al. (1999) Line profiles taken across a 1+2-DV both parallel and perpendicular to the dimer row direction are shown in Fig. 2(b). The line profile parallel to the dimer row direction agrees with previously reported profiles over 1+2-DVs and fits well with the accepted structure, Zhang et al. (1996); Owen et al. (1995) as shown by the overlayed ball and stick model. The line profile taken perpendicular to the dimer row direction, however, clearly shows that the split-off dimer of this defect is separated into two protrusions while the neighboring Si dimers are single protrusions. This is the first recognition and explanation of split-off dimers appearing as double-lobed protrusions.

To understand why split-off dimers appear as double-lobed protrusions we must consider the structure of these defects shown in Figs. 1(c) and 1(d). Normally Si(001) surface dimers appear as “bean-shaped” protrusions in STM images because the dangling bonds of each Si dimer atom mix to form a $`\pi `$-bond between the two dimer atoms. However, if we examine the split-off dimer structure closely (Figs. 1(c) and 1(d)) we see that unlike normal surface dimers, the split-off dimer has two nearest neighbor second layer atoms that each have a dangling bond. The separation distance between the split-off dimer atoms and these second layer atoms is sufficiently close to allow the formation of $`\pi `$-bonds. The resulting four-atom structure can therefore be referred to as a tetramer. We propose that the four dangling bonds of the split-off dimer tetramer interact primarily along the backbonds between the split-off dimer atoms and the second layer atoms to form $`\pi `$-bonds down the backbonds, as drawn schematically in Fig. 2(c). These two spatially separated $`\pi `$-bonds therefore lead to the double-lobed appearance of the split-off dimers under low bias filled-state tunneling conditions, which we confirm in section III with charge density calculations.

In an attempt to fully characterize the appearance of these split-off dimers in STM images, we have performed a series of experiments observing split-off dimers with changing STM sample bias. Figure 3 summarizes our results, showing images where a 1+2-DV and a 1+1-DV located next to each other are observed at four different sample biases – two filled-state images and two empty-state images. In the filled-state image of Fig. 3(a) we see that at $`0.8`$ V the split-off dimer of both the 1+2-DV and the 1+1-DV appear as double-lobed protrusions similar to those in Fig. 2(a). However, when the filled-state bias is increased in magnitude to $`2`$ V, Fig. 3(b), the split-off dimers become single-protrusions and appear very similar to the surrounding normal Si surface dimers. This is because as the bias magnitude is increased towards $`2`$ V, the dimer $`\sigma `$-bond and bulk states contribute increasingly to the tunneling current Hata et al. (1999) and the image of the split-off dimer reverts to the bean-shaped protrusion in the same manner as normal surface Si dimers. In both of the empty-state images, Figs. 3(c) and 3(d), acquired at +0.8 V and +2 V, respectively, the appearance of the split-off dimers is very similar to that of the surrounding normal surface dimers. This is because under empty-state tunneling conditions electrons tunnel into the $`\pi ^{}`$-antibonding orbitals of the dimers, resulting in the normal Si dimers appearing as double-lobed protrusions. Qin and Lagally (1999) It is therefore only under low bias magnitude filled-state tunneling conditions that split-off dimers appear significantly different to the surrounding normal Si surface dimers.

### II.2 Experimental observation of surface strain in complex defects

Another noticeable feature of Figs. 2(a) and 3(a) is the enhanced brightness of the 1+1-DV compared to the 1+2-DV. This is a reproducible effect that we attribute to an increased amount of surface strain induced by the 1+1-DV. Figure 4 shows a series of adjacent defects forming a short vacancy line channel in the surface. This channel is composed of individual 1-DV, 3-DV, 1+2-DV, and 1+1-DV defects (see figure caption). In the filled-state image, Fig. 4(a), there is a clear brightening of the dimers on one end of the 1+1-DVs and the dimers on both ends of the 1-DV, which is not present for the 1+2-DVs. In the empty-state image of the line of defect complexes, Fig. 4(b), we notice that there is a darkening of the same dimers that are enhanced in the filled-state image.

Owen et. al.Owen et al. (1995) have shown using low bias STM and first principles calculations, that the dimers neighboring a rebonded 1-DV are enhanced in low bias filled-state STM images due to the strain induced by the defect shifting the surface states upwards in energy toward the Fermi energy. This effect can be seen for the 1-DV in Fig. 4(a), where the neighboring dimers in the same row as the 1-DV are enhanced in intensity, with the magnitude of the enhancement decaying with distance from the 1-DV. A very similar enhancement can be seen around the 1+1-DV sites in this image, with the split-off dimer in particular appearing much brighter than the surrounding normal surface dimers. However, for the 1+1-DV only the dimers on one end of the defect are enhanced in intensity while the dimers on the other end of the defect are not. This observation can be readily explained since the 1+1-DV is composed of a rebonded 1-DV adjacent to a nonbonded 1-DV (Fig. 1(d)) and Owen et. al. Owen et al. (1995) have shown that while the rebonded 1-DV results in strain-induced image enhancement, the nonbonded 1-DV does not. The observation of an asymmetric strain-induced enhancement of the 1+1-DV in Fig. 4(a) can therefore be taken as an experimental confirmation of the structure of this defect (Fig. 1(d)) and the first application of the method of Owen et. al. Owen et al. (1995) for identifying strain in more complex surface defect structures.

The fact that the 1+2-DV causes no enhancement of its neighboring dimers over the surrounding normal surface dimers suggests that the 1+2-DV, unlike the 1-DV and 1+1-DVs, does not increase the strain of the surface. This at first seems strange, since the 1+2-DV involves a rebonded 1-DV similar to the 1+1-DV structure. However, Wang et. al. Wang et al. (1993) have shown, using total energy calculations, that the junction formed between the 1-DV and the 2-DV to create the 1+2-DV releases the surface strain that is present when these two defects exist separately from one another. The STM data that we have presented here is therefore the first experimental verification of this calculation. The fact that both the 1-DV and the 1+1-DV show local enhancement due to strain, while the 1+2-DV does not, indicates that the 1+2-DV structure induces less local strain than the 1-DV.

In their paper, Owen et. al. do not present empty state STM images, nor do they consider empty states in their tight binding calculations. In Fig. 4(b), we show an empty state image of the same line of defects shown in Fig. 4(a). Interestingly, in this empty state image the dimers that were enhanced in brightness surrounding the 1-DV and 1+1-DVs in the filled-state image are less bright than the surrounding Si dimers in the empty-state image. This suggests that the strain associated with these defects causes the lowest unoccupied molecular orbital (LUMO) of the adjacent dimers to also shift higher in energy, away from the Fermi energy.

## III Density functional characterization of defect complexes

To confirm the interpretation of our STM images, we have performed first principles electronic structure calculations of both the 1+2-DV and 1+1-DV complexes using the Car-Parrinello Molecular Dynamics program. Hutter et al. (1995-2002) Valence electrons were described using Goedecker pseudopotentials Goedecker et al. (1996) expanded in a basis set of plane waves with an energy cutoff of 18 Rydbergs and the exchange-correlation functional was of the BLYP form. Becke (1988); Lee et al. (1988) Slab calculations contained between 124 and 128 Si atoms in a $`31.070\times 7.675\times 19.253`$ Å3 supercell, corresponding to six layers of vacuum in the $`z`$-direction, and all calculations were performed with gamma point sampling of the Brillouin zone only. A reference calculation was performed with no surface vacancies and assuming the $`p(2\times 2)`$ structure in which the dimers buckled alternately along the row. A single 256 atom calculation with a duplication along the y-axis confirmed that the effect of dispersion across the rows is minor as has been noted elsewhere. Pollmann et al. (1987)

Both zero temperature geometry optimization and high temperature molecular dynamics calculations were used to explore a variety of surface and second-layer bonding configurations for the 1+2-DV and 1+1-DV. The results confirm the configurations in Figs. 1(c) and 1(d) are the lowest energy geometries of both defect complexes. The dimers are drawn symmetric in these schematics, however, the true minimum energy structure at zero temperature involved charge-transfer buckling of the Si dimers. It is well known that at room temperature the barrier is sufficiently small for the dimers to flip-flop between the two equivalent configurations. Ramstad et al. (1995); Wolkow (1992) Our calculations show that the split-off dimer tetramer also has two symmetrically equivalent buckling configurations, with charge transfer between the atoms of the tetramer buckling adjacent atoms in alternate directions. By analogy with the normal dimers we can expect room-temperature STM measurements of the tetramer to image the average of the two configurations. The chemical potential was determined from a 512-atom bulk calculation, which yielded a formation energy of 0.85 eV for the 1+2-DV, similar to the value of 0.65 eV computed by Wang at. alWang et al. (1993) The 1+1-DV formation energy has not been previously reported, and we found it to be 1.13 eV. We note that this value is high, but this is consistent with the rarity of observation of the 1+1-DV in STM experiments.

In Fig. 5 we present a series of calculated electron density slices through various regions of the 1+2-DV marked by (a), (b), (c), and (d) in the ball and stick schematic. The charge density shown in the figure is the sum of the occupied Kohn-Sham orbitals within 0.25 eV of the highest occupied molecular orbital (HOMO). Taking into account the $`0.5`$ eV surface band gap of Si(001) and the n-type doping of the experimental samples, these states correspond approximately to the accessible states for a $`0.75`$ V sample bias and can therefore be directly compared to the experimental data in Fig. 3(a), which was acquired with a $`0.8`$ V sample bias.

The four charge density slices in Fig. 5 show: Fig. 5(a) the 1-DV edge dimer, Fig. 5(b) the split-off dimer, Fig. 5(c) the backbond of the split-off dimer, and Fig. 5(d) the 2-DV edge dimer, as indicated schematically in Fig. 5(e). The charge densities of both buckling configurations of the dimers and backbond atoms are averaged, and the positions of the dimer and tetramer atoms are shown superimposed in both buckling configurations. In the case of the backbonds, the two configurations are not coincident, and so the atoms and bonds are shown in projection onto the plane in Fig. 5(c). The 1-DV edge dimer in Fig. 5(a) shows a clear three-lobed character with significant overlap between the up-atom charge density of the two buckling orientations, and a single lobe beneath the plane of the surface at the mid-point of the dimer. Density functional calculations by Hata et. al. Hata et al. (1999) and tight-binding Green’s function calculations by Pollman et. al. Pollmann et al. (1987) have separately identified this three-lobed feature as being characteristic of $`\pi `$-bonding in flip-flop dimers on the silicon surface, and we can therefore take this three-lobed feature as a signature of $`\pi `$-bonding in this work. The backbond of the split-off dimer in Fig. 5(c) connects a first-layer atom to a second-layer atom and also shows a three-lobed structure. By analogy with the surface dimer in Fig. 5(a) we characterize this bond as having $`\pi `$-character and have indicated this by the shaded ellipse (c) shown in Fig. 5(e). The split-off dimer itself in Fig. 5(b), however, does not exhibit three-lobed character. Instead, the split-off dimer has four lobes; two located above the up-atoms of the dimer in each buckling configuration, and a second pair of spatially separated lobes beneath the bond. The calculations thus show that $`\pi `$-bonding occurs down the backbonds of the split-off dimer, but not across the dimer itself. The absence of the $`\pi `$-bond across the split-off dimer correlates with the double-protrusions observed in the STM images. Finally, we also consider the charge density of the 2-DV edge dimer, Fig. 5(d), and note that it also exhibits three-lobed character, indicative of $`\pi `$-bonding. This gives the 2-DV edge dimer a bean-shaped appearance in the STM image, as for the 1-DV dimer in Fig. 5(a).

A similar situation exists for the 1+1-DV charge density slices shown in Fig. 6. The first three charge density slices, Figs. 6(a) – 6(c), are analogous to the slices for the 1+2-DV As was the case for the 1+2-DV, the rebonded 1-DV edge dimer, Fig. 6(a) and the split-off dimer backbonds, Fig. 6(c) exhibit three-lobed $`\pi `$-like character, while the split-off dimer, Fig. 6(b) exhibits four-lobed character, consistent with an end-on view of $`\pi `$-bonding down the backbonds. Finally, another slice is presented in Fig. 6(d), which is through the nonbonded 1-DV edge dimer as indicated schematically in Fig. 6(e). It can be seen that the nonbonded 1-DV edge dimer appears quite different to the charge density slices discussed so far. In particular, we notice that the nonbonded 1-DV edge dimer has a much reduced charge density compared to the other slices, Fig. 6(a) – 6(c). Examination of the structure identifies strain as the characteristic that differentiates the dimer in Fig. 6(d) from the other dimers. Since the dimer in Fig. 6(d) is part of a tetramer, one might expect its appearance to resemble the split-off dimer which is also part of the tetramer shown Figs. 6(b) and 6(c). However, a detailed examination of the simulated structure reveals that the nonbonded 1-DV tetramer is relaxed, since there is one adjacent dimer present, while the split-off tetramer is highly strained because of the rebonding in the second-layer. Since the nonbonded 1-DV tetramer is much less strained, its occupied states lie further from the Fermi level, explaining the charge reduction observed in calculations in Fig. 6(d). As discussed in Ref. Da̧browski and Müssig, 2000, the minimum energy arrangement of the electrons in a tetramer is one where the $`\pi `$-states are delocalized across the four atoms, to form three bonding segments, as indicated by the ellipses in Fig. 6(e). The charge density slice of Fig. 6(d) is consistent with such an arrangement where the charge density is shared between $`\pi `$-like bonds on both backbonds and across the dimer atoms. We conclude that this charge density arrangement forms for the nonbonded 1-DV tetramer because it is allowed to relax. In the case of the split-off dimer, the tetramer is constrained by the rebonding and instead forms a higher energy configuration in which the $`\pi `$-bonds conjugate to form two $`\pi `$-bonds down its backbonds.

## IV New step edge defect

Having presented a detailed understanding of the electronic structure of previously observed split-off dimer defects in the Si(001) surface using both STM and first-principles calculations, we now turn our attention to elucidating the structure of a previously unreported split-off dimer defect. In Figs. 7(a) and 7(b) we show filled- and empty-state STM images of DV defects at a single-layer SB-type step edge. At the top of these images white arrows indicate are three defects known as SB-DVs, which are rebonded 1-DVs at the step edge, which leave a single split-off dimer as the last dimer before the lower terrace begins. Koo et al. (1996) As was the case for the 1+1-DV and 1+2-DV, the split-off dimers in SB-DVs appear as double-lobed protrusions under low-bias filled-state imaging conditions, Fig. 7(a). At the bottom of Fig. 7(a) two similar DV complexes can be observed, as indicated by black arrows, however these defects have a third protrusion giving them a triangular appearance. In empty-state imaging, Fig. 7(b), however, the additional third feature is not present. These triangular-shaped defects have not been reported on the Si(001) surface before and most likely arise due to the presence of W contamination.

Our proposed structural model of the triangular-shaped defects in Fig. 7(a) is shown in Fig. 7(c). This model consists of a nonbonded 1-DV defect at an SB-type step edge, followed by a rebonded split-off dimer and a bound Si monomer. Swartzentruber has previously observed Si monomers on the Si(001) surface using high-resolution STM after depositing a few percent of a monolayer of Si atoms to the surface. Swartzentruber (1998) These monomers were bound at rebonded SB-type step edges, confirming the minimum energy binding position predicted by first principles calculations. The binding position of the monomer in our proposed structure, Fig. 7(c), is essentially the same position observed by Swartzentruber, with the difference being the presence of the DV defect adjacent to the step edge. Swartzentruber also observed that the Si monomers bound at SB-type step edges were visible in one bias polarity (empty-state) but invisible in the other (filled-state). Our images reveal a similar effect, however the feature we observe appears in filled-state images while being invisible in empty-state images.

We have performed first-principles calculations to produce charge density contours for our proposed structure. Figure 7(d) shows a constant $`z`$-height contour slice taken 1.2 Å above the monomer for occupied Kohn-Sham orbitals within 0.45 eV of the HOMO. We see in this charge density contour slice the two lobes expected for the split-off dimer as well as a third lobe due to the bound monomer. Moreover, the distance between the split-off dimer lobes and the monomer lobe is 6.4 Å in agreement with the separation seen in the STM image. In Fig. 7(e) we show an empty-state slice taken at the same $`z`$-height and summed over Kohn-Sham orbitals up to 0.45 eV above the LUMO. In this contour the double lobe of the split-off dimer is still present but the monomer lobe is significantly lessened in intensity. The results of our first-principles calculations therefore give good agreement between our proposed structure and the observed defect. The presence of the split-off dimer must therefore be responsible for the reversal of the filled- and empty-state monomer characteristics when compared to those observed for monomers bound to rebonded SB-type step edges.

## V Summary

We have investigated split-off dimers on the Si(001)2$`\times `$1 surface using high resolution STM and first principles calculations. We find that split-off dimers form $`\pi `$-bonds with second layer atoms which gives them a double-lobed appearance in low bias filled-state STM images. We apply the method of Owen et. al. Owen et al. (1995) for identifying local areas of increased surface strain to dimer vacancy defect complexes and thereby present the first experimental confirmation of the predicted strain relief offered by the 1+2-DV. Finally, we have presented a previously unreported triangular-shaped defect on the Si(001) surface and a proposed model for this structure involving a bound Si monomer.