text:
# Two-Dimensional Nucleation with Edge and Corner Diffusion

## I Introduction

Surface morphology during epitaxial growth depends on processes taking place on the substrate. hwang+91 ; michely+93 ; stroscio+93 ; brune+94 ; brune+95 ; hohage+96 ; li+97 ; brune+99  At very low temperatures as in the case of molecular beam epitaxy (MBE), deposited atoms are adsorbed and migrate on the substrate surface until they form nuclei or are incorporated into islands, but they never evaporate back to the ambient vapor. When an adsorbed atom (adatom, for short) sticks and freezes irreversibly to an island, the resulting aggregate takes an irregular dendritic form, hwang+91 ; michely+93 ; brune+94 ; brune+95 ; hohage+96  called diffusion-limited aggregate (DLA).witten+81  It has a characteristic feature of self-similarity, and is regarded as a typical example of a fractal. In actual experiments, one does not always observe fractal but also compact islands.michely+93 ; stroscio+93  Compact islands are possible if they can dissociate irisawa+95  or if atoms at the edge of them can diffuse along the periphery. michely+93 ; bales+94 ; zhang+94 ; bales+95 ; brune+99 ; ovesson+99 ; amar99 ; zhong+01  Since at low temperatures it is difficult for an edge atom to detach, the edge diffusion mainly governs the island morphology. brune+99  During the edge diffusion, however, the outer corner provides an additional barrier to be surpassed. There, an edge atom has to pass an intermediate state where the number of neighboring bonds is less than that on the straight step edge. The edge and the corner diffusion is known to influence island morphology. michely+93 ; bales+94 ; zhang+94 ; bales+95 ; hohage+96 ; li+97 ; amar99 ; brune+99 ; ovesson+99 ; zhong+01  Without edge and corner diffusion, islands are similar to DLA. With the edge diffusion they take cross shapes, with the corner diffusion in addition they become square. At the coalescence, these peripheral diffusion governs the shape relaxation. stoldt+98  Furthermore, the corner diffusion barrier is found to induce mound formation for multilayer deposition. opl+99 ; murty+99

Diffusion processes are also known to influence the number density $`n`$ of islands nucleated on a singular surface. venables73 ; stoyanov+81 ; venables+84 ; irisawa+90 ; villain+92a ; villain+92b ; tang93  When the deposition rate $`F`$ is large, the adatom density is high and the nucleation occurs frequently. Thus the island density $`n`$ is expected to be high. When an isolated adatom migrates on the surface with a large surface diffusion constant $`D_s`$, the deposited atom moves around a long distance, until it finds nucleation partners or preexisting islands and sticks to them. Thus, as the diffusion rate $`D_s`$ increases, the rate of forming new nuclei decreases, and the island density $`n`$ on a singular surface diminishes.

More quantitatively, the island density $`n`$ is defined as the total number density of immobile and undissociable clusters with sizes larger than the critical cluster size $`i^{}`$. In the classical mean-field nucleation theory, the number density of adatom clusters of size $`i`$ is analyzed by means of rate equations. venables73 ; stoyanov+81 ; villain+92a ; villain+92b ; furman+97  When clusters smaller than $`i^{}`$ can dissociate but only the mononer adatom is mobile with the diffusion constant $`D_s`$, $`n`$ is found to follow the scaling law venables73 ; stoyanov+81

$$n(\frac{F}{D_s})^{i^{}/(i^{}+2)}.$$
(1)
If the clusters with sizes $`i`$ smaller than the critical size $`i^{}`$ cannot dissociate but are mobile as a whole with their respective diffusion constant $`D_i`$, the scaling form is approximately given as villain+92a ; villain+92b ; furman+97

$$n(\frac{F^i^{}}{D_sD_2D_3\mathrm{}D_i^{}})^{1/(2i^{}+1)}.$$
(2)
By assuming that islands are point sinks to diffusing adatoms and the critical nuclear size is $`i^{}=1`$, the island density scaling $`n(F/D_s)^{1/3}`$ is confirmed by the kinetic Monte Carlo (KMC) simulation. bartelt+92  Actually, islands increase their size during growth, and the scaling relation between the density and the diffusion constant might depend on the island morphology. There are already many studies on the island density by means of KMC simulation when clustes can dissociate. bales+94 ; liu+95 ; brune+99 ; mulheran+01  There are also KMC simulations with mobile clusters. liu+95 ; bartelt+96 ; furman+97 ; mulheran+01  We study here the effect of edge and corner diffusion processes on the island density and its relation to the island morphology, systematically.

In $`\mathrm{\S }`$ 2, we explain our simulation scheme. In $`\mathrm{\S }`$ 3, the effect of peripheral diffusion on the island morphology is studied. On varying the ratio between the surface, edge and corner diffusion, morphology varies as fractal, needle, and square shapes, as previously observed. zhong+01  We point out the different roles of the edge and corner diffusions on the island symmetry. The island density is discussed in $`\mathrm{\S }`$ 4. As long as the corner diffusion is small, the density scaling is found to be almost independent of the island morphology with an exponent 1/3, in consistent to the more systematic study by Brune et al. brune+99  With a large corner diffusion constant, on the other hand, the scaling form of the island density changes drastically, because small clusters of sizes less than 3 can migrate on the surface, i.e. $`i^{}=3`$. The density is well fitted to the scaling $`nF^{3/7}`$. The increase of the critical island size $`i^{}`$ is also reflected in the island size distribution. bartelt+96 ; amar+95 ; furman+97  The result is summarized in a scaling plot of the island density, and discussed in the last $`\mathrm{\S }`$ 5. Even though the scaling of the island density to the deposition flux agrees to the rate equation analysis, the dependence on the diffusion constants fails to follow the expectation (2).

## II Monte Carlo simulation

Our main aim is to explore the systematics and the universal features of the morphology change of islands on a singular surface as the diffusion mechanism changes, rather than the precise reproduction of some specific experimental observations. For this purpose, we adopt the simplest model such that atoms are depositing on a (100) surface of a simple cubic lattice with a solid-on-solid (SOS) restriction. The SOS restriction is valid at low temperatures, common to MBE.

Atoms are deposited randomly onto the substrate with a deposition rate of $`F`$ monolayers per unit time. Here we choose the unit of time such that $`F=1`$. The temperature is assumed to be so low that all deposited atoms stick to the surface and never evaporate back into the ambient gas. An isolated adatom hops to one of its four nearest neighboring sites at a rate $`k_s`$; the waiting time of an adatom before the jump is $`(4k_s)^1`$, and the surface diffusion constant is given by $`D_s=k_sa^2`$, where $`a`$ is the lattice constant. As for the interlayer diffusion, which takes place in case a deposited atom lands on another adatom, no additional energy barrier is assumed for it to diffuse down. On the contrary, atoms at the edge are not allowed to hop up a layer. Therefore, Ehrlich-Schwoebel mound formation is supressed, and the layer-by-layer growth is expected.

When an atom touches to other adatoms or islands, it makes a nearest neighbor bond to form a cluster. At very low temperatures, this cluster will never dissociate again. We assume that adatoms attach clusters irreversibly, and no cluster dissociation takes place. This means that dimers are stable, and the critical island size is expected to be $`i^{}=1`$. Nevertherless, if an edge atom is singly-bonded in the layer, i.e. if it makes only one bond with another adatom, it is possible to diffuse along the island periphery at moderate temperatures. There are two possible diffusion jumps for a singly-bonded atom; an edge diffusion along the straight step edge, and a corner diffusion or corner crossing round the outer corner. The rate of edge diffusion is $`k_e`$ and that of corner is $`k_c`$, as is depicted in Fig.1. Then, a singly-bonded atom can be classified by its possible motions in two directions. It is classified as (a) an $`ee`$ atom on a straight step if both movements are along the step edge, (b) an $`ec`$ atom at the ridge of a corner if one is the edge diffusion and the other is the corner, and (c) a $`cc`$ atom on a tip if both movements are across the corner, as shown in Fig.1. The edge diffusion constant along the straight step edge is $`D_e=k_ea^2`$, and we may define the corner diffusion constant as $`D_c=k_ca^2`$. Precise values of diffusion constants depend on the energy barriers and the temperature. For edge and corner diffusion there may exist extra energy barriers in addition to that for the surface diffusion, and diffusion constants are probably ordered as $`D_sD_eD_c0`$. If an edge atom is incoorporated into the kink site with more than two nearest neighbors, it is assumed to cease migration, for simplicity.

Simulation starts from a clean substrate surface of a size $`L\times L`$ without adatoms. Hereafter we take the lattice constant $`a=1`$. Therefore, the jump rates $`k`$’s and diffusion constants $`D`$’s are the same. During the simulation, at some stage, there are $`N_0`$ isolated adatoms, $`N_{ee}`$ edge atoms attached to a straight step edge, $`N_{ec}`$ atoms at the corner, and $`N_{cc}`$ atoms at the tip position (see Fig.1) on the surface. These adatoms are sorted out in corresponding lists. The transition probability from this configuration to the next one is given as follows. The rate of the deposition is $`P_d=FL^2`$, of the surface diffusion $`P_s=4N_0D_s`$, of the edge diffusion $`P_e=(2N_{ee}+N_{ec})D_e`$ and of the corner diffusion $`P_c=(N_{ec}+2N_{cc})D_c`$. The total probability of a state change in a unit of time is $`P_t=P_d+P_s+P_e+P_c`$. Thus, within a time interval $`dt=1/P_t`$, one of the events takes place; the deposition with a probability $`P_ddt`$, the surface diffusion with $`P_sdt`$, the edge diffusion with $`P_edt`$ and the corner diffusion with $`P_cdt`$. For each diffusion process, the atom to be moved is picked up randomly from the list. After each state change, adatom lists have to be adjusted, for example, by eliminating those adatoms with more than two neighboring bonds from the lists, or add an adatom that happens to have a single bond to the corresponding list. The time is increased by $`dt`$. Since the time unit is chosen such that the deposition flux $`F=1`$, the monolayer is covered at a time $`t=1`$, in principle. Because of the stochastic feature of the simulation algorithm, in practice, the actual time for the completion of a monolayer fluctuates around $`t=1.`$

