text:
# Comparative study of dimer-vacancies and dimer-vacancy lines on Si(001) and Ge(001)

## 1 Introduction

Epitaxial growth of germanium on the Si(001) surface has been an active area of investigation because it represents a versatile way of fabricating quantum dots. At the same time, Ge/Si(001) provides an example of a system where the mismatch strain can be used to control the morphological features of the surface. In the early stages of growth, the surface changes its reconstruction from the $`2\times 1`$ structure to a $`2\times N`$ pattern, which consists of an array of long, parallel lines of dimer vacancies. It has recently been shown that the $`2\times N`$ reconstruction can be an excellent template for the nanofabrication of 1-dimensional structures with unique electronic transport behavior . At later stages of growth, the film breaks up into pyramidal islands (quantum dots); controlling the size, spatial distribution and electronic properties of these quantum dots represent key issues for the practical applications of the heteroepitaxial growth on Si(001). From a fundamental point of view, the current state-of-the-art in scanning tunnelling microscopy has helped renew the general interest in the $`2\times N`$ structures on Si(001). For example, Rastelli and coworkers  have recently discovered that the $`2\times N`$ reconstruction on the Ge-covered Si(001) surface disappears upon Si capping, confirming in an ingenious way the role played by the mismatch strain in determining the stability of the vacancy lines. In another recent study , Sutter et al. have used high-resolution STM imaging in conjunction with a detailed statistical analysis of morphological features of the surface to establish that the interactions between steps and vacancy lines cause the initial surface roughening of the Ge/Si(001).

Rapid progress has been made after the discovery of the $`2\times N`$ structure , and two authoritative reviews summarizing the present understanding of the growth of germanium films on Si(001) have been published . Dimer vacancies (DVs) form easily on the Si(001) surface, but they only organize into vacancy lines (VLs) when there is a sufficient amount of Germanium deposited (close to 1 ML). The mismatch strain created by epitaxial deposition increases the number of dimer vacancies per unit area, and the interactions between these vacancies become important at sufficiently high vacancy concentrations. The interaction between two DVs is long-range repulsive when the DV’s belong to the same dimer row, but weakly attractive when they belong to adjacent dimer rows . Such anisotropy of the interactions determines the alignment of the DVs in the direction perpendicular to the dimer rows. This physical picture for the origin of the $`2\times N`$ reconstruction has emerged from the experimental work of Chen et al. . The mechanism for the formation of vacancy lines is consistent not only with experiments done on other epitaxial systems (e.g., Bi/Si(001) ), but also with the observations of the Si(001) surface with metal impurities such as Ag  or Ni .

While the physical origin of the vacancy lines on Si(001) has been elucidated , it remains a puzzle as to why vacancy lines are not normally observed on the very similar Ge(001) surface. The fact that the Ge(001) surface is virtually defect free  suggests that the formation energy of the vacancies on Ge(001) is substantially larger than that of the vacancies on Si(001). Quantitative data from atomic-scale calculations on the Ge system is not available, although the need for such calculations was emphasized almost a decade ago in the work of Yang et al. . The purpose of our paper is to address the differences between Ge(001) and Si(001) concerning the formation of dimer-vacancies and their self-assembly into vacancy lines. The work presented here fits well into the current efforts  for explaining the differences between the two technologically important surfaces.

We have calculated the formation energies of point and line defects on the (001) surfaces and their dependence on applied strain, using empirical potential  models for interactions . The results obtained offer an explanation for the differences between the vacancies on Si(001) and Ge(001) and are consistent with experimental findings; moreover, these results could foster further work based on more accurate descriptions of the atomic interactions such as tight-binding or density functionals. We have considered both the heteroepitaxial system Ge/Si(001) for different values of Ge coverage, as well as pure Si and Ge surfaces under external strain. While the results concerning the epitaxial Ge/Si(001) systems have been previously obtained in different forms , the variation of the energetic parameters of DVs and VLs on the (001) surface with applied strains has not been systematically addressed in the past. Such a study of the energetics of defects on (001) surfaces as a function of applied strain is presented here. The new results obtained in this paper are:
(a) the formation energy of the dimer-vacancies on Ge(001) and its decrease under compressive, biaxial strains. In the case of unstrained surfaces, we have found that the formation energy for single DV on Ge(001) is three times larger than the corresponding value for the single DV on Si(001). Therefore, in order for the formation energy of a single-DV on Ge(001) to become negative, a large compressive strain should be applied, either biaxially or in the direction of the dimer rows.
(b) the formation energy of vacancy lines on Si(001) and Ge(001) as a function of applied strain. Using a simple model based on the interactions of nearest-neighbor vacancies in a VL, we have also identified the formation energy of a single DV and the attractive interaction (binding) between two DVs in adjacent dimer-rows as the dominant contributions to the formation energy of a VL.

The organization of this article is as follows. In section 2 we outline the competition between the formation energy of the vacancy lines and their elastic interactions as the general energetic mechanism that leads to the appearance of the $`2\times N`$ reconstructions. In contrast to previous interpretations , we have found that there is no long range attraction between the vacancy lines, as their interactions are purely repulsive. As in most other studies to date, the formation of vacancy lines is regarded here as a purely thermodynamic phenomenon, thus the kinetics concerning the diffusion of vacancies during their alignment are not considered. In section 3 we describe the computational details of the total energy calculations based on Tersoff potentials . The next two sections are devoted to the formation energies and repulsion strengths for single-vacancy lines (VLs) on Ge-covered Si(001) (section 4) and for different types of VLs on pure Si(001) and Ge(001) (section 5). We investigate the formation energies of individual dimer vacancies on (001) surfaces and their strain dependence, identifying the differences between the defects on Si and Ge surfaces. The possible relevance of these results for future experimental work is also discussed in section 5. Our conclusions are summarized in section 6.

## 2 Energy of vacancy lines on $`2\times N`$ surfaces

The $`2\times N`$ reconstruction is a structural pattern on the (001) surface of Si, characterized by the periodic lengths $`2a`$ and $`Na`$, where $`N`$ is a positive integer and $`a=3.84`$Å is the lattice constant of the unreconstructed Si(001) surface. This reconstruction is obtained by eliminating every $`N`$-th dimer from each dimer row, such that the dimer vacancies created in this manner are in registry with one another and form straight lines perpendicular to the dimer rows. To lower the number of dangling bonds at the surface, the second-layer atoms below the dimer vacancy rebond such that they are fully coordinated, at the cost of introducing some surface stress . Depending on the preparation details of the surface, other types of atomic structures of VLs can be obtained . In this section we focus on the general features of the energetics of the $`2\times N`$ surface, without reference to a specific atomic structure. We start by partitioning the surface energy into a contribution from the $`2\times 1`$ surface and a contribution from the vacancy lines:

$$\gamma _{2\times N}=\gamma _{2\times 1}+\frac{\lambda }{L_x},$$
(1)
where $`\gamma _{2\times N}`$ and $`\gamma _{2\times 1}`$ denote the surface energies of $`2\times N`$ and $`2\times 1`$, respectively, $`\lambda `$ is the energy per unit length of the vacancy lines, and $`1/L_x`$ is the density of lines in a periodic array of VLs with a spacing of $`L_xNa`$ . Eq. (1) defines the energy per unit length $`\lambda `$ of the vacancy lines as an excess energy with respect to the defect-free $`2\times 1`$ surface. We note that the energy of the vacancy lines $`\lambda `$ includes the formation energy of the lines, as well as their elastic interactions. Guided by the early work of Marchenko and Parshin , we expect the interactions to be inversely proportional with the square of the separation between VLs, so that the line energy $`\lambda `$ can be written as

$$\lambda =\mathrm{\Lambda }+\frac{\pi ^2}{6}\frac{G}{L_x^2},$$
(2)
where $`\mathrm{\Lambda }`$ is the formation energy of a VL and $`G`$ is the strength of the repulsion between lines. The numerical factor $`\pi ^2/6`$ arises since we treat a periodic array of line defects , rather than two isolated VLs. Using Eq. (2) in (1) we can write the surface energy of the $`2\times N`$ reconstruction as a function of $`N`$:

$$\gamma _{2\times N}=\gamma _{2\times 1}+\frac{\mathrm{\Lambda }}{Na}+\frac{\pi ^2}{6}\frac{G}{N^3a^3}$$
(3)
where we have used the fact that $`L_x=Na`$.

Analyzing Eq. (3), we note that:
(a) if $`\mathrm{\Lambda }`$ is positive then $`\gamma _{2\times N}>\gamma _{2\times 1}`$, so the $`2\times 1`$ structure is energetically favorable.
(b) for any $`\mathrm{\Lambda }<0`$, Eq. (3) has only one minimum located at:

$$N=N^{}\frac{\pi }{a}\sqrt{\frac{G}{(2\mathrm{\Lambda })}}$$
(4)
For the case of negative formation energy $`\mathrm{\Lambda }`$, a typical plot of the surface energy difference $`(\gamma _{2\times N}\gamma _{2\times 1})`$ as a function of $`N`$ is shown in Fig. 1.

The condition $`\mathrm{\Lambda }<0`$ is sufficient to ensure the stability of the $`2\times N`$ reconstruction (over the $`2\times 1`$ structure), because for $`N=N^{}`$ we have $`\gamma _{2\times N}\gamma _{2\times 1}=(2/3)\mathrm{\Lambda }<0`$. This analysis of Eq. (3) shows that, irrespective of the atomic-scale features of the VLs, the $`2\times N`$ reconstruction is always stable when the formation energy of a vacancy line is negative. The optimal value of $`N`$ is given by the competition between the formation energy $`\mathrm{\Lambda }`$ and the strength $`G`$ of the repulsive interactions, as shown in Eq. (4). The magnitudes of $`\mathrm{\Lambda }`$ and $`G`$ are determined by atomic-scale details, and incorporate the effects of the VL structure , wetting layer composition , as well as the anisotropy of surface stress and relaxation .

The following two sections are devoted to the computation of the line formation energies and repulsion strengths for various VL structures, and to the investigation of the behavior of these quantities as a function of germanium coverage and applied strain. Before we proceed with the numerical results, we draw a connection with previous work on the energetics of the $`2\times N`$ surface. The above theory of the energetics of the vacancy lines is similar to Refs. , since it explicitly considers the competition between the (negative) formation energy of a VL and the repulsive (positive) interaction between the VLs as the physical origin of the $`2\times N`$ structure. In Ref.  however, the emphasis falls on the surface stress anisotropy, and the treatment of the two-domain (i.e. $`2\times N`$ and $`N\times 2`$) surfaces obliterates the simplicity of the VL energetics outlined above.

A word of caution is due concerning the mechanisms of formation of the $`2\times N`$ reconstruction proposed in Ref.  and in the review article . In , the conclusions based on the energetic competition between ”trench formation” and ”strain relaxation” are fortuitous, since the $`N`$-dependence of each of these two energetic contributions is incorrect 222We believe that the energetic competition that determines the $`2\times N`$ reconstruction is described by Eq. (3).. The sign of the ”trench energy” and the magnitudes of the two energetic contributions are also incorrect, though these aspects are probably outside the scope of the simple Frenkel-Kontorova model used in . Ref.  asserts that ”in the dimer-row direction the VL-VL interaction is short-range repulsive and long-range attractive”, which seems plausible given the behavior of the $`\gamma _{2\times N}`$ in the limit of large $`N`$. However, we note that what appears to be an attractive long-range interaction between VLs is in fact the formation energy of one VL per unit cell of the $`2\times N`$ surface. As shown in the sections to follow, the VL-VL interaction takes a purely repulsive form once the vacancy line energy is extracted from the surface energy.

## 3 Computational Details

### 3.1 Structure of the vacancy lines

In the terminology of Ref. , the vacancy lines considered here are made of 1-DV, 2-DV and 1+2-DVs, as shown in Fig. 2. The choice of these structures is motivated by the experimental observations of the atomic-scale details of vacancy lines on different surfaces: the 1-DV structure is present on Ge-covered Si(001) , while lines with the 2-DV and 1+2-DV configurations have been observed on Si(001) with Ni  and Ag  impurities, as well as on radiation-quenched Si(001) .

### 3.2 Total energy calculations

Accurate calculations of the long-range behavior of the VL-VL interactions are very demanding, requiring slabs that are both long ($`N>10`$) and thick (at least 10 layers). Earlier density-functional and tight-binding calculations were restricted to thin slabs with 4–6 moving layers, due to limitations in computational power . At present, density-functional and tight-binding calculations have become possible , though they remain computationally expensive for the purpose of determining the long-range interactions between line defects. Such calculations are still prohibitive for point defects on the surface, since these defects require simulation cells that are large in all three directions. As more studies would probably become available for various types of vacancy lines in the near future, it is worth pointing out that a good test for their total-energy convergence is to check that the line energies $`\lambda `$ can be fitted well with the inverse-square power law given by Eq. (2). Such a test would also indicate whether or not the chosen thickness of the computational slab is sufficient. This type of analysis is routinely performed for steps on crystal surfaces (see  and references therein) in order to determine the step formation energy and step-step interactions, but it has almost never been applied to the case of vacancy lines.

Since we will investigate vacancy lines as well as individual vacancies (point defects) on the surface, we employ empirical potentials to model the atomic interactions. We use the Tersoff potential for Si and Ge , since it satisfactorily describes the structure and energetics of steps on (001) surfaces , as well as the surface energies of the $`2\times N`$ reconstructions . Despite known limitations, the empirical potential described in  represents a widely accepted standard for Si-Ge systems with large numbers of atoms. Our aim is restricted to a meaningful comparison (rather than an accurate determination) of the properties of dimer vacancies on Si(001) and Ge(001). For this purpose, the use of the Tersoff potential is appropriate, and even desirable as a starting point.

For simulation cells with a thickness of $`200`$Å and values of $`N`$ in the range $`4N36`$, we have performed conjugate-gradient structural optimizations. From the minimized total energy $`E`$ of the atoms in the simulation cell we compute the surface energy $`\gamma _{2\times N}`$. When there is only one atomic species in the simulation cell (e.g., Si), the surface energy is given by:

$$E=ne_{Si}+\gamma _{2\times N}L_xL_y,$$
(5)
where $`e_{Si}`$ is the bulk cohesion energy per Si atom and $`n`$ denotes the number of atoms in the simulation cell that are allowed to relax. For the case of two atomic species in the simulation cell, the surface energy is given by:

$$E=n_{Ge}\mu _{Ge}+n_{Si}\mu _{Si}+\gamma _{2\times N}L_xL_y,$$
(6)
where $`n_{Ge}`$, $`\mu _{Ge}`$ ($`n_{Si}`$, $`\mu _{Si}`$) are the number of atoms and chemical potential of Ge (Si). Since we are considering deposition on Si(001), the chemical potential of Si is set to the bulk cohesion energy per atom $`e_{Si}`$. For low Ge coverages (1–3ML), the chemical potential of Ge ($`\mu _{Ge}`$) depends on the deposition conditions and is not known, although its value is usually set to the bulk cohesion energy of (unstrained) Ge, $`e_{Ge}`$ . Exceptions are found in the work of Oviedo et al.  and Li et al. , who carefully define a reference state from which the surface energies are measured; the procedure in  is equivalent to computing the chemical potential of Ge based on a wetting layer composition assumed to be thermodynamically stable. From the results of Refs.  we infer that for a pure Ge wetting layer (no intermixing), setting $`\mu _{Ge}=e_{Ge}`$ is a poor approximation for monolayer and sub-monolayer coverages, but it becomes reasonable for 2 ML and 3 ML of Ge coverage. Though we use $`\mu _{Ge}=e_{Ge}`$ in the present work, we will also address the effects of changing the chemical potential of Ge.

## 4 Vacancy lines on the Ge-covered Si(001) surface

In this section, we calculate the surface energy of the $`2\times N`$ surface with 0–3 ML Ge coverage, then compute the line formation energy and the VL-VL interactions. Our results for the surface energies are very close to the ones originally obtained by Liu and Lagally ; the deviations are due to the slightly different versions of the Tersoff potentials 333The parameter $`\lambda _3`$ of the Tersoff potential given in  is set to zero in , thus modifying the original T3 parametrization for pure Si. Here, we have found it convenient to use a form of the interactions that reduces to  in the limit where only one atomic species is present in the simulation cell.. What differentiates the results of this section from those of Ref.  is that we pursue the determination of the line formation energy and VL-VL interaction strength, whereas Liu and Lagally address the interplay between surface stress and stoichiometry .

Using Eq. (1), we have computed the line energy $`\lambda `$ and plotted it as a function of the inverse square of the line spacing ($`1/L_x^2`$) in Fig. 3. The line energy $`\lambda `$ is a linear function of $`1/L_x^2`$ for each integer value of the Ge coverage between 0ML and 3ML (Fig. 3). This behavior confirms the expectation (Eq. (2) ) that the interactions between the vacancy lines are dipolar in nature . The parameters of the linear fits shown in Fig. 3 give the line formation energy $`\mathrm{\Lambda }`$ (intercept) and the line-line repulsion strength $`G`$ (slope), and are listed in Table 1. The table also includes the fitting parameters for two cases of intermixed wetting layer compositions. We have verified that the the dipolar character of the line-line interactions is not changed by intermixing in the wetting layer. These intermixing cases are considered here only for the purpose of illustration, as the actual structure of the wetting layer at given Ge coverage is under debate .

The data in Table 1 shows that for both of the wetting layer compositions considered (pure Ge or intermixed), the formation energy decreases strongly with increasing Ge coverage, which is consistent with previous findings (refer to ). The repulsion strength $`G`$ is very sensitive to the WL composition, as the elastic interactions are largely transmitted through the first few layers: different intermixed compositions affect the elastic properties of the film and thus the strength of the repulsive interactions.

The trends shown in Table 1 are well understood  and physically meaningful, although we emphasize that the exact values of the formation energy $`\mathrm{\Lambda }`$ are influenced by the choice of the chemical potential for the deposited Ge atoms. By taking finite differences in Eq. (6), we note that a change $`\delta \mu _{Ge}`$ in the chemical potential of Ge introduces a change $`\delta \gamma _{2\times N}`$ in the surface energy of the $`2\times N`$ structure given by:

$$0=n_{Ge}(\delta \mu _{Ge})+(2Na^2)(\delta \gamma _{2\times N}),$$
(7)
where we have used the values for the periodic lengths, $`L_x=Na`$ and $`L_y=2a`$. Similarly, for a $`2\times 1`$ surface with the same area, thickness and WL composition, we have

$$0=(2+n_{Ge})(\delta \mu _{Ge})+(2Na^2)(\delta \gamma _{2\times 1}),$$
(8)
where the number of Ge atoms was increased by 2, as no dimers are missing from the $`2\times 1`$ structure. Subtracting Eq. (8) from Eq. (7) and using Eq. (2) we obtain

$$\delta \mu _{Ge}=Na^2(\frac{\delta \mathrm{\Lambda }}{Na}+\frac{\pi ^2}{6}\frac{\delta G}{N^3a^3}),$$
(9)
where $`\delta \mathrm{\Lambda }`$ and $`\delta G`$ are the changes in the line formation energy and repulsion strength generated by the change in chemical potential $`\delta \mu _{Ge}`$. From Eq. (9) in the limit of large $`N`$, we find that $`\delta \mathrm{\Lambda }=\delta \mu _{Ge}/a`$ and $`\delta G0`$. Therefore, modifications of the Ge chemical potential amount to shifting the line formation energy while leaving the VL-VL interactions unchanged –a conclusion which is independent of the model used to describe the atomic interactions.

In summary, we have computed the line formation energies and the repulsion strengths for certain Ge coverages and WL compositions. As shown in Table 1, the line formation energy $`\mathrm{\Lambda }`$ is sensitive to the WL thickness and composition, while the VL-VL interactions are mostly sensitive to composition. To gain further insight into how strain affects the interactions and self-organization of vacancy lines, we will now consider $`2\times N`$-reconstructed surfaces of pure species (either Si or Ge) under external applied strain.

## 5 Defects on Si(001) and Ge(001) under applied strain

The study of vacancy lines and dimer vacancies on uniformly strained Si(001) and Ge(001) surfaces has two technical advantages: (a) the ambiguity regarding the value of chemical potential of Ge is removed, and (b) the uncertainties related to intermixing in the wetting layer obviously do not arise when studying systems with only one atomic species. Therefore, the study of defects on clean surfaces under external strain allows for a more direct comparison between Si(001) and Ge(001).

### 5.1 Vacancy lines on strained surfaces

The surface energy of the $`2\times N`$ structure under external strain can be computed using formula (5) in which the bulk cohesion energy is also a function of the applied strain. The in-plane deformations are applied along the dimer rows ($`ϵ_x\epsilon `$, $`ϵ_y=0`$), perpendicular to the dimer rows ($`ϵ_x=0,ϵ_y\epsilon `$), or equibiaxially ($`ϵ_x=ϵ_y\epsilon `$). For each type of deformation we vary the strain in the range $`2\%ϵ2\%`$, by appropriately scaling the atomic coordinates. In this range of deformation, we find that the line formation energy $`\mathrm{\Lambda }`$ is a linear function of strain:

$$\mathrm{\Lambda }(ϵ_x,ϵ_y)=\mathrm{\Lambda }_0+\mathrm{\Lambda }_\alpha ^{}\epsilon ,$$
(10)
where $`\mathrm{\Lambda }_0`$ is the line formation energy in the absence of strain, and the subscript $`\alpha `$ denotes the uniaxial strain along different directions ($`\alpha =x`$, or $`\alpha =y`$) or the equibiaxial strain ($`\alpha =xy`$). A consequence of the linear relation described in Eq. (10) is that the biaxial strain sensitivity $`\mathrm{\Lambda }_{xy}^{}`$ is the sum of the two sensitivities $`\mathrm{\Lambda }_x^{}`$ and $`\mathrm{\Lambda }_y^{}`$ that correspond to the uniaxial deformations:

$$\mathrm{\Lambda }_{xy}^{}=\mathrm{\Lambda }_x^{}+\mathrm{\Lambda }_y^{}.$$
(11)
For the VL structures shown in Fig. 2, we have computed all the strain sensitivities $`\mathrm{\Lambda }_\alpha ^{}`$ ($`\alpha =x,y,xy`$) independently, and found that formula (11) holds to within 1–2% of $`\mathrm{\Lambda }_{xy}^{}`$.

These strain sensitivities are listed in Table 2, along with the formation energies $`\mathrm{\Lambda }_0`$ and repulsion strengths $`G_0`$ computed in the absence of strain. The strain-dependence of the line formation energy $`\mathrm{\Lambda }`$ is strongly anisotropic. As shown in Table 2, the strain sensitivity $`\mathrm{\Lambda }_x^{}`$ along the dimer row is positive, while $`\mathrm{\Lambda }_y^{}`$ is negative and several times smaller (up to one order of magnitude for 1-DV vacancy lines) in magnitude than $`\mathrm{\Lambda }_x^{}`$. This indicates that for all the vacancy structures shown in Fig. 2, the rebonding in the second layer plays a crucial role in determining the strain dependence of the line formation energy. In the absence of strain, the bonds that ”bridge” the second-layer atoms in the vacancy lines are stretched compared to their bulk values, and are especially sensitive to strains applied parallel to them. Therefore, compressive strains along the dimer rows will decrease the amount of stretch in those bonds, lowering the formation energy of the vacancy lines. This correlation between atomic bonding and formation energy is consistent with the ones given previously for the case of 1-DV lines (see, for example ). On the other hand, an uniaxial strain applied perpendicular to the dimer rows only indirectly affects the formation energy of the vacancy lines via a Poisson-type effect: a stretch along the $`y`$-direction creates a small compressive stress in the $`x`$-direction, which in turn lowers the line formation energy.

While the formation energy $`\mathrm{\Lambda }`$ of the vacancy lines varies significantly with applied strain, the VL-VL repulsion strength $`G`$ shows a weak strain dependence, changing by less than ten percent for the values of strain considered here. The reason for the weak strain-dependence of the VL interactions is that the relaxations that occur after applying external strains to an already relaxed $`2\times N`$ surface are rather small, since no bonds are broken or formed. Consequently, the VL–VL interactions are not substantially affected by external strains as long as these strains remain in the small-deformation regime given by Eq. (10).

By comparing the Si and Ge data (Table 2) for each type of VL structure, we find several systematic trends: (a) the repulsion strengths $`G`$ are stronger for Si than for Ge, (b) the formation energies $`\mathrm{\Lambda }_0`$ are smaller for Si than for Ge and (c) the strain sensitivities along dimer row ($`\mathrm{\Lambda }_x`$) and biaxial ($`\mathrm{\Lambda }_{xy}`$) are larger for Si than for Ge. Most notably, the formation energy of a single-DV line on Si(001) is negative, while for the case of Ge(001) it is positive. While there is a possibility that this sign difference between Si(001) and Ge(001) might be due to artifacts of the empirical potential, it is consistent with the observation of self assembly of DVs on Si(001)  and not on Ge(001) . Some of these trends are related to the formation and binding of individual dimer vacancies, as will be shown in the next section.

The occurrence of the vacancy lines is a consequence of the strongly anisotropic interactions between the individual dimer vacancies, coupled with the preference of the DVs to diffuse along the dimer rows. As discovered by Chen et al. , two DVs that belong to the same dimer row repel one another, while DVs in adjacent dimer rows experience a weak attraction at short distances. This finding was confirmed by other independent experiments , as well as by theoretical studies . Thus the physics of vacancy lines is intrinsically related to the formation, interactions and diffusion of individual dimer-vacancies. In what follows, we show how the strain dependence of the formation energy of the VLs is linked to the behavior of individual dimer vacancies under strain.

### 5.2 Dimer vacancies on Si(001) and Ge(001)

In this section, we compute the formation energy of three types of DVs on Si(001) and Ge(001), as well as the attractive energy of two vacancies placed in adjacent dimer rows. In the case of Si(001), a similar calculation was performed by Wang and coworkers , and has served as useful reference for many subsequent investigations. In contrast, we have found no similar studies for the formation energies of defects on Ge(001). Here we perform such a study using the Tersoff potential , and identify the differences between the vacancies on the two surfaces, Si(001) and Ge(001). We denote the formation energy of a DV complex by $`u`$, and that of an isolated pair of complexes lying in adjacent dimer rows (called a bound pair form here on) by $`u_p`$. The binding energy of the pair is defined 444With this definition, the negative sign is retained to indicate the stability of the bound pair. as the difference between the pair’s formation energy and the formation energies of two isolated DVs:

$$wu_p2u.$$
(12)
We have calculated the formation and binding energies for different vacancy complexes, and have listed the results in Table 3.

We shall first discuss the formation energy $`u`$, and compare our results with those of Wang et al. , which were obtained using density-functional calculations coupled with an empirical potential. The advantages and disadvantages of such a coupling scheme are clearly discussed in the original paper . In the case of 1-DV vacancy on Si(001), we find the formation energy of single DV’s $`u0.19`$ eV, in good agreement with the value of 0.22 eV reported in Ref. . For the other vacancy complexes, the agreement with Ref.  is poor because the empirical model  does not accurately capture the relaxations of the under-coordinated second-layer atoms that are present in the 2-DV and 1+2-DV complexes. Due to this limitation, the Tersoff potential  predicts that the 2-DV complex is unstable with respect to the formation of two isolated single DVs, $`u_{2\mathrm{DV}}>2u_{1\mathrm{DV}}`$. We find however that the empirical potential is still able to predict that a 2-DV and a 1-DV on the same dimer row can bind to form a 1+2-DV complex (i.e. $`u_{1+2\mathrm{D}\mathrm{V}}<u_{2\mathrm{DV}}+u_{1\mathrm{DV}}`$), though the binding energy for the 1+2DV complex is about half the value obtained from . In the case of Ge(001) surface, the potential  fares better than the silicon parametrization  in describing the stability of DV complexes. Both the 2-DV and 1+2-DV vacancies on Ge(001) are stable against breaking into isolated smaller complexes, as indicated by their respective formation energies in Table 3. This result is in agreement with the observations of Yang et al., who report that 1+2-DVs and 2-DVs are predominant on the Ge(001) surface .

We now turn to a discussion of the binding energies of DV pairs $`w`$ given in Table 3. There are two main issues that complicate the comparison of theoretical results in Table 3 with the experimentally determined binding energies. First, we note that the binding energy of an isolated bound pair computed here and in Ref.  are defined in a different manner than what is usually obtained from experiments. Experiments measure the binding energy when the DVs are part of a vacancy line, while calculations address an isolated pair of vacancies. We have estimated that the binding of two vacancies (on an average) within a VL is up to $`40\%`$ stronger than the binding within an isolated bound pair. Secondly, the experiments report the binding energies for DVs on Si(001) with a certain amount of Ge coverage  or metal contamination , as opposed to clean surfaces. Given these two concerns related to experimental conditions, as well as the use of an empirical potential in this study, discrepancies between the binding energies $`w`$ listed in Table 3 and the ones reported experimentally are expected. However, we have found that there is reasonable agreement with the experimental data. For the single DVs on a Ge-covered Si(001) surface (1.5ML) Chen et al.  estimate the binding energy in the range $`0.18\mathrm{eV}<w<0.25\mathrm{eV}`$, while our result for 1-DVs on pure Ge(001) is $`0.24`$ eV and $`0.28`$ eV on pure Si(001) (refer to Table 3). On the Ni-contaminated Si(001) surfaces, Zandvliet and coworkers  report an average attractive interaction per unit length of $`0.11`$ eV/$`a`$ ; since the distance between two neighboring dimer rows is $`2a`$, we infer that the binding energy of a pair within a VL is $`0.22`$ eV. Chang et al. report a similar figure for the binding energy ($`0.24`$ eV) of DVs on Ag-contaminated Si(001) surfaces . On both Ni- and Ag- contaminated Si(001) the vacancy lines contain almost equal proportions of 2-DVs and 1+2-DVs  and we note that the binding energies reported in  lie between the calculated (Table 3) binding energy of a 2-DV pair ($`0.11`$ eV) and that of an 1+2-DV pair ($`0.32`$ eV) on Si(001).

The vacancy formation energy $`u`$ of a vacancy complex and the binding energy $`w`$ of a pair determine the formation energy of the vacancy lines to a large extent. Chen et al.  and Weakliem et al.  have found that interactions between single DVs that belong to different dimer rows which are not adjacent can be neglected. This approximation implies that the DVs that make up a vacancy line can be treated as a 1-dimensional system with nearest-neighbor interactions. For such a system, the formation energy per unit length can be written as:

$$\mathrm{\Lambda }\frac{u+w}{2a},$$
(13)
where $`2a`$ is the distance between two DVs in a straight vacancy line. The data in Table 3 shows that the formation energies given by Eq. (13) are in good agreement with the line formation energies obtained from the $`2\times N`$ surface energies (Table 2).

The differences between the line formation energies shown in the last two columns of Table 3 are due to the presence of interactions between a given vacancy with higher-order neighbors in the same VL, which are neglected in formula (13). While such effects could be considered small, note from Table 3 that they systematically lead to a lowering of the formation energy of vacancy lines below the estimate given by (13). As seen in the next section, the DVs are more stable as part of a vacancy line than as members of a bound pair, which is due to the presence of interactions with higher-order neighbors.

### 5.3 Experimental implications

The Si(001) and Ge(001) have been investigated intensely due to their technological and fundamental importance. As noted by Zandvliet  there are still gaps in understanding the Ge(001) surface in comparison to Si(001), despite the similarity of the two surfaces. One aspect where the two surfaces show remarkable differences is that the vacancy lines appear easier to form on Si(001) rather than on Ge(001), as illustrated by the much larger number of reported observations in the case of silicon. To our knowledge, the only observation of vacancy lines on Ge(001) is for the case of Bi deposition (1ML) . No vacancy lines have so far been observed on clean Ge(001) , while two different groups report the observation of vacancy lines on the clean Si(001) surface .