## III Island Morphology

Simulations are performed for systems of sizes $`200\times 200`$, and $`1000\times 1000`$ under various combinations of diffusion constants, and some typical island morphology is shown in Fig.2 at a coverage of $`\theta =`$0.1 monolayer(ML).

From Fig.2(a) to (e) the edge and corner diffusion constants, $`D_e`$ and $`D_c`$, are equally increased from 0 to $`D_s`$, by keeping the surface diffusion constant $`D_s`$ fixed at $`D_s/F=10^9`$. Without peripheral diffusion, the island takes an irregular dendritic form or a fractal DLA-like shape, as shown in Fig.2(a). On increasing the edge and corner diffusion to $`D_e=D_c=10^6D_s`$, dendrite branchs thicken and they extend mainly in the $`<11>`$ direction, as shown in Fig.2(b). Further increase of $`D_e=D_c=10^5D_s`$ let the islands compact in a hopper shape with a corner instability due to Berg effect,berg38  as in Fig.2(c). Since the corner is extending in the region with high monomer concentration, the corner captures many adatoms and the edge diffusion transfers them to the island center. But the edge diffusion is still weak and insufficient, and the growth of the central part is too slow to catch up the fast corner growth. Thus, the island grows mainly in the $`<11>`$ direction, as in Fig.2(b) and (c). With sufficient edge and corner diffusion as $`D_e=D_c=10^4D_s`$, islands take square shape surrounded by rather straight {10} steps, as in Fig.2(d). So far island density seems independent on the island morphology. When all the diffusion constants are equal, as in Fig.2(e), islands are square and very large, but their density decreases drastically. On this density change we shall discuss later in §4. The sequence of simulations from Fig.2(b) to (d) almost corresponds to the parameter range studied by Bales and Chrzanbales+94  by changing the temperature. The morphology variation is of course quite analogous to what they have obtained.

We now study the effect of corner diffusion more precisely, in a sequence of pictures shown from Fig.2(e) to (i). There, the edge diffusion constant is kept as large as $`D_e=D_s`$, and only the corner diffusion constant $`D_c`$ is decreased from $`D_s`$ to 0. In all these cases, islands are not irregular. With a large corner diffusion constant, $`D_c10^4D_s`$, island are square with {10} steps, as shown in Fig.2(e) and (f). At $`D_c=10^5D_s`$ in Fig.2(g), they look round and rather isotropic, and for a still smaller $`D_c=10^6D_s`$ in Fig.2(h) some islands are in diamond shape and some others are in a cross shape with thick arms. In the extreme case of Fig.2(i) without the corner diffusion, $`D_c=0`$, islands consist four sharp needles growing in the $`<10>`$ direction. The same four-needle shape is observed previously by Zhong et al. zhong+01  in a submonolayer growth and by Caspersen et al.caspersen+01  in a multilayer study.

In a series of simulations Fig.2(i)-(l) the corner diffusion is completely suppressed ($`D_c=0`$), and the effect of the edge diffusion $`D_e`$ is exemplified. The series, in fact, continues back to Fig.2(a). With $`D_e=D_s`$ in Fig.2(i) needles in cross-shape islands are rather stable and sidebranches are rare. On decreasing the edge diffusion to $`D_e=10^4D_s`$ in Fig.2(j), needles have many sidebranches growing in the $`<10>`$ direction. Further decrement to $`D_e=10^5D_s`$ leads to tip splitting in Fig.2(k). With still smaller $`D_e=10^6D_s`$ clusters look like DLA aggregates. Comparison with Fig.2(b) and (l) reveals that the corner diffusion produces $`<11>`$ dendrite in Fig.2(b), whereas without corner diffusion aggregates have a preference in the $`<10>`$ direction in Fig.2(l). Also without corner diffusion, the dendrite tip is pointed, whereas with a corner diffusion, the tip becomes fat. At $`D_e=D_c=0`$, the DLA looks isotropic. Of course, it is known to be anisotropic asymptotically by reflecting the underlying square symmetry of the lattice, but the shape anisotropy appears only when the DLA grows very large. In the present simulation, islands coagulate to acomplish a layer growth before the asymptotic anisotropy appears.

From these variety of island morphology the edge and corner diffusions are found to play different roles on the island symmetry. The edge diffusion smears out shot noise introduced by the irreversible attachment, and permits the formation of regular shape. The corner diffusion varies the orientational preference. Without the corner diffusion, those atoms on the tip positions or on the narrow step edges cannot cross to the other sides of an island, and thus the supersaturation on the narrow {10} side increases. This leads to the enhancement in the growth velocity normal to the narrow {10} side, and due to the strongly anisotropic kinetic effect a needle tip is eventually formed into the $`<10>`$ direction. A corner diffusion allows the transfer of adatoms from one step edge to the other, and islands tend to be square. If the edge diffusion is weak, however, the adatom supply at the center of the step edge is meager, and the corner instability takes place. This leads to the dendrite growth in the diagonal $`<11>`$ direction, as in Fig.2(b) and (c).

## IV Island Density

In Fig.2 at a fixed surface diffusion constant $`D_s`$, we observe that the island density is almost independent of the edge and corer diffusion constants, $`D_e`$ and $`D_c`$. We said ”almost”, because in Fig.2(e) the island density is clearly different from the others. We now study the effect of peripheral diffusion on the island density.

The rate equation analysis predicts that the island density $`n`$ depends on the surface diffusion constant as in eqs.(1) and (2). venables73 ; stoyanov+81 ; villain+92a ; villain+92b  Experimentally, temperature dependence of the island density is used to estimate values of various energy barriers by assuming the critical island size $`i^{}`$. venables+84  Simulations where the temperature is varied also has to assume $`i^{}`$ to interpret their results. venables+84  To study the relation of island morphology and its density by fixing $`i^{}=1`$, we control directly on the rate of edge and corner diffusion, rather than tune many parameters as bond energies, energy barriers and the temperature.

The island density $`n`$ at the coverage $`\theta =`$0.1ML is plotted in Fig.3 as a function of the surface diffusion constant $`D_s/F`$ for various combinations of the edge and corner diffusion constants. The system sizes plotted are $`200^2`$ (open circles) and $`1000^2`$ (filled circles).

For the case without edge and corner diffusion, $`D_e=D_c=0`$, islands are fractal as shown in Fig.2(a), at least at low densities with a large $`D_s`$. The island density is shown in Fig.3(b). Data fits well with the scaling relation $`n=0.23(D_s/F)^{1/3}`$. The result is consistent to that obtained by Venables et al. venables+84  $`n=0.25(D_s/F)^{1/3}`$ from the rate equation analysis and that by Brune et al. brune+99  $`n=0.27(D_s/F)^{1.027/3}`$ for the DLA cluster in kMC simulation at the coverage 0.12ML. For the DLA, in fact, another scaling law $`n(D_s/F)^{2/(4+D_f)}`$ is proposed, where $`D_f=1.71`$ is the fractal dimension. villain+92a ; villain+92b ; tang93  From the present simulations of sizes up to $`1000^2`$, we cannot descriminate which of the two scaling laws is superior.

For the case with a large edge diffusion $`D_e=D_s`$ but without corner diffusion $`D_c=0`$, islands have cross shapes as in Fig.2(i). Even for this situation, the density $`n`$ is compatible with a scaling law $`n=0.25(D_s/F)^{1/3}`$, as is shown in Fig.3(a). The plot is shifted by a factor 10 to the right for the visual purpose. Otherwise, data overlaps almost completely with those for fractal islands, Fig.3(b).

With or without edge diffusion, the island density is found to follow the same scaling relation with the surface diffusion constant $`D_s/F`$, obtained by a simple rate equation theory.venables73 ; stoyanov+81  Inclusion of the corner diffusion, however, drastically changes the situation. By keeping $`D_e=D_c`$ but at different ratio to the surface diffusion as (c) $`D_e/D_s=D_c/D_s=10^4`$, (d) $`10^3`$, (e) $`10^2`$, (f) $`10^1`$, and (g) 1, the relation between the island density $`n`$ and the surface diffusion constant $`D_s/F`$ varies as shown in Fig.3(c)-(g). For (c) the island density at low $`D_s/F10^{10}`$ is higher than that for (b) with fractal islands by about 20 percent, $`n=0.28(D_s/F)^{1/3}`$. The density increase is due to the compactness of the islands, as shown in the surface configuration in Fig.4(b) for $`D_e/D_s=10^4`$, compared to the fractal islands in Fig.4(a) for $`D_c=0`$. There are more space for adatoms to nucleate new compact islands. As for the $`D_s`$ dependence, the island density follows the same scaling law with an exponent 1/3, unless $`D_s`$ takes the largest value $`D_s/F=10^{11}`$. There, the density becomes less than that of fractal islands, Fig.3(b).

As $`D_e`$ and $`D_c`$ increase further from Fig.3(d) to (g), the island density starts to deviate from the curve Fig.3(c) at smaller $`D_s/F`$. For the case Fig.3(d) with $`D_e/D_s=D_c/D_s=10^3`$, deviation takes place at about $`D_s/F=10^9`$, for (e) at about $`D_s/F=10^7`$. For (f) and (g) with larger $`D_e=D_c`$, the island density remains always less than that of fractal islands with the same $`D_s`$. This tendency of density decrease is obvious by plotting the surface configurations in Fig.4. There, the coverage $`\theta =0.1`$ML and the surface diffusion constant $`D_s/F=10^9`$ are the same, but the edge and corner diffusion constants, $`D_e/D_s=D_c/D_s`$, are varied. On increasing $`D_c`$ from Fig.4(b) to (f), the number of islands decreases and the island size increases.

From Fig.3(c) to (g), the island density $`n`$ seems to satisfy a new scaling relation to the surface diffusion constant $`D_s/F`$ as $`n(D_s/F)^\chi `$ with an exponent $`\chi 0.43`$, when the density $`n`$ is lower than that of fractal islands. In Fig.3 we depict fitting lines with an exponent $`\chi =3/7`$, which is expected from the critical mobile island size $`i^{}=3`$ and from eq.(2). By fitting the island density in the scaling form with an exponent 3/7 as

$$n=A(\frac{D_s}{F})^{3/7}$$
(3)
the coefficient $`A`$ depends on the peripheral diffusion constant $`D_c=D_e`$. By plotting $`A`$ versus $`D_c/D_s`$ in double logarithmic way as in Fig.5 the coefficient $`A`$ is found to be proportional to $`(D_c/D_s)^{0.195}`$.

From Fig.4, one further notices that as the corner diffusion increases, not only the total number of islands but also the island size distribution changes. bartelt+96 ; amar+95  This tendency is apparent by plotting island size distribution in Fig.6. Among the total $`N_i=nL^2`$ islands, there are $`N(s)`$ islands of a size $`s`$. Since the average island size $`s`$ changes as $`D_c`$ varies, the distribution of islands with a normalized size $`x=s/s`$ has the relevant concern. We define the probability $`P(x)\mathrm{\Delta }x`$ as the number of islands with sizes between $`x`$ and $`x+\mathrm{\Delta }x`$ devided by the total number of islands $`N_i`$ with $`\mathrm{\Delta }x=0.1`$. For each $`D_c`$, size distribution of 10 samples are averaged. It is normalized as $`_0^{\mathrm{}}𝑑xP(x)=_0^{\mathrm{}}𝑑xxP(x)=1`$, and corresponds to the scaled number density as $`N(s)=(\theta /s^2)P(s/s)`$ at a coverage $`\theta `$. bartelt+96 ; amar+95  At low $`D_c`$’s as $`D_c=0`$ and $`10^4D_s`$ in Fig.6(a), island sizes are widely spread with a broad peak such that there are very large islands with sizes about 2.5 times the average size. This is the typical size distribution for $`i^{}=1`$. bartelt+96 ; amar+95  On the contrary, at high $`D_c`$’s the island size distribution has a narrow peak around the average size, and extends less than twice the average size in the region of large sizes. This is the typical size distribution for $`i^{}>1`$. bartelt+96 ; amar+95 ; furman+97  When the clusters smaller then $`i^{}`$ are dissociable, the island size distribution is proposedamar+95  to have the form

$$P_i^{}(x)=C_i^{}x^i^{}\mathrm{exp}(i^{}a_i^{}x^{1/a_i^{}})$$
(4)
with

$$\frac{\mathrm{\Gamma }[(i^{}+2)a_i^{}]}{\mathrm{\Gamma }[(i^{}+1)a_i^{}]}=(i^{}a_i^{})^{a_i^{}},\text{and}C_i^{}=\frac{(i^{}a_i^{})^{(i^{}+1)a_i^{}}}{a_i^{}\mathrm{\Gamma }[(i^{}+1)a_i^{}]},$$
(5)
where $`\mathrm{\Gamma }`$ means the Gamma function. The probabilities corresponding for $`i^{}=1`$ and $`i^{}=3`$ are drawn in Fig.6(a) and (b) by continuous curves, respectively, and fit quite well with the simulation data, even though the dissociation is not allowed there. In fact, the similar size distributions with mobile clusters are obtained in the previous Monte Carlo simulation. furman+97  There, however, no peripheral diffusion is granted, and the grown islands had fractal structures, in constrast to our square shaped clusters.

## V Discussions and conclusions

We have seen that a large corner diffusion influences the island number density and the size distribution drastically. Without or weak corner diffusion, the island density varies in proportion to $`(D_s/F)^\chi `$ with an exponent $`\chi =`$1/3, irrespective of the island morphology. This means that the nuclei are formed at the initial stage of deposition, as soon as the two adatoms collide with each other. They are immobile and act as a center of nucleation. Therefore, the randomness of the two-adatom encounter is frozen in. After sufficient density of nuclei is formed, further deposited adatoms are incorporated into the preexisting clusters. They develope into various shapes, depending on the edge and corner diffusion constants. Since the nucleation takes place at random, the island sizes are distributed rather broad, reflecting the randomness in the island separation. In any case, this senario permits the island density to be independent of the morphology: The island density is determined long before the island shape appears.

With a large corner diffusion, on the other hand, inspections of the time evolution reveal that clusters of small sizes $`i3=i^{}`$ become mobile. Since atoms in a cluster smaller than the critical size $`i^{}=3`$ can move around each other via the corner diffusion, small clusters migrate around randomly on the substrate surface. During this cluster diffusion, they coagulate with isolated atoms and as a result the island density decreases. Of course, islands larger than the size $`i^{}=`$3 can migrate in some configurations, but when compact islands with doubly or more bonded atoms are formed, they cease migration on the substrate surface, and start to act as a nucleation center. During the further stage of cluster growth, the adsorbed atoms form small mobile clusters and are incorporated into islands in group.

This situation is treated by Villain et al. villain+92a  in terms of the rate equation, and the island density $`n`$ is expected to follow the scaling relation eq.(2). Since the trimer can migrate on the surface, $`i^{}=3`$, $`nF^{3/7}`$ is expected. This is the reason why we fitted the density with a scaling exponent $`\chi =3/7`$ in Fig.3. This explains the $`F`$-dependence of $`n`$, but as for the dependence on diffusion constants, we have to know $`D_2`$ and $`D_3`$ in terms of the present edge and corner diffusion constants, $`D_e`$ and $`D_c`$.

There is an exact calculation by Sanchez et al. sanchez+99  of diffusion constants of small clusters mediated by a peripheral diffusion. The result is simple for a dimer and a trimer on a square lattice as

$$D_2=\frac{D_c}{2},D_3=\frac{D_eD_c}{3(D_e+D_c)}.$$
(6)
For a dimer to diffuse, one of the adatoms which compose dimer has to cross round the corner. For a trimer, for example, in an extended form to move an edge atom has to cross the corner and to slide on the edge to another corner. Without edge diffusion, the center of mass of a trimer cannot move. In the case shown in Fig.3, $`D_e=D_c`$ and thus $`D_3=D_c/6`$. Inserting these exact results into the expression (2) obtained by rate equation, the island number density is expected to follow the relation

$$n\left(\frac{F^3}{D_sD_c^2}\right)^{1/7}$$
(7)
If this scaling holds, the coefficient $`A`$ in eq.(3) should be proportional to $`A(D_c/D_s)^{2/7}=(D_c/D_s)^{0.286}`$, instead the exponent -0.195 found in Fig.5. Therefore, the dependence of $`n`$ on cluster diffusion constant does not seem to be properly given by the rate equation.

To check the scaling behavior of the island density with the crossover from the monomer- to the trimer-limited form, let us study the scaled island density $`\stackrel{~}{n}=n(D_s/F)^{1/3}`$. It should be constant around 0.23 for small edge and corner diffusion constnat $`D_e=D_c`$. For large $`D_e=D_c`$, it should be proportional to $`(D_s^2F/D_c^3)^{2/21}`$, according to the result of the rate equation analysis (2). Therefore, we try the plot $`\stackrel{~}{n}`$ versus $`\stackrel{~}{F}_1=D_s^2F/D_c^3`$, but the data does not collapse on a single universal curve.

From our simulation result, the island density behaves as $`n=A(D_s/F)^{3/7}`$ with an $`D_c`$-dependent coefficient $`A=0.38(D_c/D_s)^{0.195}`$ for large $`D_c`$, as depicted in Fig.4. Since the exponent 0.195 is approximately equal to 4/21, the density is expected to be scaled by the variable $`x=D_c^2/FD_s`$ at a large $`D_c`$ more properly. The scaling plot is shown in Fig.7 for a wide range of combinations of diffusion constants as $`D_s/F=10^510^{11}`$ and $`D_e=D_c=10^{10}D_sD_s`$, and one observes a nice data collapse on a universal curve. It is fitted as