We look for the differences between DVs on clean Si(001) and Ge(001) by comparing their formation and binding energies in the presence of strain. Focusing on single DVs, we compare the formation energy per vacancy for isolated vacancies, bound pairs and vacancy lines under biaxial strain. This comparison is performed separately for the case of Si and Ge surfaces. The results plotted in Fig. 4 which show that the formation energy per DV, to a good approximation, is a linear function of strain for strains in the range -2%–2%. In contrast, we have found that the binding energy $`w=u_p2u`$ exhibits strong nonlinear character in this strain range; the strains at which $`w`$ remains linear is limited to a much smaller range, $`\pm 0.2\%`$. Compressive strains lower the formation energy per DV both for the case of Si and Ge, consistent with the previous results on vacancy lines (i.e. with the positive sign of the biaxial strain sensitivities shown in Table 2. We can now identify the difference between DVs on Si(001) and Ge(001). From Table 3 note that the formation energy $`u`$ of a single DV is positive both for Ge(001) and Si(001), but in the case of Ge(001) it is about 3 times higher than that of an isolated DV on Si(001). Therefore, in the absence of any kinetic effects, a compressive strain of $`0.5\%`$ would be required to create stable DVs ($`u<0`$) on Si(001), which is readily achievable as result of quenching . In contrast, the critical value of compression needed to have single DVs on Ge(001) is four times larger, $`1.9\%`$. As seen in Fig. 4, when DVs are part of a bound pair the compression at which the formation energy per DV becomes negative decreases due to the attractive interaction between DVs: for Si(001) the critical compression becomes $`0.1\%`$, while for DVs on Ge(001) it remains rather large, $`1.5\%`$.

Based on these results (Fig. 4), we propose that organization of the DVs into vacancy lines on Ge(001) could be observed provided that the surface is under compressive strains of $`2\%`$. Biaxial strains are somewhat difficult to apply and control accurately in practice. A simple way to provide some compression would be quenching, as done in . As far as lowering the formation energy of DVs is concerned, uniaxial strains along the dimer row direction are more efficient than biaxial deformations. Therefore, the use of ultraflat, single-domain surfaces under uniaxial compression along the dimer row can lead to observations of vacancies organizing into lines on Ge(001).

## 6 Conclusions

In conclusion, we have performed an investigation of the VLs and DVs on strained Si(001) and Ge(001). We have found that the dimer-vacancies on the two surfaces exhibit differences in terms of the formation energy, their strain sensitivity, and the binding energy of two DVs in adjacent dimer rows. With the present model of atomic interactions , each of these energetic parameters appear to favor the self-organization of DVs into vacancy lines on Si(001), but not on the Ge(001) surface.

The most important difference concerns the formation energy of single DVs, which was found to be 0.55 eV for Ge(001), but only 0.19 eV in the case of Si(001). The formation energy difference between single DVs on Si and Ge surfaces seems robust with respect to the empirical potential used: we have recalculated these energies with the Stillinger-Weber model , and found similar results, though the difference between the two formation energies predicted by the Stillinger-Weber potential was larger. This finding cannot be explained by scaling the DV formation energies with the bulk cohesive energy (or melting temperature) of Si and Ge. Rather, it could likely be due to subtle details in the relaxation of the atoms that surround the vacancy, particularly to the interplay between the lateral and vertical relaxation; further studies are underway in order to elucidate this point. While it is well known that the stability of DVs on (001) surfaces is strongly enhanced by kinetic effects (i.e. two atoms are unlikely to find the DV and fill it at the same time ), it is conceivable that the large gap between the formation energies of DVs on Si(001) and Ge(001) determines their different physical behavior.

Based on studies of how the formation energy varies with external strains, we have proposed that the self-assembly of vacancies on clean Ge(001) would be possible if compressive strains are applied either biaxially or along the dimer row direction. Pending an increase in the available computational power and methodologies, the results of this article warrant future studies using higher-level total-energy methods such as tight-binding or density functional calculations. The use of such methods is expected to bring quantitative improvements of the energetic properties computed here. Furthermore, higher-level methods can also address other important aspects of the self-assembly of vacancies such as the energetic barriers associated with the creation and diffusion of DVs on Si(001) and Ge(001). Future investigations regarding the kinetics of dimer vacancies would be invaluable for fully elucidating the differences between Si(001) and Ge(001).

## 7 Acknowledgments

Support from the MRSEC at Brown University (DMR-0079964), National Science Foundation (grants CMS-0093714, CMS-0210095), and the Salomon Research Award from the Graduate School at Brown University is gratefully acknowledged. We thank H.J.W. Zandvliet for useful discussions on the Ni-induced vacancy lines , and for bringing Ref.  to our attention.


format: text

text:
# Role of the initial conditions on the enhancement of the escape time in static and fluctuating potentials

## 1 Introduction

In the two last decades a remarkable increase has been done in the study of noise induced effects on nonlinear nonequilibrium systems, with attention also to complex and biological systems . The introduction of noise and of a deterministic driving force in such a systems gives rise to counterintuitive effects which explain apparently anomalous behaviors in experiments. Nonlinear relaxation decay of physical systems from an initial unstable or metastable state involves fundamental aspects of non-equilibrium statistical mechanics. Examples of these resonance-like phenomena are stochastic resonance (SR) , resonant activation (RA) , noise enhanced stability (NES) , noise induced phase transitions (see ref.), etc… A nonmonotonic behavior of the average escape time as a function of the noise intensity was revealed in a numerical study of a Brownian particle moving in a periodic fluctuating cubic potential . This is the NES phenomenon: the stability of an otherwise unstable system can be enhanced by the presence of a finite amount of noise. The nonmonotonic behavior that often occurs in NES effect contradicts standard Kramers-like behavior, i.e. an exponential or monotonic decrease of the mean escape time with noise intensity . The NES phenomenon was experimentally detected in the transient dynamics of an unstable physical system  and observed in different physical systems . Recently some works have been done studying piecewise linear potentials with a metastable state, and exact evaluation of the decay time was obtained . In this work we present a study of the average decay time of an overdamped Brownian particle subject to a cubic potential with a metastable state at different unstable initial conditions. A relevant result of this work is the existence, for fixed potential, of a peculiar initial position $`x_c`$ of the Brownian particle, which corresponds to the intersection point between the potential profile and the $`x`$-axis (see Fig. 1). For all the initial unstable states between the maximum of the potential and $`x_c`$, the average escape time diverges, while for initial positions $`x_o>x_c`$ we find a nonmonotonic behavior of the same quantity. Our results are consistent with that obtained in the case of piecewise linear potential profile . For a periodic fluctuating potential we recover the NES phenomenon for different values of the amplitude and of the frequency of the periodical driving force and large fluctuations of the average escape time as a function of the amplitude and the frequency at $`x_o=x_c`$.

## 2 The Model

The starting point of our study is the following Langevin equation

$$\dot{x}=\frac{U(x,t)}{x}+\sqrt{D}\xi (t)$$
(1)
where $`\xi (t)`$ is the white Gaussian noise with the usual statistical properties: $`\xi (t)=0`$ and $`\xi (t)\xi (t+\tau )=\delta (\tau )`$, and

$$U(x,t)=ax^2bx^3xAcos(\omega t),$$
(2)
is the potential shown in Fig. 1, with $`A`$ and $`\omega `$ respectively the amplitude and the frequency of the driving force, and $`a=0.3`$, $`b=0.2`$. The potential profile has a local stable state at $`x_0=0`$ and an unstable state at $`x_0=1`$. For fixed potential $`A=0`$ the average escape time of a particle starting from $`x_o`$ and reaching a final position $`x_F`$ is given by

$$\tau (x_0,x_F)=\frac{2}{D}_{x_0}^{x_F}e^{2u(z)}_{\mathrm{}}^ze^{2u(y)}𝑑y𝑑z.$$
(3)
where $`u(y)=(0.3y^20.2y^3)/D`$ is a dimensionless potential profile. We evaluate this double integral partly analytically and partly numerically, splitting it as

$$\tau (x_0,x_F)=\frac{2}{D}_{x_0}^{x_F}e^{2u(z)}𝑑z\left[_{\mathrm{}}^0e^{2u(y)}+_0^ze^{2u(y)}𝑑y\right],$$
(4)
and evaluating analytically the first term inside the square parenthesis:

$$_{\mathrm{}}^0e^{2u(y)}=0.6046e^z[I_{1/3}(z)+I_{1/3}(z)]\frac{1}{2}{}_{2}{}^{}F_{2}^{}(\frac{1}{2},1;\frac{2}{3},\frac{4}{3};2z),$$
(5)
where $`z=1/(10D)`$, $`I_n(z)`$ is the modified Bessel function of the first kind and $`{}_{p}{}^{}F_{q}^{}(a_1,a_2;b_1,b_2;z)`$ is the generalized hypergeometric function.

## 3 Results and Comments

For all the initial unstable states beyond the potential barrier (see Fig. 1) we find an enhancement of the average escape time with respect to the deterministic time as a function of the noise intensity (Fig. 2). Specifically, for high values of the noise intensity with respect to the height of the barrier, we recover the Kramers behavior. For intermediate values of noise intensity, we observe an increase of the MFPT whose maximum goes to infinity when the initial position is in the range $`x_{max}<x_o<x_c`$. These results are quite consistent with those obtained in the case of piecewise linear potential profiles , where a divergent behavior for $`D0`$ is found. On the other hand for very low values of noise intensity we obtain the deterministic decay time. This effect, which appears in contradiction with the divergence for $`D0`$, can be explained with the finite ensemble of particles in numerical experiments. Indeed, for low noise intensity the probability that some particle is trapped is very low and decreases exponentially to zero, that is only very few particles will be trapped. Only these particles, which represent very rare events for very small noise values, contribute to the enhancement of the escape time and give rise to the divergence. Therefore, because of this, we do not observe such particles in simulation and the average escape time becomes equal to the deterministic time. In Fig. 2 we can clearly see this effect. For initial positions approaching the maximum of the potential this effect is less pronounced, because the trapping probability of the particle is proportional to $`exp(\mathrm{\Delta }u(x_o)/D)`$, where $`\mathrm{\Delta }u(x_o)`$ is the potential barrier “seen” from the particle at $`x(0)=x_o`$. When the particle is near the maximum a relatively low amount of noise can push the particle back into the “stable” state, where it remains trapped for a long time, because of the low noise intensity. For higher values of initial position the amount of noise to put the particle into the potential well must be higher. This explains why the maximum of the average decay time is shifted towards higher values of the noise intensity. The observation time in digital simulation is finite, therefore the simulated escape time takes the maximum at this limiting point, and this is what we find (see Fig. 2a). For all the initial conditions within the range: $`x_m<x_o<x_c`$, we find that the average escape time can be more and more increased by the noise, while for initial conditions $`x_o>x_c`$ we obtain nonmonotonic behavior with a maximum and a finite average escape time when the noise intensity goes to zero. We calculate also the standard deviation ($`\sigma _\tau `$) of MFPT as a function of the noise intensity for all the initial positions investigated (see Fig.2b). For $`x_o>x_c`$ we obtain a nonmonotonic behavior like the MFPT behavior. For $`x_o<x_c`$ the standard deviation grows towards divergent values for $`D\mathrm{}`$, and this corresponds to a big tail in the first passage time distribution. From the above analysis it appears that a peculiar initial position for the potential profile of Eq.(2) is $`x_o=x_c=1.5`$. We calculate then the theoretical average escape time from Eqs.(4,5) for an initial position just after the cross point $`x_c`$ and we show the result in Fig. 3. The agreement between the numerical simulation of the SDE of the Eq. (1) and the theoretical evaluation is very good.

We calculate the average escape time at the static cross position $`x_c=1.5`$ for a periodical fluctuating cubic potential and the results are shown in Fig.5 and Fig.5. Here we recover the enhancement of stability of the metastable state for all values of the amplitude and of the frequency of the periodical driving force investigated, with a nonmonotonic behavior of MFPT (a-figures). We find also that (b-figures): (i) the MFPT values are almost independent from the amplitudes and the frequencies of the driving force, except for the case $`D=0.05`$; (ii) very large fluctuations of MFPT as a function of both the amplitude and the frequency of the periodic driving force at $`D=0.05`$ are present. The origin of these large fluctuations can be ascribed to the particular initial position represented by $`x_c`$, which is one of the two boundaries delimiting the parameter region (A,$`\nu `$) where NES effect can be observed . Because of the oscillation of the potential between two limiting curves the intersection point is a function of the time $`x_c(t)`$ and as a consequence the Brownian particle experiences, during each period of the driving force, the two different dynamical regimes of the case of the static potential discussed above for $`D0`$. Specifically the escape time passes from a divergent regime ($`x_c(t)<x_c`$) to a convergent one ($`x_c(t)>x_c`$), where $`x_c`$ is the intersection point for the static potential. This effect becomes more important for low noise intensity, i.e. for noise intensities less than the barrier height, as we can see in Figs.5 and 5 for $`D=0.05`$. We obtain the same large fluctuations for other values of $`D`$ less than the barrier height (0.1), not reported here for sake of clearness of the figures. For values of amplitude and frequencies of the driving force higher then those presented in the figures, up to $`\nu =20`$ and up to $`A=3`$ respectively, the calculations confirm these large fluctuations. We also investigated the influence of different initial phases on the MFPT behavior as a function of the frequency and we recover large fluctuations again for $`D=0.05`$.

## 4 Conclusions

Nonmonotonic behavior of the mean escape time as a function of noise intensity is a noise-induced effect for nonlinear nonequilibrium systems with metastable states. In this work we analyzed the role of the initial conditions on the enhancement of the escape time from initial unstable states for a cubic potential. We obtain NES effect for static and periodical fluctuating potential and an enhancement of the NES effect for initial positions between the maximum of the potential well and the cross point $`x_c`$. We find also large fluctuations of the average escape time as a function of both the amplitude and the frequency of the sinusoidal force for low noise intensity, due to the different dynamical regimes experienced by the Brownian particle. Our results obtained for a particle moving in a cubic potential are quite general, because we always obtain NES effect when a particle is initially located just to the right of a local potential maximum, with a local minimum in its left side and the global escape region in its right side.

## 5 Acknowledgement

We thank Dr. A. La Barbera for useful discussion and suggestions concerning the numerical simulations. This work was supported by INFM, MIUR and by INTAS Grant 01-0450.


text:
# Volume-limited SDSS/FIRST quasars and the radio dichotomy

## 1. Introduction

Quasars were first found as optical identifications of luminous radio sources. However, only about 10% of optically identified quasars were radio-luminous, leading to a division of quasars into radio-“loud” and radio-“quiet” objects. Radio observations of optically selected quasars (Strittmatter et al. 1980 and papers citing them) found a bimodal distribution of the radio flux, radio luminosity, or ratio of radio to optical flux. This has resulted in many claims that there are two distinct populations, although it is still unclear whether they should be distinguished by considering the radio luminosity or the radio-optical flux ratio, and what the physical origin of this bimodality or dichotomy might be. More recently, the FIRST radio survey has filled gaps in the radio-to-optical flux ratio distribution found in earlier, shallower surveys (White et al 2000), suggesting that the radio-loud and radio-quiet objects are instead the extremes of a continuum of sources.

The SDSS quasar survey targets both optically selected quasar candidates and point-like optical counterparts of FIRST sources. Together with the large number of sources compared to previous surveys, this allows to construct a series of volume-limited samples. This avoids any potential biases arising from the rapid evolution of quasar counts and luminosities with redshift. Therefore, the SDSS allows to take a different look at the bimodality question. Once selection effects are taken into account, this approach amounts to constructing the bivariate radio-optical quasar luminosity function.

## 2. Results

We use data (including rest-frame $`i`$-band luminosities) from the quasar catalog (Schneider et al. 2003; also Schneider et al. in this volume) created from the SDSS Data Release 1 (Abazajian et al. 2003). We restrict ourselves to FIRST-selected quasar candidates and those identified by the “low-redshift UV excess” ($`ugri`$) part of the color selection algorithm (Richards et al. 2002). This algorithm selects targets with non-stellar colors down to $`i=19.1`$ and has been sufficiently uniform and complete for DR1 data so that first results can be obtained without any corrections.

We construct a series of volume-limited samples of radial extent $`\mathrm{\Delta }z=0.2`$, retaining only objects with luminosities between that corresponding to the faint flux limit at the upper redshift boundary ($`i=19.1`$ for SDSS, $`f_{1.4\mathrm{GHz}}=1\mathrm{mJy}`$ for FIRST) and that corresponding to the bright limit at the lower redshift boundary ($`i=15`$ for SDSS; FIRST practically has no bright limit).

### 2.1. Relation between optical and radio luminosity

We construct volume-limited quasar samples between redshifts of 0.2 and 3.2. This retains a few hundred objects in each redshift shell out of originally 10,000. Figure 1 shows an example of the distribution of quasars making the luminosity cut of both surveys in the luminosity-luminosity plane. They straddle the dividing line between radio-loud and radio-quiet objects as defined by a radio-optical luminosity ratio of 10 (dashed line), with many objects lying close to the line. A caveat in the interpretation of this result is that two simplifications have been made: the $`K`$-correction did not take into account the actual spectral index of the quasar, but assumed a spectral shape $`f_\nu \nu ^{0.5}`$ for all objects at all wavelengths (compare the contribution about quasar bimodality by Ivezić et al. in these proceedings), and the radio sources include both flat- and steep-spectrum sources, while in principle we need to remove intrinsically fainter flat-spectrum sources which have been beamed above the survey limit. These corrections will be made in future work.

In this particular case, the most luminous radio source is also the most luminous optical source, but this is not generally true - the maximum optical and radio luminosity in each bin are not well correlated. Similarly, the radio and optical luminosities of objects within each bin are not strongly correlated, with correlation coefficients of the logarithmic luminosities typically in the range 0-0.5. This rules out any strict power-law dependence between radio and optical luminosity with a scatter smaller than the luminosity range in each redshift shell (which spans, however, only one decade or so in optical luminosity). Instead, it requires a hidden parameter governing the radio luminosity relative to the optical.

In addition to the objects shown in Fig. 1, each redshift shell also contains quasars in the same optical luminosity range but which have a FIRST luminosity below the cutoff. It is possible to compare the radio properties of all SDSS-detected quasars by constructing optical luminosity functions cumulated up to different maximum radio luminosities. The optical luminosity function becomes steeper as the maximum radio luminosity limit is lowered from the brightest detected source to the luminosity cutoff, i.e., a luminous radio source is more likely (but not *required*) to be a luminous optical source as well.

### 2.2. Radio-loud fraction as function of redshift

In view of the strong evolution of quasar counts with redshift, it is of considerable interest whether the fraction of radio-loud objects changes with redshift. (The magnitude of the fraction depends on the ratio of the radio and optical flux limits and cannot in itself be interpreted easily.) In the case of SDSS, it is not sufficient to consider merely the fraction of FIRST detections among all quasars in each volume slice without completeness corrections because the survey employs sparse sampling of the color-selected candidates in certain redshift ranges (Richards et al. 2002), while the number of quasars targeted as point-source identification of FIRST sources remains unchanged (lower panel in Fig. 2). Ignoring this would lead to the erroneous conclusion that the fraction of radio-detected quasars increases drastically beyond $`z=2`$.

Instead, it is useful to consider the fraction of color-selected quasars (both resolved and point sources) which have a FIRST detection (upper panel in Fig. 2). This fraction is remarkably constant, implying a similar evolution of radio and optical luminosity with redshift. This finding is in apparent contradiction to the preceding section, where we found no good radio-optical correlation on an object-by-object basis. The reconciliation lies in a radio-optical correlation with a scatter larger than the luminosity range we have considered in each redshift shell. It remains to investigate these correlations in more detail and to compare them to radio-optical correlations such as those found in the most powerful radio sources (see, e.g., Willott et al. 1999).

We are grateful to Željko Ivezić for challenging discussions, to Matt Jarvis for valuable suggestions, and to the SDSS Quasar Working Group which has laid the foundation for this work. Funding for the Sloan Digital Sky Survey (SDSS) has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/.

## References

Abazajian, K., et al. 2003, AJ, 126, 2081

Richards, G. T., et al. 2002, AJ, 123, 2945

Schneider, D. P., et al. 2003, AJ, accepted \[astro-ph/0308443\]

Strittmatter, P. A., et al. 1980, A&A, 88, L12

White, R. L., et al. 2000, ApJS, 126, 133

Willott, C. J., et al. 1999, MNRAS, 309, 1017


text:
# First observation of 𝝍(𝟐𝑺)→𝑲_𝑺^𝟎𝑲_𝑳^𝟎

## Abstract

The decay $`\psi (2S)K_S^0K_L^0`$ is observed for the first time using $`\psi (2S)`$ data collected with the Beijing Spectrometer (BESII) at the Beijing Electron Positron Collider (BEPC); the branching ratio is determined to be $`(\psi (2S)K_S^0K_L^0)=(5.24\pm 0.47\pm 0.48)\times 10^5`$. Compared with $`J/\psi K_S^0K_L^0`$, the $`\psi (2S)`$ branching ratio is enhanced relative to the prediction of the perturbative QCD “12%” rule. The result, together with the branching ratios of $`\psi (2S)`$ decays to other pseudoscalar meson pairs ($`\pi ^+\pi ^{}`$ and $`K^+K^{}`$), is used to investigate the relative phase between the three-gluon and the one-photon annihilation amplitudes of $`\psi (2S)`$ decays.

preprint: Draft-PRL

It has been determined that for many two-body exclusive $`J/\psi `$ decays, like vector pseudoscalar (VP), vector vector (VV), and pseudoscalar pseudoscalar (PP) meson decays and nucleon anti-nucleon (N$`\overline{\text{N}}`$) decays, the relative phases between the three-gluon and the one-photon annihilation amplitudes are near $`90^{}`$ suzuki ; dm2exp ; mk3exp ; a00 ; a11 ; ann . For $`\psi (2S)`$ decays, the available information about the phase is much more limited because there are fewer experimental measurements. It has been argued that the relative phases in $`\psi (2S)`$ decays should be similar to those in $`J/\psi `$ decays suzuki ; gerard , but the analysis of $`\psi (2S)`$ to VP decays in Ref. suzuki  indicates this phase is likely to be around $`180^{}`$. Another analysis of this mode though shows the relative phase observed in $`J/\psi `$ decays could also fit these decays wymphase , but it could not rule out the $`180^{}`$ possibility due to the big uncertainties in the experimental data. Therefore it is important to measure phases in other $`\psi (2S)`$ decay modes.

In $`\psi (2S)\text{PP}`$, the currently available measurements on $`\pi ^+\pi ^{}`$ and $`K^+K^{}`$ are from DASP dasp , with huge errors. These do not provide enough information to extract the phase since there are three free parameters in the parametrization of the $`PP`$ amplitudes a00 ; a11 ; haber . The result for $`\psi (2S)K_S^0K_L^0`$ is also needed to get the phase phase\_pp .

Furthermore, there is a longstanding “$`\rho \pi `$ puzzle” between $`J/\psi `$ and $`\psi (2S)`$ decays in some decay modes; many $`\psi (2S)`$ decay channels compared with the corresponding $`J/\psi `$ decays are suppressed relative to the perturbative QCD predicted “12% rule” beswangwf . It is of great interest to check this in more channels.

In this letter, the first observation of $`\psi (2S)K_S^0K_L^0`$ is reported, and its branching ratio is used to determine the phase between the three-gluon and one-photon annihilation amplitudes and to test the “12% rule” between $`J/\psi `$ and $`\psi (2S)`$ decays. The data used for the analysis are taken with the Beijing Spectrometer (BESII) detector at the Beijing Electron-Positron Collider (BEPC) storage ring at a center-of-mass energy corresponding to $`M_{\psi (2S)}`$. The data sample contains a total of $`14(1\pm 5\%)\times 10^6`$ $`\psi (2S)`$ decays, as determined from inclusive $`\psi (2S)`$ hadronic decays pspscan .

BES is a conventional solenoidal magnet detector that is described in detail in Ref. bes , BESII is the upgraded version of the BES detector bes2 . A 12-layer vertex chamber (VC) surrounding the beam pipe provides trigger information. A forty-layer main drift chamber (MDC), located radially outside the VC, provides trajectory and energy loss ($`dE/dx`$) information for charged tracks over $`85\%`$ of the total solid angle. The momentum resolution is $`\sigma _p/p=0.017\sqrt{1+p^2}`$ ($`p`$ in $`\text{GeV}/c`$), and the $`dE/dx`$ resolution for hadron tracks is $`8\%`$. An array of 48 scintillation counters surrounding the MDC measures the time-of-flight (TOF) of charged tracks with a resolution of $`200`$ ps for hadrons. Radially outside the TOF system is a 12 radiation length, lead-gas barrel shower counter (BSC). This measures the energies of electrons and photons over $`80\%`$ of the total solid angle with an energy resolution of $`\sigma _E/E=22\%/\sqrt{E}`$ ($`E`$ in GeV). Outside of the solenoidal coil, which provides a 0.4 Tesla magnetic field over the tracking volume, is an iron flux return that is instrumented with three double layers of counters that identify muons with momentum greater than 0.5 GeV/$`c`$.

A Monte Carlo simulation is used for the determination of the detection efficiency. For the signal channel, $`\psi (2S)K_S^0K_L^0`$, the angular distribution of the $`K_S^0`$ or $`K_L^0`$ is generated as $`\mathrm{sin}^2\theta `$, where $`\theta `$ is the polar angle in the laboratory system. The $`K_L^0`$ is allowed to decay and interact with the material in the detector, and for the $`K_S^0`$, only $`K_S^0\pi ^+\pi ^{}`$ is generated. The detector response is simulated using a Geant3 based Monte Carlo program, SIMBES. Reasonable agreement between data and Monte Carlo simulation has been observed in various channels tested, including $`e^+e^{}(\gamma )e^+e^{}`$, $`e^+e^{}(\gamma )\mu ^+\mu ^{}`$, $`J/\psi p\overline{p}`$ and $`\psi (2S)\pi ^+\pi ^{}J/\psi ,J/\psi \mathrm{}^+\mathrm{}^{}`$ $`(\mathrm{}=e,\mu )`$.

Candidate events are required to satisfy the following selection criteria:

1. The number of charged tracks is required to be two with net charge zero.
2. Each track should satisfy $`|\mathrm{cos}\theta |<0.8`$, where $`\theta `$ is the polar angle in the MDC, and should have good helix fit so that the error matrix from track fitting is available for secondary vertex finding.
3. $`E_\gamma ^{tot}<1.0`$ GeV, where $`E_\gamma ^{tot}`$ is the total energy of the neutral clusters in the BSC which are not associated with the charged tracks.

The two tracks are assumed to be $`\pi ^+`$ and $`\pi ^{}`$. To find the intersection of the two tracks, an iterative, nonlinear least squares technique is used. The intersection is taken as the $`K_S^0`$ vertex, and the momentum of the $`K_S^0`$ is calculated at this point. Figure 1 shows a scatter plot of the $`\pi ^+\pi ^{}`$ invariant mass versus the decay length in the transverse plane ($`L_{xy}`$) for events that satisfy the above selection criteria and have $`\pi ^+\pi ^{}`$ momentum greater than 1.7 GeV/$`c`$. The cluster of events with mass consistent with the nominal $`K_S^0`$ mass and with a long decay length indicates a clear $`K_S^0`$ signal. The lack of events at $`L_{xy}>0.1`$ m is due to the trigger, which will be discussed later.

Fits of the $`\pi ^+\pi ^{}`$ invariant mass distributions (not shown) indicate good agreement between data and Monte Carlo simulation in both mass and the mass resolution. After requiring $`L_{xy}>1`$ cm and $`\pi ^+\pi ^{}`$ mass within twice the mass resolution around the $`K_S^0`$ nominal mass, the $`K_S^0`$ momentum distribution, shown in Figure 2, is obtained. In the plot, there is a clear peak at around 1.77 GeV/$`c`$ with low background, as indicated by the $`K_S^0`$ mass side band (three sigma away from the $`K_S^0`$ nominal mass on both sides) events. The excess at lower momentum, which is not explained by the side band background, is due to the contribution of the background channels with $`K_S^0`$ production.

The main $`K_S^0`$ production backgrounds near the signal region are from $`\psi (2S)K^0(892)\overline{K^0}+c.c.`$ and $`\psi (2S)\gamma \chi _{c1},\chi _{c1}K^0(892)\overline{K^0}+c.c.`$ where the $`K^0`$ decays into $`K^0`$ and $`\pi ^0`$ and one of the $`K^0`$s becomes a $`K_S^0`$ and the other one becomes a $`K_L^0`$. The background from $`\psi (2S)\gamma \chi _{cJ}(J=0,2),\chi _{cJ}K_S^0K_S^0`$ where one of the $`K_S^0`$ decays into $`\pi ^+\pi ^{}`$ and the other decays into $`\pi ^0\pi ^0`$ is removed almost entirely (more than 95%) by the $`E_\gamma ^{tot}`$ cut, according to the Monte Carlo simulation. The decay $`\psi (2S)J/\psi +X`$ with $`J/\psi `$ decaying into $`K_S^0X`$ has a big branching ratio, but the $`K_S^0`$ momentum is much lower. The light shaded histogram in Figure 2 shows the contribution of the background channels normalized to the known branching ratios beschic ; pdg . It can be seen that the agreement between the background estimation and data is good near the $`K_S^0K_L^0`$ peak, indicating that the estimation of the background under the $`K_S^0K_L^0`$ peak is reliable.

Under SU(3) symmetry, $`K_S^0K_L^0`$ production via virtual photon annihilation is forbidden. This is checked by applying the same selection criteria to the data sample taken below the $`\psi (2S)`$ peak, at $`\sqrt{s}=3.650`$ GeV, with an integrated luminosity of about one-third of that at the $`\psi (2S)`$ peak. Figure 3 shows the $`K_S^0`$ momentum spectrum of the selected events; the events in the signal region agree well with expectation from the $`K_S^0`$ mass side band events. Taking all four candidates with momentum between 1.7 and 1.9 GeV/$`c`$ as signal, the upper limit of the production cross section at the 90% C. L. is measured to be $`\sigma <5.9\text{pb}`$. The background from the continuum contribution is thus neglected in the folowing analysis since no evidence for $`K_S^0K_L^0`$ production via the virtual photon process is observed.

The $`K_S^0`$ momentum spectrum of the selected events is fitted from 1.45 to 2.0 GeV/$`c`$ with a Gaussian distribution for the signal and an exponential for the background using the unbinned maximum likelihood method. The result is shown in Figure 2. The backgrounds from the $`K_S^0`$ mass side bands and the simulated background channels are also shown, and they agree with the fitted background reasonably well near the signal region, considering the uncertainties in the global normalization of the background channels. The peak $`K_S^0`$ momentum is $`(1775.0\pm 3.3)`$ MeV/$`c`$, which is in good agreement with the Monte Carlo expectation. The momentum resolution also agrees with the Monte Carlo simulation well. The fit yields $`n^{obs}=156\pm 14`$ events, and the efficiency for detecting $`\psi (2S)K_S^0K_L^0`$, with $`K_S^0\pi ^+\pi ^{}`$ is $`\epsilon _{MC}=(41.59\pm 0.48)\%`$ from the Monte Carlo simulation, where the error is due to the limited statistics of the Monte Carlo sample.

Due to the long decay length of the high momentum $`K_S^0`$ particles and the trigger requirement for hits in the VC, the trigger efficiency of $`K_S^0K_L^0`$ events is very different than for normal hadronic events. Since the trigger system is not included in the Monte Carlo simulation, the trigger efficiency is measured using $`K_S^0K_L^0`$ events by comparing the number of events beyond and within the outer radius of the VC with what would be expected for an exponential decay, which yields a trigger efficiency of $`\epsilon _{trig}=(76.0\pm 1.8)\%`$. The systematic error in the branching ratio measurement from this source and all other sources are listed in Table. 1.

The efficiency of the secondary vertex finding, $`\epsilon _{2nd}`$, is studied using $`J/\psi K^{}(892)\overline{K}+c.c.`$ events. It is found that the Monte Carlo simulates data fairly well. Extrapolating the difference between data and Monte Carlo simulation to the $`K_S^0`$ momentum range under study and correcting by the polar angle dependence of the efficiency, a correction factor of $`(98.1\pm 4.0)\%`$ is obtained to the Monte Carlo efficiency.

The effect of the cut on the total energy of the photon candidates is checked with $`J/\psi K_S^0K_L^0`$ events. The efficiency difference between data and Monte Carlo simulation is measured to be $`0.99\pm 0.01`$; data is slightly lower than Monte Carlo simulation. No correction to the final efficiency is made, and 2% is taken as the systematic error on the efficiency associated with this cut.

The simulation of the tracking efficiency agrees with data within 1-2% for each charged track as measured using channels like $`J/\psi \mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ and $`\psi (2S)\pi ^+\pi ^{}J/\psi `$, $`J/\psi \mu ^+\mu ^{}`$. The systematic error for the channel of interest is taken conservatively as 4% .

The Monte Carlo simulated mass resolution and momentum resolution of the $`K_S^0`$ agree with those determined from data within the statistical uncertainties. The cut on the $`K_S^0`$ mass distributions at 2 standard deviations introduces a very small systematic bias and is neglected.

Varying the lower and upper bounds of the fitting range results in a 2.4% change in the number of the events; using a second order polynomial for the background parametrization causes a 3% change in the number of events. These are taken as systematic errors. The systematic error in the total number of $`\psi (2S)`$ events, $`N_{\psi (2S)}`$, which is measured using inclusive hadrons in the same way as in Ref. pspscan , is taken as 5%. The systematic error on the branching ratio $`(K_S^0\pi ^+\pi ^{})`$ is obtained from the PDG pdg  directly.

Figure 4 shows the cosine of the $`K_S^0`$ polar angle for $`K_S^0K_L^0`$ events from $`\psi (2S)`$ decays; agreement between data and Monte Carlo simulation is observed. This distribution is also checked with a larger sample from $`J/\psi K_S^0K_L^0`$, where the Monte Carlo simulation agrees with data very well. This indicates that the angular distribution used in the Monte Carlo generator is correct.

The branching ratio of $`\psi (2S)K_S^0K_L^0`$ is calculated with

$$(\psi (2S)K_S^0K_L^0)=\frac{n^{obs}/(\epsilon _{MC}\epsilon _{trig}\epsilon _{2nd})}{N_{\psi (2S)}(K_S^0\pi ^+\pi ^{})}.$$

Using numbers from above (listed in Table. 2), one obtains

$$(\psi (2S)K_S^0K_L^0)=(5.24\pm 0.47\pm 0.48)\times 10^5,$$

where the first error is statistical and the second systematic.

Comparing with the corresponding branching ratio for $`J/\psi K_S^0K_L^0`$ from BESII ($`(1.82\pm 0.14)\times 10^4`$bes2ksklj , one gets

$$Q_h=\frac{(\psi (2S)K_S^0K_L^0)}{(J/\psi K_S^0K_L^0)}=(28.8\pm 3.7)\%,$$

where the common errors in the $`J/\psi `$ and $`\psi (2S)`$ analyses are removed in the calculation. The result indicates that $`\psi (2S)`$ decays are enhanced by more than 4$`\sigma `$ relative to the “12% rule” expected from pQCD, while for almost all other channels where the deviations from the “12% rule” are observed, $`\psi (2S)`$ decays are suppressed.

The branching ratio of $`K_S^0K_L^0`$, together with branching ratios of $`\psi (2S)\pi ^+\pi ^{}`$ and $`\psi (2S)K^+K^{}`$, can be used to extract the relative phase between the three-gluon and the one-photon annihilation amplitudes of the $`\psi (2S)`$ decays to pseudoscalar meson pairs. It is found that a relative phase around $`\pm 90^{}`$ can explain the experimental results phase\_pp .

In summary, the flavor SU(3) breaking process $`K_S^0K_L^0`$ is observed for the first time in $`\psi (2S)`$ decays with the BESII $`\psi (2S)`$ data sample, and the branching ratio is determined to be $`(\psi (2S)K_S^0K_L^0)=(5.24\pm 0.47\pm 0.48)\times 10^5`$. Compared with the branching ratio of $`J/\psi K_S^0K_L^0`$, $`\psi (2S)`$ decays are enhanced relative to the “12%” pQCD prediction. The phases of the three-gluon and the one-photon annihilation amplitudes of $`\psi (2S)`$ decays to pseudoscalar meson pairs are found to be nearly orthogonal.

The BES collaboration thanks the staff of BEPC for their hard efforts. This work is supported in part by the National Natural Science Foundation of China under contracts Nos. 19991480, 10225524, 10225525, the Chinese Academy of Sciences under contract No. KJ 95T-03, the 100 Talents Program of CAS under Contract Nos. U-11, U-24, U-25, and the Knowledge Innovation Project of CAS under Contract Nos. U-602, U-34(IHEP); by the National Natural Science Foundation of China under Contract No. 10175060 (USTC); and by the Department of Energy under Contract No. DE-FG03-94ER40833 (U Hawaii).


text:
# XDR Chemistry in the Circumnuclear Disk of NGC 1068

## References

Fuente, A., Rodríguez-Franco, A., García-Burillo, S., et al. 2003, A&A, 406, 899

Lepp, S., & Dalgarno, A. 1996, A&A, 306, L21

Schinnerer, E., Eckart, A., Tacconi, L. J. et al. 2000, ApJ, 533, 850

Smith, M. A., Schlemmer, S., von Richthofen, J., & Gerlich, D. 2002, ApJ, 578, L87

Sternberg, A., Genzel, R., & Tacconi, L. J. 1994, ApJ, 436, L131

Tacconi, L. J., Gallimore, J. F., Genzel, R., et al. 1997, Ap&SS, 248, 59

Tacconi, L. J., Genzel, R., Blietz, M., et al. 1994, ApJ, 426, L77

Usero A., García-Burillo, S., Fuente, A., et al. 2003, A&A, submitted.


text:
# Strongly Coupled Plasmas in High-Energy Physics

## I Basic Facts of Quark, Gluons, and Quantum Chromodynamics

One of the main tasks of modern nuclear and high-energy physics is the quest of the quark-gluon plasma (QGP), a new state of matter, predicted to exist at high temperatures and densities. Quarks have been proposed in the middle of the sixties as elementary substructure particles of the strongly interacting particles, the so called hadrons, in order to explain the large spectrum of hadrons discovered at this time . Quarks are fermions, carrying the spin 1/2. Nowadays, 6 different quark species, characterized by the quantum number flavor (up, down, strange, charm, bottom, top) are known. These quark flavors cover a wide regime of masses, from about 5 MeV for the up-quark to 170 GeV for the top-quark. In addition, quarks have an electric charge, either 2/3 or -1/3 of the electron charge.

There are two different combinations of quarks for building hadrons: baryons, i.e. hadrons which are fermions, e.g. protons and neutrons, are fomed by 3 quarks, whereas mesons, which are bosons, e.g. pions, consist of a quark and an anti-quark. Atomic nuclei contain only protons and neutrons, which are build up of up- and down-quarks. The other (heavier) quark flavors exist only in short-living hadrons, created in accelerator experiments or cosmic rays. So far no free quarks have been observed in nature, a phenomenon called confinement.

For explaining the hadron spectrum as well as cross sections for the hadron production in electron-positron annihilations, additional quark states had to be introduced. These states are characterized by the quantum number color, which exists in three different possibilities (red, blue, green). In the middle of the seventies, it was realized, that this quantum number can be regarded as the charge of the strong interaction between the quarks. This led to the formulation of a theory for the quark interaction, called Quantum Chromodynamics (QCD) , analogously to QED, describing the electromagnetic interaction. For example, in QED the scattering of two electrons is mediated by the exchange of a photon. The exchange particles of QCD, mediating the strong interaction, are called gluons. The main difference between QED and QCD is the fact that quarks come in three different colors, whereas electrons carry only one charge. Mathematically this implies that QED is based on the abelian $`U(1)`$ gauge group, while QCD on the non-abelian $`SU(3)`$. As a consequence, gluons carry a color charge, in contrast to photons which have no electric charge. Therefore, gluons can interact directly with other gluons. A consequence of this gluon self-interaction is the fact that, in contrast to QED, the interaction strength of two color charges increases with the distance between these partons. (A parton is either a quark or gluon.) Then, for instance, a meson can be pictured as a rubber band with a quark on one end and an anti-quark on the other. If the quarks are close together, they do not feel each other. This phenomenon is called asymptotic freedom, which can be derived directly from QCD. If the distance between the quark and the anti-quark gets larger, however, the force between them increases, until the rubber band breaks into two new mesons, e.g. in high-energy scattering experiments. In this way, one can explain confinement: the interaction between hadrons can lead only to new hadrons but not to freely propagating quarks. In other words, only colorless states, the hadrons, can exist as free particles. However, it should be noted that confinement has not been proved rigorously from first principles.

Following this picture, hadrons can be considered as quark bags in which the quarks are confined by their interaction with each other. Since the interaction is mediated by gluons which also can split into virtual quark anti-quark pairs, the hadrons are complicated systems containing quarks and gluons (see Fig.1). A nucleus is now represented by a dense system of individual quark bags, the protons and neutrons.

## II Quark-Gluon Plasma: Experiment and Theory

Shortly after quarks were proposed the question arose, what will happen if we compress the nucleus, i.e. increase the nuclear or baryon density, or if we heat up the nucleus in a fixed volume to such high temperatures that pair creation (e.g. of pions) occurs in scattering events of hadrons with large thermal energies. Then the system becomes denser and the bags should start to overlap and to dissolve, transforming the nucleus into a “soup” of freely propagating quarks and gluons as sketched in Fig.2, which was called a quark-gluon plasma . Theoretical estimates predicted a critical baryon density of about 10 times normal nuclear matter density, $`\rho _0=0.125`$ Gev/fm3 corresponding to $`2.2\times 10^{17}`$ kg/m3. The critical temperature $`T_c`$ of this deconfinement phase transition from hadronic matter to the QGP phase should be of the order 150 to 200 MeV (1.8 to $`2.4\times 10^{12}`$ K).

The phase diagram of nuclear matter is pictured in Fig.3, where the phase transition line in the plane given by the temperature and baryon density is shown. Also three different possibilities to cross this border line are indicated. Until about a few $`\mu `$s after the Big Bang the Early Universe had a temperature above the critical temperature and should have existed in the QGP phase. In the center of neutron stars, on the other hand, the baryon density might exceed the critical one at low temperatures. Then the interior of these compact stars might actually consist of degenerate quark matter instead of neutron matter111As a matter of fact, recent investigations revealed that there might be a quark pairing (color superconductivity) at high quark densities, which could have an influence on the properties of neutron stars containing quark matter (see e.g. Ref.).. Finally in high-energy nucleus-nucleus collisions, also called relativistic heavy-ion collisions, a dense and hot fireball is created for a short time, which might be also in the QGP phase .

Different accelerator experiments have been and will be performed for searching for the QGP. For producing a large and hot fireball, heavy-ion collisions with Au-, Pb-, and U-beams at center of mass energies from 4 to 5500 GeV per nucleon have been or will be conducted at the Brookhaven National Laboratory on Long Island (AGS, RHIC) and at CERN in Switzerland (SPS, LHC). At the accelerator SPS first indications for the discovery of the QGP plasma have been found and announced at CERN in 2000 . Also in 2000 the dedicated heavy-ion experiment RHIC started its operation with energies up to 200 GeV per nucleon, which is about 10 times the energy of SPS. Due to the huge amount of data - thousands of particle tracks are produced in a single collision - the analysis of the data is still in progress, although preliminary results confirm the conclusions from SPS. The high density phase at moderate temperatures will be addressed at a future accelerator facility at the GSI Darmstadt, which eventually may allow to study the superconducting phase of QCD.

Theoretical estimates for these collisions predicted a maximum temperature clearly above the critical at RHIC ($`T_{max}>230`$ MeV), a maximum volume of the hot fireball of about 3000 fm3, containing about 10000 partons, and a life-time of the QGP of about 5-10 fm/c (1 fm/c corresponds to $`3\times 10^{24}`$ s). This short life-time is caused by the fast expansion of the fireball leading to a fast reduction of the maximum temperature and the phase transition to the hadronic state. The pre-equilibrium phase, i.e., the time needed to form an equilibrated fireball after the collision of the cold nuclei, was estimated to be of the order of only 1 fm/c due to the strong interaction between the particles. Hence, the chances to create an equilibrated QGP in relativistic heavy-ion collisions seem to be good.

However, the main problem lies in the detection of this state since the QGP, living only for an extremely short period in a small spatial volume, cannot be observed directly. Only by comparing theoretical predictions for signatures of the QGP, such as the spectra of the emitted particles (photons, leptons, and hadrons) which reach the detector after the decay of the expanding fireball, with experimental data, the QGP can be discovered. Unfortunately, so far no unambiguous signature has been proposed due to the complexity of the strongly interacting many-body system. Therefore, the final identification of the QGP might not be possible by a unique signature but rather by circumstantial evidence owing to a combination of all possible signatures.

The theoretical description of the QGP has to take into account that it is a strongly coupled system at temperatures reachable in heavy-ion collisions. At $`T200`$ MeV the strong coupling constant $`g`$, corresponding to the electromagnetic coupling $`e`$, assumes values of $`g1.5`$ \- $`2.5`$ in natural units ($`\mathrm{}=c=k_B=1`$). The strong coupling constant is temperature dependent and decreases with increasing temperature and density because of asymptotic freedom, i.e. the QGP becomes an ideal plasma in the infinitely high temperature limit. At $`T200`$ MeV, however, the QGP is a non-ideal plasma. Its Coulomb coupling parameter is given by the chromoelectrostatic interaction energy divided by the thermal parton energy222Of course, in the relativistic case, the (chromo)magnetic interaction becomes equally important.

$$\mathrm{\Gamma }=\frac{Cg^2}{\mathrm{\Delta }T},$$
(1)
where $`\mathrm{\Delta }`$ is the distance between the partons, following approximately from the number density $`nT^3`$ assuming the Stefan-Boltzmann limit (ideal gas of massless particles). At $`T=200`$ MeV the distance is of the order $`\mathrm{\Delta }=0.5`$ fm. The constant $`C=3`$ for gluons and 4/3 for quarks come from the different color charges of gluons and quarks, respectively. Typical values of $`\mathrm{\Gamma }`$ for temperatures, which can be achieved in heavy-ion collisions, are between 10 and 30, corresponding to a strongly coupled plasma. For $`T\mathrm{}`$ the coupling parameter approaches zero, although $`\mathrm{\Delta }1/T`$ since $`g(T=\mathrm{})=0`$, which demonstrates the ideal plasma limit at high temperatures.

Furthermore, the QGP is a relativistic system as the masses of its components are much smaller than their thermal energies which are of the order of the temperature. Gluons are massless, while the up- and down-quarks have masses of a few MeV. Only the mass of the strange-quark ($`m_s150`$ MeV) is comparable to the temperature. The other quark flavors have much higher masses and are therefore not thermally exited, i.e., they are suppressed compared to up-, down, and strange-quarks in the QGP.

Finally, at temperatures not too far from the critical, quantum effects play an important role in the QGP. Hence a theoretical description has to start from a statistical quantum field theory, i.e. thermal QCD. Due to the large coupling constant non-perturbative methods are required. The only well-established non-perturbative approach for quantum field theories are lattice simulations, where the quantum theoretical equations of motions are solved numerically on a 4 dimensional space-time grid using Monte-Carlo techniques . This method has the advantage that it works for all temperatures below as well as above the phase transition, allowing, for example, to determine the critical temperature, the order of the phase transition, or the equation of state of the QGP. However, it is very difficult to compute dynamical quantities such as most of the signatures of the QGP which are based on particle production rates or similar dynamical processes. Hence, perturbative methods, extended to finite temperatures and improved by resummation techniques, are also adopted to estimate dynamical properties and signatures of the QGP . Although they might be valid only at temperatures far above the critical one, these methods give at least important insights in the dynamical structure of the QGP.

In the high temperature limit ($`TT_c`$ corresponding to $`g1`$), perturbation theory can be shown to be equivalent to the semiclassical Vlasov equation , generalized to the relativistic case and the strong instead of electromagnetic interaction. The only quantum input comes from using the Bose distribution for gluons and the Fermi distribution for quarks instead of the Maxwell distribution. As in usual plasma physics one can derive, for instance, the dielectric tensor from the Vlasov equation. In the case of an isotropic and homogeneous system this tensor has only two independent components (longitudinal and transverse dielectric functions), which depend on the energy (frequency) $`\omega `$ and momentum (wave number) $`k`$:

$`ϵ_L(\omega ,k)`$ $`=`$ $`1+{\displaystyle \frac{3m_g^2}{k^2}}\left(1{\displaystyle \frac{\omega }{2k}}\mathrm{ln}{\displaystyle \frac{\omega +k}{\omega k}}\right),`$ (2)
$`ϵ_T(\omega ,k)`$ $`=`$ $`1{\displaystyle \frac{3m_g^2}{k^2}}\left[1\left(1{\displaystyle \frac{k^2}{\omega ^2}}\right){\displaystyle \frac{\omega }{2k}}\mathrm{ln}{\displaystyle \frac{\omega +k}{\omega k}}\right],`$
where $`m_g^2=g^2T^2/3(1+N_f/6)`$ is called the effective gluon mass with $`N_f`$ thermally excited quark flavors in the plasma. The same result can be obtained from the gluon self-energy in lowest order perturbation theory in the high temperature limit. Going beyond this approximation, perturbative QCD allows a systematic treatment of higher order effects, such as damping rates coming from collisions between the particles. However, in order to avoid infrared singularities due to the exchange of massless gluons and violations of the gauge invariance, resummation techniques have to be adopted at finite temperature . The development of such methods is a very active research topic.

From the chromoelectromagnetic Maxwell equations one obtains the constraints on the dielectric functions, $`ϵ_L(\omega ,k)=0`$ and $`ϵ_T(\omega ,k)=k^2/\omega ^2`$, which describe the propagation of collective plasma modes. From these equations dispersion relations $`\omega _{L,T}(k)`$ for longitudinal and transverse plasma modes in the QGP follow as shown in Fig.4. Both branches start at $`k=0`$ with the same energy $`\omega _{L,T}(k=0)=m_g`$, which can be identified with the plasma frequency. The longitudinal branch is known as plasmon. In the case of relativistic plasmas both branches are equally important, whereas in non-relativistic plasmas the transverse mode is usually neglected. Both modes approach the vacuum dispersion $`\omega =k`$ for large momenta. However, the spectral strength of the plasmon mode vanishes exponentially for large momenta, $`kgT`$. Debye screening follows from the static limit ($`\omega =0`$) of the longitudinal dielectric function describing the modification of the Coulomb potential to a Yukawa potential in the presence of a medium. The Debye screening length is given by $`\lambda _D=1/(\sqrt{3}m_g)`$. The dielectric functions of (2) also contain Landau damping for $`\omega ^2<k^2`$, for which the logarithms become imaginary. Although the plasma waves in the semiclassical approximation are located in the region $`\omega ^2>k^2`$ (see Fig.4), where no Landau damping takes place, Landau damping of virtual gluons (plasma waves) in the QGP with $`\omega ^2<k^2`$ plays an important role in higher order processes.

A further very interesting property of ultrarelativistic plasmas is the existence of collective fermion modes in analogy to the plasma modes, i.e. collective quark modes in the QGP. They can be derived within perturbation theory from the in-medium quark self-energy .

Finally, it is interesting to note that in the high temperature limit the non-abelian effects are negligible. Hence, apart from numerical factors, the above results hold also for an ultrarelativistic QED plasma, i.e. an electron-positron plasma in which the temperature is much higher than the electron rest mass. Such a situation is realized in Supernovae explosions, where temperatures up to 30 MeV are possible. The dielectric functions of such a system are given by the above equation (2), in which one only has to replace the effective gluon mass $`m_g`$ by an effective photon mass $`m_\gamma =eT/3`$.

## III Thermophoretic Flow in the Mixed Phase

As a last point we want to consider a possible mixed phase of a QGP and hadrons. Such a mixed phase would occur if the phase transition from the QGP to the hadronic phase were of first order. During the expansion of the fireball in a heavy-ion collision or in the Early Universe, the mixed phase is reached, once the temperature has dropped to the critical. Then bubbles of hadrons will form (hadronization) and grow in the QGP until the entire phase consists of hadrons. QCD lattice simulations show that the order of the phase transition depends sensitively on the values of the up-, down-, and strange-quark masses. The latest calculations show that the QCD phase transition is either of first order or continuous, with a preference for the latter case .

Assuming a first order phase transition, the mixed phase can be considered as a complex plasma, in which extended and composite objects, the hadrons, are embedded in the QGP containing point-like partons, similar as dust grains in a low-temperature plasma. The main difference lies in the fact that the dust particles are highly charged, whereas the hadrons are color-neutral objects. Now, we assume in addition to the mixed phase a temperature gradient in the fireball from the center to the surface and, for simplicity, a radial expansion of the fireball. Then, similar as in complex plasmas, a thermophoretic force will act on the hadrons, pushing them from the center to the surface as shown in Fig.5. The thermophoretic force in a complex plasma has been estimated from simple kinetic arguments and shown to be in excellent agreement with experimental results . We have generalized the kinetic derivation of the thermophoretic force to the relativistic case and calculated the final velocity of the hadrons from this effect. Even for small temperature gradients of a few MeV/fm and small life-times of the mixed phase (1-2 fm/c), we obtained large outward flow velocities of the order of 0.8 c for pions and 0.4 c for protons . These values are larger than the typical flow velocities of 0.5 c, observed in relativistic heavy-ion collisions at SPS. We take this as an indication for the absence of a mixed phase in these experiments, either due to a continuous phase transition instead of first order or due to an absence of the QGP phase at all.