$$n\left(\frac{D_s}{F}\right)^{1/3}=\{\begin{array}{cc}0.23+0.05x^{0.1}\hfill & \text{for }x0\hfill \\ 0.38x^{2/21}\hfill & \text{for }x\mathrm{},\hfill \end{array}$$
(8)
as shown in Fig.7.

The effect of the edge and corner diffusion during the homoepitaxial nucleation growth on a singular surface is now summarized as follows. As for the island morphology, these peripheral diffusion makes the island regular and compact. With only the surface diffusion, the irreversibly nucleated islands have a fractal shape characteristic to DLA. With the edge diffusion, islands become thick and compact. The corner diffusion affects the shape anisotropy. One may typically say that without the corner diffusion dendritic needle pointing in the $`<10>`$ direction appear, whereas with it squares with {10} faces appear.

The density of islands is mainly determined before the shape appears. Without or with a weak corner diffusion, the critical island size is $`i^{}=1`$, and the density $`n`$ is related to the surface diffusion constant $`D_s`$ normalized by the deposition rate $`F`$ as $`n(D_s/F)^{1/3}`$, irrespective of island morphology. With a large corner diffusion constant, on the other hand, small clusters $`i3`$ can migrate randomly on the substrate surface, and another scaling relation $`nF^{3/7}/D_s^{5/21}D_c^{4/21}=(F/D_s)^{3/7}(D_s/D_c)^{4/21}`$ is found appropriate. The flux $`F`$-dependence agrees with the rate equation theory with $`i^{}=3`$, but the dependence on diffusion constants $`D_c`$ contradicts the theory. The reason of this discrepancy is not clear. But, since the rate equation is a mean-field type approximation which neglects spatial correlation, it may be fortuitous that $`F`$ dependence agrees in the simulation and the rate equation theory.

Recently, the mobility of large clusters is proposed to be very effective in selecting the island size and in arranging islands in order in heteroepitaxial growth.liu+01  Our model allows the motion of only small clusters less than the size $`i^{}=3`$ by the corner diffusion, and some size selection is observed. Due to the absence of any direct interaction between islands, however, the ordering in island positions cannot be achieved. Our model neither contains the size-limiting effect, and the coarsening takes place as the corner diffusion increases.

###### Acknowledgements.

This work is supported by Grant-in-Aid for Scientific Research from Japan Society of the Promotion of Science. The author is benefited from the inter-university cooperative research program of the Institute for Materials Research, Tohoku University.


format: text

text:
# Chemically active substitutional nitrogen impurity in carbon nanotubes

## Abstract

We investigate the nitrogen substitutional impurity in semiconducting zigzag and metallic armchair single–wall carbon nanotubes using *ab initio* density functional theory. At low concentrations (less than 1 at. %), the defect state in a semiconducting tube becomes spatially localized and develops a flat energy level in the band gap. Such a localized state makes the impurity site chemically and electronically active. We find that if two neighbouring tubes have their impurities facing one another, an inter–tube covalent bond forms. This finding opens an intriguing possibility for tunnel junctions, as well as the functionalization of suitably doped carbon nanotubes by selectively forming chemical bonds with ligands at the impurity site. If the intertube bond density is high enough, highly packed bundle of interlinked single–wall nanotubes can form.

Ever since their discovery1976nanotube ; iijima , carbon nanotubes (CNT) have been heralded as the new wonder–material of the future. Their remarkable mechanical and electronic properties destine them to play a major role in all kinds of nanotechnologies and molecular electronics. At least two major hurdles have to be overcome in order to fulfill this potential. First, manipulation of individual tubes is at best difficult today, which prevents mass production of devices. Second, the ability to fine–tune the various properties of the material to suit particular applications has to be achieved.

Carbon nanotubes exhibit a variety of electronic properties, for example, depending on their diameter and chirality, they vary from being metallic to semiconductingDresselhaus . The presence of defects and impurities that are electronically or chemically active can change these properties and thus have a significant bearing on a broad range of applications. In this paper we investigate the nitrogen substitutional impurity in single–wall carbon nanotubes (SWNT). In the case of semiconducting tubes, doping with electrons or holes is the principle route to making electronic devices. Alternatively, introducing new levels in the band gap with associated electronic states that are spatially localized, can create chemically active impurity sites.

Several groups have produced nanotubes containing nitrogen Yudasaka1997 ; Sen98 ; Terrones1999 ; Grobert2000 ; Han2000 ; Kurt2001 ; Czerw2001 ; Terrones2002 ; Trasobares2002 ; Wang2002 ; Lee2002  or nitrogen and boronSen98 ; Stephan1994 ; GolbergSWBCN . All of these were multi wall tubes, except for GolbergSWBCN , where nitrogen was always accompanied by boron. The measured atomic concentrations of the impurities were relatively high in all cases, above 1%, and all the CNx tubes were found to be metallic, as predicted by theoretical modelsTerrones2002 ; Kaun2002 ; Czerw2001 . Indeed, one other reason for the particular interest in CNx tubes was that they might make even initially semiconducting tubes metallic.

Various structural models have been proposed for the incorporation of nitrogen into the carbon network. Direct substitution was first studied by Yi and BernholcYi  using density functional theory at an atomic impurity concentration of about 1%, where they found an impurity state lying 0.27 eV below the bottom of the conduction band at the $`\mathrm{\Gamma }`$ point and a significant overlap of the states associated with the adjacent nitrogens (8.4 Å apart). Lammert *et al.*Lammert2001  considered disordered substitution of boron and nitrogen at a concentration of 3.5% and 1.4%, and observed both donor and acceptor levels in the gap. For an isolated nitrogen or boron impurity, they found the impurity state to be localized within about 10 Å for an (8,0) tube.

Another atomic arrangement for the experimentally observed nitrogen impurity was proposed by Terrones *et al.*Terrones1999 , where they suggest that at high concentrations (above about 10%) nitrogens are divalent, pyridine–like. Later publications from the same group refine this model and correlate the calculated density of states of metallic tubes with experimental observations Czerw2001 ; Terrones2002 .

To study the isolated substitutional nitrogen defect, the simulation cell has to be as long as possible. The existing constraints of computational power mean that the diameter of the tube should therefore be as small as possible. In all calculations we used the achiral zigzag (8,0) and armchair (5,5) nanotubes. Our simulation cell has periodic boundary conditions, the tubes are arranged parallel to the z axis and form a triangular lattice in the perpendicular x–y plane. This arrangement is experimentally seen in carbon nanotube ropesDresselhaus . We studied a range of nitrogen concentrations from 0.26% to 1% by introducting one impurity atom and changing the size of the simulation cell along the axis of the tube. The largest simulation cell contained 12 unit cells along the z axis for the zigzag tube and 14 unit cells for the armchair tube, respectively, with one impurity atom in it.

Throughout, we used *ab-initio* plane–wave pseudopotential density functional theorypayne  as implemented in the CASTEP codecastep . The generalized gradient approximation was used to account for exchange and correlation in the Perdew–Burke–Ernzerhof formgga-PBE . All calculations were spin-polarized and ultrasoft pseudopotentials were used for carbon and nitrogen with a cutoff energy of 300 eV. We assessed the adequacy of this cutoff by performing calculations at 350 eV for the smallest unit cells and verified that none of properties reported in this work changed perceptibly. For each concentration, the ionic positions were optimized by using one k–point at the $`\mathrm{\Gamma }`$ point for the largest simulation cells, and correspondingly more k–points for the smaller cells. After the geometry optimization has converged, the band structure was calculated using a finer mesh of k–points.

We find that the equilibrium position of the nitrogen ion is almost unchanged with respect to the corresponding C atom in the undoped nanotube, being moved by at most 0.01 Å, similarly to the case of azafullerene, C59NC59N ; gabor\_C59N . A slice of the calculated spin–density for both zigzag and armchair cases is shown in Fig. 1. The spin density is maximal on the N ion and the neighbouring C ions and shows oscillatory behaviour as it reaches a local maximum on every other ion in the network, again similarly to that observed in C59NC59N ; gabor\_C59N . The shape of the spin density isosurfaces is reminiscent of the $`\pi `$–orbitals in graphite with their axes of symmetry perpendicular to the sheet.

The localization of the extra electron associated with the isolated impurity shows a striking difference between the metallic and semiconducting cases, despite the fact that the band gap of the latter is quite small (0.6 eV according to the DFT calculations). In the armchair tube, the unpaired electron is almost completely delocalized along the tube, as can be seen from the plot of the envelope of the spin density in the lower panel of Fig. 1. However, for the zigzag tube, it is localized near the impurity site, but the localization region is fairly large, about 30 Å. This region contains about 90% of the unpaired spin density. The decaying envelope of the spin density is exponential.

The oscillations in the spin density can be understood as follows. Regarding the nitrogen atom as perturbation, the response of the system is characterized by the static susceptibility $`\chi (r)`$. The characteristic period in $`\chi (r)`$ corresponds to a maximum in the Fourier transform, $`\chi (q)`$, which happens for $`q`$ values that connect points in the Brillouin zone with the minimum energy difference. For the 2D graphene lattice, the energy gap is zero at the high symmetry K points, so for this, and similar systems, $`\chi (r)`$ will oscillate with a period corresponding to the momentum transfer between different K points, which is what is observed. The well known oscillations in the susceptibility in the RKKY formalismRKKY  are essentially a special case of the above considerations.

The decay of the envelope of the spin density is also analogous to the RKKY case. For metals in 1D, the result is $`\chi (r)1/r`$ aristov97 ; litvinov98 . For semiconductors, an analogous coupling to excited states takes place, and assuming a parabolic conduction band, it was shown rowland55  that $`\chi (r)\mathrm{exp}[(\alpha E_{\mathrm{gap}}^{1/2})r]`$, where $`\alpha =(2m_{\mathrm{eff}})^{1/2}/\mathrm{}`$. For the zigzag nanotube, using the calculated value of $`m_{\mathrm{eff}}/m_e=0.2`$ this crude estimate gives a decay constant of about 6 Å, close to the observed 7 Å.

Note that in all the graphene-based systems (C59N, nitrogen doped graphene and nanotubes), the spin density is positive almost everywhere. At first sight, this might seem to contradict the results from electron spin resonance hyperfine measurements, which indicate that the spin density oscillates symmetrically around zero. In fact, as was mentioned above, the spin density is very strongly dominated by the $`p`$ component, whereas the Fermi contact interaction in the hyperfine coupling is only sensitive to the $`s`$ component, so its observed symmetrical oscillation is consistent with the calculations.

A crucial difference between the armchair and the zigzag tubes is that the introduction of the nitrogen atom breaks the left–right mirror symmetry in the latter case. As two impurity atoms approach each other in an armchair tube, their interaction will show oscillatory behaviour depending on whether the spin density maxima induced by the respective nitrogens coincide or not. But in the zigzag tube the local maxima of the spin density occur on different sublattices on the two sides of the nitrogen atom. As two impurities approach, their respective spin density maxima never coincide, irrespective of whether the nitrogens occupy the same sublattice or not.

Figure 2 shows the band structure of the nitrogen-doped nanotubes along with the corresponding densities of states. The band structure of undoped tube (folded in such a way, that it fits into the shrinked Brillouin zone corresponding to the larger supercell) is shown for comparison. For the zigzag case, the extra band of the impurity falls somewhere near the bottom of the conduction band, and a complicated hybridization takes place with the existing unoccupied bands. The end result is that a new flat level appears about 0.2 eV below the bottom of the conduction band. In the armchair tube, the impurity introduces a peak at about 0.6 eV below the first van Hove singularity above the Fermi level, as well as a number of additional peaks, which are formed because of the hybridization and level repulsion taking place between the energy levels of impurity and pristine nanotube.

The localized electronic state of the semiconducting tube offers the intriguing possibility of being chemically active and forming covalent bonds with ligands or other tubes. Figure 3 shows two zigzag nanotubes each doped with nitrogen (one per four unit cells) in such a way that two carbon atoms with originally significant spin density on them face each other. As the charge density slice shows, a covalent bond is formed between these atoms, and thus between the tubes. The two carbon atoms pop out of the tangent plane of the tubes, and take up an $`sp^3`$ configuration, the bond length between them being 1.65 Å(note that a similar geometry was discussed in the context of carbon nitride materialssandre ). In order to obtain the geometry shown, the two tubes were pressed against each other, which is realistic when modelling a nanotube lying on top of another on a surfacetubetubetunneling . The further corroboration of the formation of a true covalent bond between the carbon atoms comes from the analysis of the charge density originating from different orbitals. The largest contribution to the total charge density between the newly bonded carbon atoms comes from the HOMO orbital (Fig. 3b), which is occupied by about 0.3 of an electron.

On the other hand, for the intertube bonds to be realized in a bundle of SWNTs, they have to be stable under zero external pressure. We found that this is possible if the number of intertube bonds is increased to one per two unit cells (8.4 Å), so that enough energy is gained in the new bonds to overcome the intertube repulsion. Here, the bonded configuration is a local minimum and the tube-tube distance is reduced to 2.5 Å, from 3.4 Åin the case of the unbonded bundle. This finding opens the possibility of producing highly packed crosslinked bundles of doped nanotubes; we obtained a crude estimate of 80 GPa for the bundle shear modulus (as a result of just one interlink for every 2 unit cells), which is substantially larger than 1-6 GPa, reported for ordinary bundlesForro .

Panel a) of Fig. 4 shows the $`p`$ projected local density of states (DOS) of various carbon atoms in different sections of the nanotube with one intertube bond per 2 unit cells, as compared to that of the unbonded system. The dashed curve shows the fine structure of the $`\pi ^{}`$ peak, which is almost completely absent for the $`sp^3`$ carbon atom (solid curve) that takes part in the intertube bond. As a first approximation, the $`p`$ projected total DOS in the conduction band, when suitably shifted, is proportional to the electron energy loss spectrum (EELS)egerton\_book ; rez99 , and shown on panel b). In the bonded system the $`\pi ^{}`$ peak is somewhat diminished, but more importantly, the peak marked with an arrow is absent. Since the proportion of $`sp^3`$ bonded carbon atoms is expected to be quite small, a local probe such as EELS would be most suitable for the detection of the intertube bond.

We performed *ab initio* calculations of zigzag and armchair carbon nanotubes substitutionally doped with nitrogen. In the case of the armchair nanotube, the impurity state is totally delocalized and the corresponding energy level falls inside the unoccupied bands of the pristine nanotube. In the case of the semiconducting zigzag nanotube, however, the isolated nitrogen impurity forms a flat energy level lying inside the band gap, 0.2 eV below the bottom of conduction band. It hybridizes with the $`\pi `$ orbitals to create a spatially localized state decaying exponentially with an overall extent of about 30 Å. This singly occupied state is chemically active, and we demonstrate the possibility of the formation of a covalent bond between two nanotubes with impurities that face one another. This opens the intriguing general possibility for chemical bonding between suitably doped carbon nanotubes or indeed with other ligands. Such tube–tube junctions could be interesting from a device perspective as the intertube bond will probably change the tunneling properties between the tubes considerably. If the density of intertube bonds is high enough, a bundle of SWNTs can become highly interlinked and packed, substantially enhancing its mechanical properties. The ligand–binding properties of the localized impurity site, of interest in chemical sensing applications, are currently under investigation.

###### Acknowledgements.

The authors would like to thank Peter Littlewood, Ben Simons, David Khmelnitskii, Malcolm Heggie, Bill Allison and Howard Hughes for valuable discussions. Computational work was carried out at the HPCF, University of Cambridge and supported by grant EC HPRN-CT-2000-00154.


text:
# Observing the dark matter density profile of isolated galaxies

## 1 Introduction

Measuring the distribution of mass around galaxies at large radii ($`100\text{kpc}`$) provides a critical test for cosmological models. By studying the mass distribution at these distances, we directly address one of the most intriguing questions of modern cosmology – the nature of dark matter. Specifically, we wish to determine the extent of galaxy dark matter halos and their density profile. While the inner parts of galaxy halos have density distributions that yield approximately flat rotation curves, our understanding of the profile at larger distances is much poorer. The issue is critical because the density profile that gives rise to a flat rotation curve ($`\rho r^2`$) is different from that predicted by cosmological models (NFW Navarro, Frenk, & White, 1997) at larger distances ($`\rho r^3`$).

The challenge in measuring the mass distribution at large radii arises from the difficulty in finding a visible tracer to probe the mass. Historically, neutral hydrogen (HI) has been used to study the outer parts of rotation curves (e.g., Roberts & Rots, 1973; Bosma, 1981; van Albada et al., 1985), since HI is detected well beyond the optical boundaries of spiral galaxies. The observations demonstrating that HI rotation velocities of field spirals are roughly constant at large galactocentric distances provides one of the primary pieces of evidence that galaxies are embedded in massive dark matter halos (e.g., Ostriker, Peebles, & Yahil, 1974; Faber & Gallagher, 1979; Sancisi & van Albada, 1987; Sofue & Rubin, 2001). However, HI emission is detected only out to 30–50 kpc. X-ray emission from the diffuse hot gas in isolated elliptical galaxies also indicates the presence of massive dark matter halos (e.g., Buote & Canizares, 1994; Romanowsky & Kochanek, 1998; Buote et al., 2002), but, just like HI in spiral galaxies, the X-ray emission does not extend far enough to distinguish isothermal from cosmological $`\rho r^3`$ profiles (Buote et al., 2002). Strong gravitational lensing places an important upper limit on the amount of dark matter but only in the inner (few tens of kpc) regions of galaxies (e.g., Keeton, Kochanek, & Falco, 1998; Keeton, 2001; Maller et al., 2000).

Weak gravitational lensing provides a more promising method of studying the outer parts of galaxies. Individual field galaxies produce a small distortion of background galaxies, and this distortion can be used to measure the galaxy–mass correlation function. Unfortunately, the existing data (e.g., Tyson et al., 1984; Brainerd, Blandford, & Smail, 1996; Hoeskstra, Yee, & Gladders, 2001; Fischer et al., 2000; Smith et al., 2001) do not distinguish between a singular isothermal sphere (SIS; Schneider & Rix, 1997) and a cosmological NFW profile: both models provide good fits to the galaxy–mass correlation function in the outer regions where the error bars are large. Moreover, the galaxy–mass correlation function can be affected by neighbors of the lens galaxy (e.g., Smith et al., 2001; Guzik & Seljak, 2002). Once the contamination from the neighboring galaxies is taken into account, as done by, for example, Fischer et al. (2000); Hoeskstra, Yee, & Gladders (2001), McKay et al. (2001) place only a lower limit on the size of the halo.

Velocities of satellites of galaxies provide another way to probe the mass distribution at large radii. These have been used to constrain the mass of the Milky Way (Zaritsky et al., 1989; Lin, Jones, & Klemola, 1995; Kochanek, 1996; Evans & Wilkinson, 2000) and M31 (Evans et al., 2000). However, because the number of detectable satellites around galaxies outside of the Local Group is small, one needs to study many galaxies to accumulate enough statistics to study the profiles of other galaxies. As a result, observational efforts to study the dynamics of satellites have been somewhat limited (Erickson, Gottesman, & Hunter, 1987; Zaritsky et al., 1993; Zaritsky & White, 1994; Zaritsky et al., 1997; McKay et al., 2002). Important early results were obtained by Zaritsky et al. (1993), Zaritsky & White (1994) and Zaritsky et al. (1997), who compiled and studied a sample of about 100 satellites of nearby isolated spiral galaxies with an average of 1–2 satellites per primary galaxy. It was found that the line-of-sight velocity dispersion of the satellites does not decline with the projected distance to the primary galaxy. This result has been generally considered a strong argument for the presence of dark matter at large distances ($``$200–400 kpc) around galaxies. Zaritsky & White (1994) also found that the satellite velocity dispersion does not correlate with the luminosity of the primary galaxy.

McKay et al. (2002) used the SDSS to select and study a much larger sample of 1225 satellites. The mean number of satellites they found around each host was only two. The new analysis confirmed that the velocity dispersion of satellites does not decline with distance from the primary, as observed by Zaritsky et al. (1993); Zaritsky & White (1994); Zaritsky et al. (1997). This implies that halos of isolated galaxies extend to distances of several hundred kpc. In contradiction with the earlier results, however, McKay et al. (2002) find that the velocity dispersion of satellites, $`\sigma `$, increases with the luminosity, $`L`$, of the primary galaxy as $`\sigma L^{0.5}`$ and claim that this is consistent with the predictions of the semi-analytical models of galaxy formation of Kauffmann et al. (1999) when the sample of primary galaxies in the models is defined as it is in the observational analysis.

Cosmological models made definite predictions for the mass profile of dark matter halos at large distances. Although SIS halo models with density $`\rho r^2`$ and $`\sigma =\mathrm{constant}`$ are often assumed, there is little justification for using SIS at these distances. The main motivations for using SIS are simplicity and extrapolation of flat rotation curves to much larger radii. However, during the last two decades of intensive numerical modeling of galaxy formation, not a single model has produced an SIS. Every cosmological model studied so far (and there have been plenty) has $`\rho r^3`$ at large radii. The slope does not depend on the mean density of matter: CDM with $`\mathrm{\Omega }_{\mathrm{matter}}=1`$ (Ghigna et al., 2000) and $`\mathrm{\Lambda }`$ CDM with $`\mathrm{\Omega }_{\mathrm{matter}}=0.3`$ (Klypin et al., 2001) have the same slope. It does not depend on the nature of the dark matter: warm dark matter (Avila-Reese et al., 2001), self-interacting dark matter (Colín et al., 2000), and cold dark matter all make the same prediction. The slope does not depend on the halo mass: halos ranging from cluster masses to dwarf masses all have $`\rho r^3`$. In other words, the declining velocity dispersion of satellites is not a test for the parameters of cosmological models, but rather of the hierarchical scenario itself. It is ironic that the same argument – the constant velocity dispersion of satellites – which just a few years ago provided one of the strongest cases for the existence of dark matter, is now an argument against dark matter. The only model that predicts a constant velocity dispersion is Modified Newtonian Dynamics (MOND).

As previous authors have noted (Zaritsky, 1992; Zaritsky & White, 1994), a comprehensive measurement of the satellite velocity dispersion profile around isolated galaxies should ideally combine two main elements: (1) a large number of primary galaxy and satellite (candidates) to estimate the satellite kinematics in separate bins of primary mass (or luminosity) and projected separation, and (2) a model-independent, robust, and testable rejection of interloper galaxies – dwarfs with large physical distances from the primary galaxies but small projected and velocity differences – that are not bound to the primary. The effect of the interlopers is to make the halo mass profile difficult to measure, and that may have led to systematic effects in the few observations that suggest a constant velocity dispersion at large radii. Understanding their effect is crucial for any interpretation of the halo mass distribution. Results from weak lensing can be affected by the presence of nearby projected neighbors, leading to similar systematic problems (see Smith et al., 2001; Guzik & Seljak, 2002).

In this paper we use the vast data base of SDSS to study the motion of satellites around a carefully selected sample of isolated nearby galaxies. We also develop an improved, model-independent approach for interloper rejection, which, together with state-of-the art cosmological $`N`$-body simulations of galaxies and their satellites, allows us to measure the dark matter halo profile.

In Section 2, we describe the observational data and the criteria for selecting our samples of primaries and satellites. The numerical simulations are briefly described in Section 3, where we also discuss the use of simulations to test prescriptions to remove interlopers. Our observational results are presented in Section 4, and our conclusions are given in Section 5.

## 2 Observational Data: Selection of primaries and satellites

The SDSS (York et al., 2000; Stoughton et al., 2002) is a survey that images up to $`10^4`$ $`deg^2`$ of the the northern Galactic cap in five bands $`ugriz`$ (Fukugita et. al., 1996; Pier et al., 2002; Hogg et. al., 2001; Smith et al., 2002) using a drift-scanning mosaic CCD camera (Gunn et al., 1998), down to a limiting magnitude of $`r^{}23`$ (Lupton et al., 2001). Approximately 900,000 galaxies down to $`r^{}17.7`$ will be targeted for spectroscopic follow-up using two fiber-fed spectrographs on the same telescope (Strauss et al., 2002; Blanton et al., 2002).

For the current analysis, we use data on galaxies for which spectra were obtained before August 2002. We use “survey quality” redshifts derived from these spectra by D. Schlegel. The recessional velocity errors are always less than 20 km s-1 . Our SDSS redshift sample consists of 254,073 galaxies distributed in several strips on the sky; the total area covered by our data is about 2,500 deg2. Photometric parameters were taken from the PHOTO measurement (Lupton et al., 2001) of the brightest object within 3 arcsec that was loaded into the SDSS collaboration database in August 2002 (but see also discussion below); magnitudes were corrected for foreground extinction using the (Schlegel, Finkbeiner, & Davis, 1998) values loaded in the SDSS database. The sky coverage of the resulting SDSS imaging and spectroscopic data is shown in Figure 1.

We supplemented the SDSS data with data from the RC3 catalog (de Vaucouleurs et al., 1991) for two reasons. First, as we are identifying isolated galaxies (see below), we want to avoid choosing galaxies that have bright neighbors outside the area covered by the SDSS at the time we made our sample; the RC3 provides a catalog of galaxies over the whole sky to make this possible. While it is true that the RC3 catalog is not complete, especially at larger distances, missing companions in the SDSS data are much more likely for the nearest galaxies for which the area on the sky corresponding to our isolation criterion is large. For these, the RC3 catalog should be nearly complete. Another, possibly more important, reason for using the RC3 is that the version of PHOTO that we were using did not provide accurate photometry for large and bright galaxies owing to problems with the deblending algorithm. As a result, we used RC3 $`B_T`$ magnitudes (or photographical magnitudes if $`B_T`$ was not available), corrected for foreground reddening, for all of the SDSS galaxies for which we found a match in the RC3 catalog (within 0.5 arcmin on the sky and 100 km s-1 in velocity); for all of the other galaxies (which are the smaller galaxies for which the SDSS photometry is accurate), we transformed the SDSS $`g`$ magnitude to $`B`$ magnitudes using the relation given by Yasuda et al. (2001).

Our main galaxy list was selected from the full redshift sample by taking all galaxies with $`500<cz<60,000\text{km s}\text{-1}`$ and apparent magnitude $`r^{}<17.7`$; the total number of selected galaxies is 206,352. Distances are derived assuming a smooth Hubble flow with $`h=0.7`$ ($`H_0=`$ 100 $`h^1`$ km s-1 Mpc-1); the heliocentric velocities were converted to the Local Group Standard of Rest before computing distances.

For the study of satellite dynamics, we define primary galaxies to be those with absolute blue magnitude brighter than $`M_B=19.0`$. We put no restriction on Hubble type. We use three different isolation and selection criteria to define three different samples of isolated primaries. To be considered as isolated, a galaxy must have no other galaxies within a magnitude difference $`\mathrm{\Delta }M`$, projected separation $`\mathrm{\Delta }R`$, and velocity separation $`\mathrm{\Delta }V`$. Satellites are defined as all objects within a projected distance $`\delta r`$ and velocity difference $`\delta v`$ and being at least $`\delta m`$ magnitudes fainter than the primary. Table 1 gives the parameters used to define each sample, and the number and median distance of primaries, and the number of satellites.

The first two samples have identical isolation criteria, but Sample 2 is defined for primaries out to a larger distance, hence only bright satellites can be seen for the more distant galaxies. The shallower sample has more satellites per primary, but the more distant sample has an overall larger number of total satellites. The third sample mimics conditions used by McKay et al. (2002), in which the isolation criterion is significantly relaxed. It has a very large search radius (almost 3 Mpc), but the luminosity of the neighbors can be quite large: half that of the luminosity of the primary. In other words, this prescription can pick a bright member of a group of galaxies and erroneously treat it as an “isolated” galaxy. As discussed below we obtain similar results for all of the samples, suggesting that our results are statistically significant and do not depend on the exact isolation criterion. Our main results are based on the first two samples.

Figure 2 shows distributions of primary and satellite brightnesses, brightness ratios, and number of satellites per primary for Sample 1. Considering only primaries for which there is at least one satellite, we find on average two satellites per primary in our first sample, which is similar to that found by Zaritsky et al. (1997) and McKay et al. (2002). A visual morphology classification of our primaries in Sample 1 indicates that about 30% of the galaxies are ellipticals and S0s (E–SO). This fraction of early-type galaxies is simiar as seen in the RC3 catalog and in the SDSS (Nakamura et al., 2002). A fraction of our primaries is close to the edges of the strips and therefore we do not cover their entire halos; for 84$`\%`$ of the primaries the whole volume is covered.

Given the limiting magnitude of the SDSS spectroscopic survey, we expect to find satellites down to $`M_B14`$ for a primary with $`cz=2000`$ km s-1 and down to $`M_B17`$ for a primary with $`cz=7100`$ km s-1, the mean of our primary recessional velocities in Sample 1. These limits would yield about four satellites per primary for our nearest primaries and 1–2 satellites for the further ones if we assume that all galaxies have satellite luminosity functions similar to that of the Local Group (see Pritchet, & van den Bergh, 1999; Grebel, 2000). Samples 2 and 3 include more distant primaries, and only brighter satellites for the more distant galaxies (see Table 1 for limiting magnitudes).

The main physical properties of satellites and primaries, such as their color distribution, spectral properties, and their subdivision into early- and late-type galaxies will be discussed in Vitvitska et al. (2003).

## 3 Numerical models

Cosmological simulations of galactic dark matter satellites in the hierarchical model of structure formation predict a remarkably large number of dark matter satellites orbiting around a Milky Way-size halo (Klypin et al., 1999; Moore et al., 1999). In many respects the simulations are quite realistic; we use them to develop a method for removing the effects of interlopers on the amplitude and shape of the r.m.s. velocity profile of the satellites. Note that the results of the simulations are not used in the reduction of the real astronomical data.

We use one of the simulations presented by Klypin et al. (2001). The simulation is done for the standard $`\mathrm{\Lambda }`$ CDM cosmological model with $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`h=0.7`$; the spectrum normalization is $`\sigma _8=0.9`$. The Adaptive-Refinement-Tree (ART) code (see Klypin et al., 2001, for details) is used to run the simulation. The simulation box was $`25h^1\text{ Mpc}`$. The formal force resolution was $`100h^1`$ pc; the mass resolution is $`1.2\times 10^6h^1M_{}`$.

The simulation has three galaxy-size halos. All halos are relaxed and have cusps. The outer parts of the density profiles of the halos are well approximated by an NFW profile. The distance between two of the halos is 600 kpc; we decided not to use this pair because the halos are not isolated. We selected the third halo, which does not have a massive companion within 3 Mpc radius. This halo has a virial mass of $`1.5\times 10^{12}h^1M_{}`$, a virial radius of 235 $`h^1`$ kpc, and a maximum circular velocity of 214 km s-1. There are 200 dark matter satellites inside a 350 $`h^1`$ kpc radius from the center of the halo. Most of the satellites are very small with circular velocity of 10–15 km s-1 and mass $`5\times 10^7h^1M_{}`$.

The simulations were done in such a way that only regions around the three large halos ($``$1/2 Mpc) have high resolution. The rest of the computational volume was simulated with significantly lower resolution. No halos outside the regions of high resolution are used in our analysis.

### 3.1 Satellites in cosmological simulations: removal of interlopers

In any sample of isolated bright galaxies and their satellites there is a fraction of galaxies that are probably not related to the host, but which are included in the sample because of projection effects. We call these objects interlopers. The interlopers tend to have large velocity differences and large projected separations. The fraction of interlopers in the samples of Zaritsky & White (1994) and Zaritsky et al. (1997) was $`15\%`$. We expect a similar level of contamination in our observational samples. If the interlopers are not removed, the r.m.s. velocity of satellites, and the inferred mass of the primaries, will be overestimated. We use the simulations described in the last section to test prescriptions for removing the interlopers.

To understand how interlopers might affect the results we introduce them into the simulation. We have to do this because the $`N`$-body simulation resolves only the region very close to the central halo; no interlopers are present in the simulation. To make a sample that realistically mimics the observational data we start by making an estimate of the number of expected interlopers in the observed sample. Using the SDSS field luminosity function (Blanton et al., 2001) we find that 0.43 galaxies between $`18<M_B<16`$ are expected inside a volume $`R<350h^1\text{ kpc}`$ and with a velocity dispersion of $`|\mathrm{\Delta }V|<500\text{km s}\text{-1}`$ at a distance of $``$5000 km s-1. This is the estimated number of interlopers per primary, if galaxies are randomly distributed in space. Neglecting the large-scale correlations of galaxies appears to be a reasonable assumption at this point, because our primaries are selected to avoid clusters and groups of galaxies. The latter are responsible for most galaxy clustering on a scale of a few megaparsecs. Thus, the clustering of galaxies (with the exception of true satellites) around our primaries is small and can be neglected. Using simulations this point should be quantified in more detail in the future.

In the observed samples we find on average two satellites per primary. Thus, we expect 0.215 interlopers for each true satellite. In our numerical simulation we have 200 true satellites. This means that we expect about 43 interlopers. We randomly place 43 interlopers into a cylinder of radius $`R=350h^1\text{ kpc}`$ and length $`1000\text{km s}\text{-1}`$.

Figure 3 shows the resulting velocity differences for both satellites and interlopers as a function of projected distance to the parent dark matter halo, as well as the r.m.s. velocities for the true satellites and for the combined sample of the satellites and interlopers. Since the interlopers are not physically associated with the primaries, their number depends only on the volume. Therefore, their effect on the r.m.s. velocity is stronger at larger distances. Neglecting interlopers clearly leads to an apparent increase in the r.m.s. velocity with distance, while the r.m.s. velocity of the true satellites actually decreases with distance.

To reliably estimate the satellite velocity dispersion, one needs to find a way to account for the interlopers. We start by characterizing the velocity distribution of true satellites. Motivated by Figure 3, we consider a Gaussian velocity distribution, but allow the width of the distribution to change with radius. Figure 4 shows the resulting distribution of normalized (to the r.m.s. velocity in different radial bins) velocities, and demonstrates that the distribution is very close to Gaussian. In reality, the distribution cannot be completely Gaussian because satellites with very high velocities will not be bound to the halo, and will escape. Using simple arguments based on the virial theorem, one naively expects that satellites with velocities twice the r.m.s. velocity should not be bound. Indeed, we find an indication of this: in Figure 4 the points are systematically below the Gaussian fit for $`\mathrm{\Delta }V>2\sigma `$. However, it is difficult to use the simple argument because we measure only one component of the velocity and because the escape velocity changes with distance. In any case, the Gaussian fit provides a reasonably accurate approximation for the distribution of the line-of-sight velocities of true satellites.

A simple prescription far removing the interlopers was used by McKay et al. (2002). The velocity difference distribution of the entire sample of neighbors (including objects at all projected distances) was fitted by a Gaussian + constant function. The idea is that the Gaussian describes the velocity distribution of satellites, and that the constant takes into account the distribution of interlopers; the width of the Gaussian component gives a direct measurement of the r.m.s. velocity of the satellites. Our simulations show that the velocities of the true satellites are indeed well approximated by the Gaussian distribution, but that the width of the distribution (i.e., the velocity dispersion) is a function of radius. We find that the McKay et al. (2002) prescription has two problems that result from making a single fit over all radii:

(1) The level of interloper contamination (the constant in the fit) does not depend on the distance to the primary. This is inaccurate because of a purely geometrical factor: the volume and thus, the number of interlopers, is larger for larger radii. Neglecting this dependence produces a bias: the r.m.s. velocity is overestimated at large radii and underestimated at small radii.

(2) There is no allowance for dependence of the velocity dispersion on projected distance. Because this dependence is in fact what we wish to study, we do not want to impose the assumption of a constant velocity dispersion on the prescription for the removal of interlopers.

In fact, if we apply the McKay et al. (2002) prescription to the simulations, we find that the recovered profile of the (interloper-removed) r.m.s. velocity is flat; the prescription fails to recover the true dispersion profile of the satellites. In this article, we apply an improved interloper prescription by making a separate “Gauss + constant” fit for separate radial bins. In this case the increase of the number of interlopers with the projected distance is taken into account (at least approximately). This also allows for a radially dependent velocity dispersion of true satellites. In practice we use $`100\text{kpc}`$ radial bins (70 $`h^1`$ kpc in numerical simulations) and 100 km s-1 velocity bins. We do not include satellites with projected distance less than $`20\text{kpc}`$. The binned data are used to minimize the sum over all bins of normalized deviations $`(V_{i,\mathrm{obs}}V_{i,\mathrm{model}})^2/N_i`$, where $`N_i`$ is the number of velocities in the $`i`$th bin, and $`V_{i,\mathrm{obs}}`$ and $`V_{i,\mathrm{model}}`$ are the observed and modeled velocities in the bin.

Figure 5 shows the accuracy of the recovery of the true r.m.s. velocity for the simulation; the short-dashed line shows the recovered velocity while the solid line shows the true velocity distribution of the satellites. The procedure works remarkably well. It recovers the value of the velocity dispersion and the number of interlopers. In the next section, we apply the same procedure to the observational samples.

The theoretical prediction for the velocity dispersion profile for an NFW profile with the same virial mass as that of the simulated dark matter halo is shown by a long-dashed line. In order to make this prediction for the r.m.s. line-of-sight velocity, we use the NFW profile and assume that the satellites are in equilibrium. We then use the Jeans equation to find the velocity dispersion. The velocity anisotropy $`\beta (r)1\sigma _{}^2/2\sigma _r^2`$, which is needed for the equation, was taken from cosmological simulations (Colín, Klypin, & Kravtsov, 2000; Vitvitska et al., 2002). (Here $`\sigma _{}`$ and $`\sigma _r`$ are the tangential and radial r.m.s. velocities). The velocity anisotropy changes from very small values $`\beta 0`$ close to the center of a halo to $`\beta 0.50.6`$ at the virial radius. We use this dependence of $`\beta `$ on radius to make analytical estimates of the projected velocity dispersion. Specifically, we integrate the r.m.s. velocities, which we obtain from the Jeans equation, along a line of sight. The integration is truncated at 1.5 virial radii. The velocity anisotropy enters these calculations twice: through the Jeans equation and during the integration along the line of sight. For a fixed mass profile $`M(r)`$ the effects of the velocity anisotropy are rather mild. The line-of-sight r.m.s. velocity in the realistic non-isotropic halo declines with distance slightly more rapidly as compared with the isotropic case ($`\beta =0`$).

## 4 Results

### 4.1 Radial dependence of satellite velocity dispersion

Figure 6 shows the relative line-of-sight velocities for a subset of Sample 2; only companions of galaxies with $`19.5<M_B<20.5`$ are shown. Figure 7 shows a comparable plot for Sample 3, but only for the brighter primaries, since the statistics are poor for the fainter primaries. The raw r.m.s. velocity is shown by the dashed curve. Both figures show that the r.m.s. velocity of the raw (uncorrected) data either does not change or slightly increases with the distance from the primary. A few years ago this would have been a strong argument in favor of dark matter in the outer parts of galaxies, but as discussed in the introduction, a flat r.m.s. velocity curve extending to 300 kpc strongly contradicts any existing dark matter model. Taken at face value, the raw r.m.s. velocities presented in Figures 6 and 7 contradict the dark matter scenario, and actually favor MOND.

However, if we remove interlopers using the Gauss + constant technique discussed in the previous section, as shown by the full lines with error bars in Figures 7 and 8, we find that the velocity dispersion of satellites declines with distance. We note that the results for Sample 1, which extends to fainter satellites, are similar to those shown here for the other two samples. Using Monte Carlo simulations, we reject a hypothesis that the r.m.s. velocity is constant at the 97% level corresponding approximately to the 3$`\sigma `$ level for Gaussian statistics. For the simulations we combine Samples 2 and 3, which provide independent measurements. We assume that deviations of the r.m.s. velocities are Gaussian in each of four bins for each sample. The average value of the r.m.s. velocity is defined by four random points. The probability of having $`\chi ^2`$ larger than the observed $`\chi ^2`$ and declining curves in two independent realizations is $`2.9\times 10^2`$, based on 10,000 realizations. Most of the constraints come from Sample 2, where the first and the last points deviate from the average by 1.5$`\sigma `$ and 1.8$`\sigma `$ correspondingly. Sample 3 alone is consistent with the constant r.m.s. velocity. Still, even in this case the velocity declines. When we combine both samples, the constant r.m.s. velocity becomes very improbable.

The number of estimated interlopers is quite modest. For Sample 1, which has 132 satellites for primaries in the range $`19.5<`$MB$`<20.5`$, we estimate that 23 are interlopers: 17% of the sample. The numbers are similar for the other samples: 19% in Sample 2 and 20% in sample 3. As expected from geometrical arguments, the number of interlopers increases with the radius. For example, in Sample 3 the number of estimated interlopers is 16.7 for radii in the range 100–200 kpc, and 54 for 200–360 kpc, which is roughly consistent with a constant number of interlopers per unit volume.

The observed decline in velocity dispersion with radius closely matches that predicted by an NFW dark matter density profile. Figures 7 and 8 give the predicted NFW velocities for halos of mass $`6\times 10^{12}\text{M}\text{}`$ (for primaries with $`20.5<M_B<21.5`$) and $`M_{\mathrm{vir}}=1.5\times 10^{12}\text{M}\text{}`$ (for primaries with $`19.5<M_B<20.5`$).

If we remove interlopers using the method of McKay et al. (2002), the mean velocity dispersion for each magnitude bin is smaller than the raw value, but we do not get a declining velocity curve. This explains the difference between our conclusions and those of McKay et al. (2002), who do not find a radial dependence on velocity dispersion; as discussed above, we feel that our procedure of allowing for the radial dependence of interloper contamination is better motivated physically.

### 4.2 Relation of satellite velocity dispersion with host luminosity

A comparison of Figures 7 and 8 clearly indicates that the r.m.s. velocity increases with the luminosity of the primary. This contradicts the earlier results of Zaritsky et al. (1997). Because the NFW profile provides an accurate fit to the observed velocity dispersion profile, we can use it to estimate the virial mass of each group of galaxies. For galaxies with $`19.5<M_B<20.5`$ we find $`M_{\mathrm{vir}}=1.5\times 10^{12}\text{M}\text{}`$. For the brighter magnitudes, $`20.5<M_B<21.5`$, the virial mass is $`M_{\mathrm{vir}}=6\times 10^{12}\text{M}\text{}`$. These values give us virial mass-to-light ratios $`M/L=100`$ and $`M/L=150`$ for the dimmer and brighter galaxies correspondingly. Thus, the $`M/L`$ ratio increases with the luminosity roughly as $`M/LL^{0.5}`$. We believe that this change in $`M/L`$ with luminosity is real and does not reflect a change in the stellar populations for both magnitude bins. A reason for this is that we did not see a dependence of the fraction of E–S0 with luminosity for the primaries in Sample 1, at least from $``$19 to $``$22 blue absolute magnitudes. A similar increase of $`M/L`$ with luminosity ($`M/LL^{0.4\pm 0.2}`$) has been seen from weak lensing studies by Guzik & Seljak (2002). See, however, McKay et al. (2001) who, also from weak lensing, found no dependence of the $`M/L`$ on luminosity.

Instead of trying to estimate masses of the galaxies (which is bound to be model-dependent), we can study the more direct dependence of the r.m.s. velocity on the luminosity of the primary galaxy. This can be viewed as an analog of the Tully–Fisher relation. Figure 9 presents the dependence of the satellite r.m.s. velocity, computed within 120 kpc as a function of galaxy luminosity. We calculate the r.m.s. velocity relatively close to the primary galaxy because within this distance the velocity dispersion depends only weakly on the distance to the primary for all primaries considered here. Hence, we can meaningfully compare results from galaxies that might have different sizes. The estimated value of the velocity dispersion for an $`L_{}`$ galaxy from weak lensing studies (Hoeskstra, Yee, & Gladders, 2001; McKay et al., 2001) is also shown; the results of lensing studies are very close to those derived from the SDSS satellite data.

The results are consistent with a slope of $`\sigma L^{0.3}`$. We can write this dependence using a form analogous to the Tully–Fisher relation:

$$M_B=0.028.7\mathrm{log}(2\sigma ).$$
(1)
This relation is shown in Figure 9. The slope of this relation is the same as for spiral galaxies (Verheijen, 2001), which is shown as a long-dashed line in the figure. The zero point for the spiral relation was shifted down by factor of 1.6 to match the SDSS data points; this factor is close to the naively expected $`\sqrt{3}`$ correction from rotational velocity to 1D line-of-sight velocity. The remaining difference of 8% is difficult to assign to any particular effect; there are several effects that could produce the small difference (e.g., projection effects, and non-isotropic orbits).

However, when we compute the r.m.s. velocity within a larger radius of 350 kpc, we find $`\sigma L^{0.5}`$. The reason for the steeper slope is a combination of the very large radius and declining $`\sigma (r)`$, which will occur if the density profiles of halos at different mass are not homologous, e.g., if the concentration of the halo decreases with increasing mass. This is consistent with theoretical prediction.

It is rather straightforward to find the $`\sigma `$–virial mass relation for the NFW halos. Using the model described in Section 3.1 (equilibrium NFW with slightly anisotropic velocities), we find the line-of-sight r.m.s. velocity averaged for radii 20–100 kpc and the line-of-sight velocity for projected distance 350 kpc. The results are shown in Figure 10. For the range of masses considered here, the $`\sigma `$$`M`$ relations are power laws. Just as the observed slope of the $`\sigma `$$`L`$ relation, the slope of the $`\sigma `$$`M_{\mathrm{vir}}`$ relation increases with the radius within which the r.m.s. velocity is calculated. For the 350 kpc radius the slope 0.5 is the same as the slope of the $`\sigma `$$`L`$ relation for the same radius. This would imply a constant mass-to-light ratio. This result, being formally correct, is a bit misleading. Measured in terms of a more physically motivated quantity, the virial radius, the $`M/L`$ ratio increases with luminosity.

Comparison with theoretical predictions of the luminosities is more complicated and will be deferred to another paper. We use only the virial $`M/L`$ for two magnitudes $`M_B=20`$ and $`21`$. In a recent paper, Yang et al. (2002) ask which model parameters are required for the cosmological models to explain the observed luminosity function and the Tully–Fisher relation. The main uncertainty is the $`M/L`$ ratio. Yang et al. (2002) predict an $`M/L(L)`$ that is needed to account for observations. We compare our two $`M/L`$ ratios with the Yang et al. (2002) predictions and find them remarkably close. However, this comparison actually required numerous corrections, including a correction for overdensity of 180 to the virial overdensity, correction for galactic absorption, correction for different bands, and scaling with the Hubble constant. The corrections were provided by Rachel Somerville. Preliminary comparison with semi-analytical models (Somerville, private communication) indicate a problem: observed galaxies are $``$0.75 mag too bright as compared with theoretical predictions.

## 5 Conclusions

We construct samples of isolated galaxies and associated satellites using data from the SDSS. We detect about 3000 satellites with absolute blue magnitudes going down to $`M_B=14`$. Analysis of the satellite velocities clearly indicates that galaxies are embedded in large halos extending to distances up to $`350\text{kpc}`$. Our measured velocity dispersion for $`L_{}`$ galaxies is in agreement with measurements from weak lensing studies.

Our analysis of galaxies selected from the SDSS database clearly indicates that the line-of-sight r.m.s. velocity of satellites declines with the distance to the primary galaxy. Similar results were found for each of our three samples, suggesting that it is robust to statistical effects and the exact definition of an isolated galaxy. This decline agrees remarkably well with $`\rho r^3`$ predicted by all cosmological models. Both the isothermal and MOND profiles contradict the observational results. Observations of weak lensing from the SDSS combined with the Tully–Fisher and Fundamental Plane relations (Seljak, 2002) also finds that the galactic halo profiles has to be steeper than isothermal at large radii.

Interlopers play an important role. We believe that the decline in satellite velocity was not seen previously (Zaritsky & White, 1994; Zaritsky et al., 1997; McKay et al., 2002), mainly because of either lack of statistics or insufficient removal of interlopers. We use numerical simulations to test prescriptions for correcting the effects of interlopers and isolation criteria. Both effects – interlopers and isolation criteria – have a tendency to overestimate of the mass of the dark matter halo and imply a velocity profile that is too flat.

We find that the satellite velocity dispersion inside a projected 100 kpc radius and the absolute blue magnitude of the galaxy are related as $`M=0.028.7\mathrm{log}(2\sigma )`$ ($`\sigma L^{0.3}`$). This dependence is in agreement with the slope and the zero point of the Tully–Fisher relation for spiral disks. If we calculate the r.m.s. velocity within a radius of $`350\text{kpc}`$, the slope of the velocity–luminosity relation is visibly steeper: $`\sigma L^{0.5}`$. The likely explanation for the difference in the slopes is an interplay between the shape of the velocity dispersion and the virial radius of galaxies.

The mass-to-light ratio increases with luminosity as $`M/LL^{0.5}`$. For $`M_B=20`$ (-21) we find that the virial $`M/L=100`$ (150).

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the US Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the participating institutions. The participating institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. FP acknowledges the hospitality of the NMSU Astronomy Department, where part of this work was done, and to Anatoly Klypin for financial support during this visit. MV acknowledges the hospitality and financial support of the MPIA, where the work was started. MV was supported by the Los Alamos National Laboratory. JH acknowledges the hospitality of the Instituto de Astrofísica de Andalucía, where some of this work was completed.