text:
# On resonance line profiles predicted by radiation driven disk wind models.

## 1 Introduction

It appears that mass accretion onto compact objects through accretion disks is often accompanied by mass outflow from these disks. For example, disk outflows are observed in active galactic nuclei; non-magnetic cataclysmic variables (CVs); and young stellar objects. In the case of CVs, key evidence for outflows comes from P-Cygni profiles of strong UV lines such as C IV$`\lambda `$1549. However, the evidence for the outflows is not limited just to the strong UV lines (e.g., Long & Knigge 2002 and references therein). Understanding the outflows in CVs is important because they have been the best observed outflows from compact objects and promise to provide us with insights into all disk outflows. Until recently the interpretation of data was limited to fitting observed profiles to synthetic profiles calculated from kinematic models (e.g., Mauche & Raymond 1987; Drew 1987; Shlosman & Vitello 1993; Knigge, Woods & Drew 1995; Long & Knigge 2002).

The most plausible mechanism for driving the CV outflows is radiation pressure on spectral lines. Recently, numerical hydrodynamic models of 2.5-D, time-dependent radiation driven disk winds have been constructed for application to CV disk winds (Pereyra, Kallman & Blondin 1997; Proga, Stone & Drew 1998, hereafter PSD 98; Proga, Stone & Drew 1999, hereafter PSD 99). These initial studies confirm viability of the radiation driving mechanism. However, the same studies reveal some difficulties in obtaining line profiles matching the observed profiles (Proga et al. 2002, Paper 1 hereafter). The latter is likely true even when we allow for the magnetic field effects (Proga 2003).

We focus here on re-examining pure radiation-driven disk wind models for CVs. In particular, we examine the size of the region where the line is formed and show that what is required to improve the model predictions is to capture the entire line formation region by increasing the radial range of the models. In Section 2, we summarize the key elements of our calculations. The results are presented in Section 3. In Section 4, we discuss key aspects of the synthetic line profiles and make comparisons with recent observations.

## 2 Method

For this paper, we recalculated dynamical models presented by PSD 99 and calculated a few new models. We refer a reader to PSD 99 for details on wind models and to Paper 1 for details on line-profile calculations. The computations of line profiles are exactly as in Paper 1. We made only two minor modifications to the models: (i) we increased the outer radius of the computational domain, $`r_o`$, from 10 to 100 white dwarf radii and (ii) as in PSD 99, we computed models with the radiation force using the intensity of the radiation integrated over all wavelengths, $`\lambda `$, but also using the intensity of the radiation integrated over the UV band only (i.e., $`200\AA \lambda 3200\AA `$). We call the former $`I`$ and the latter $`I_{\mathrm{UV}}`$. We found that the second modification can change the wind solution (e.g., there is no outflow originating from large radii on the disk) but typically this is not much reflected in the predicted line profiles. On the other hand, the first modification does not change the wind solution but has a significant effect on the predicted line profiles.

## 3 Results of Profile Calculations

Our basic set of line profiles is based on four hydrodynamical models computed on the larger computational domain with the model parameters as in PSD 99’s models A, B, C, and D (see Table 1). We call our models recalculated with the intensity $`I`$ on the larger grid A1, B1, C1, D1, respectively, while models recalculated with the intensity $`I_{\mathrm{UV}}`$ on the larger grid are labeled A2, B2, C2, D2. To examine the effect of viewing angle on the line profiles, for each disk wind model we compute line profiles for five inclination angles: $`i=5^{},30^{},55^{},70^{}`$, and $`85^{}`$. Examining effects of inclination angle is important because of the strong inclination angle dependence observed in CV spectra. Figure 1 shows some of the model line profiles obtained after scaling to unit continuum level.

To illustrate the key differences between profiles from Paper 1 and our profiles, we present results from the B models and the C models. We mention however, that the changes in the line profiles for the A models and the D models are qualitatively similar to the changes for the B models and the C models.

Table 1 summarizes the main input parameters of the hydrodynamical disk wind models including the mass accretion rate, $`\dot{M}_a`$, total stellar luminosity relative to the disk luminosity, $`x`$, and resultant total system luminosity, $`L_{tot}=L_D+L_{}=(1+x)L_D`$. Additionally, the table lists, the gross properties of the disk winds including the mass-loss rate, $`\dot{M}_w`$, characteristic velocity of the fast stream at 10$`r_{}`$, $`v_r(10r_{})`$ for model A, B, C, and D and their counterparts for the models based on the larger grid, $`v_r(100r_{})`$, and flow opening angle, $`\omega `$.

First, we present our results for a wind for which $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ and the white dwarf (WD) is assumed to be dark $`x=0`$ (model B, B1, and B2). We note that the B models have similar properties, except for $`\dot{M}_w`$ which is lower for B2 than for models B and B1 (see Table 1). The lower $`\dot{M}_w`$ in model B2 is due to the fact that the radiation force in model B2 is lower than for models B and B1 as it is computed using $`I_{\mathrm{UV}}`$, which is lower than $`I`$. The absence of a wind model A2 is caused by the same effect, at a more extreme level. The A models are for a relatively low luminosity and a reduction of the driving radiation flux resulting from taking into account only the UV photons makes the total radiation force too weak to overcome the gravity and drive a wind.

The five top panels of Figure 1 (Figs. 1a-1e) compare the profiles for model B as a function of $`i`$ from Paper 1 (thin solid lines) with the profiles predicted by models B1 and B2 (dashed and thick solid lines, respectively). The line profiles from Paper 1 for $`i=55^{}`$ and $`70^{}`$ show two maximum absorptions almost equally shifted from the line center to the blue and red. This effect is due to the fact that the rotation dominates over expansion in shaping the line profile for $`r<10r_{}`$. There is little absorption near the line center. The most obvious change in the profiles caused by taking into account the wind at larger radii is the increase of the absorption near the line center. Our detailed analysis of the line-forming region shows that this central absorption is due to the wind originating at $`3<r<10r_{}`$ which has optically thick resonance surfaces at $`r<30r_{}`$ to the photons close to the line center. Thus, the ‘third’ component of the absorption due to the wind at large radii changes qualitatively the shape of the total line profile for $`i=55^{}`$ and $`70^{}`$ and resolves the problem of the two maximum absorptions. For $`i<30^{}`$, the outer wind significantly contributes the line absorption at all wavelengths.

We note that the red and blue edges of the lines do not change when the wind at larger radii is taken into account (compare thick solid lines with dashed and thin solid lines in Fig. 1). This reflects the fact that the old and new wind solutions are very similar within $`10r_{}`$ and that the old models capture well the acceleration zone of the fastest part of the wind (e.g., PSD 98). Note that the fastest wind is responsible for the position of the blue edge of the absorption at the lowest inclination.

The line profiles based on our new models have a similar dependence on the inclination angle as the line profiles from Paper 1. Namely, the absorption weakens whereas the scattered emission strengthens as inclination angle increases. Additionally, the red edge is at the line center for $`i=5^{}`$ and moves to the red as $`i`$ increases for $`i<85^{}`$.

The five bottom panels of Figure 1 (Figs. 1f-1j) compare the profiles for model C as a function of $`i`$ (thin solid lines) with the profiles predicted by models C1 and C2 (dashed and thick solid lines, respectively). Model C is an example of a strong steady state wind where $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ and WD radiates at the same rate as the disk, $`x=1`$. In analogy with the B models, the C models have similar properties except for $`\dot{M}_w`$, which is lower for C2 than for model C and C1 (see Table 1).

The line profiles for model C (fig.1f-1j) are different from those for model B for $`i>30^{}`$. For model B, the line develops a strong, more or less symmetric, double-peaked emission when seen nearly edge-on whereas for model C, the line develops a typical P Cygni profile with blueshifted absorption and redshifted emission. Additionally, the line is significantly broader for model C than for model B. The rotation of the wind does not dominate over expansion in shaping the line. In particular, the shape and strength of the blueshifted absorption is due to expansion. For $`i<55^{}`$, there is a rotationally broadened redshifted absorption. The redshifted absorption is weaker and narrower than the expansion-dominated blueshifted absorption.

The contribution from the outer wind causes similar changes to the line profiles for model C as it does to the profiles for model B. The line profiles for model C from Paper 1 show a deep narrow absorption at $`v1000\mathrm{km}\mathrm{s}^1`$ for $`i=70^{}`$ (Fig. 1i) and at $`v200\mathrm{km}\mathrm{s}^1`$ for $`i=85^0`$ (Fig. 1j). The central absorption due to the outer wind broadens and strengthens the absorption for $`i=70^{}`$ and eliminates the central emission which has a peak at $`v500\mathrm{km}\mathrm{s}^1`$ in the line profile from Paper 1. Additionally, the outer wind strengthens the red-shifted emission for $`v<1500\mathrm{km}\mathrm{s}^1`$. As a result, for $`i=70^{}`$, the line profiles predicted by models C1 and C2 does not exhibit the double-humped structure seen in the line profile predicted by model C. For $`i=85^{}`$, the line profiles predicted by models C1 and C2 develop a double-peaked emission without the deep narrow absorption at $`v1000\mathrm{km}\mathrm{s}^1`$ which is present in the profile from Paper 1. This narrow absorption disappeared in the profiles for model C1 and C2, mainly because the blue-shifted absorption is weaker in our new model than in the old model.

For $`i<70^{}`$, the line profiles for model B and C are similar in three respects: (i) the scattered emission is weak, (ii) the contribution from the outer wind is negligible in the red part but it is significant in the blue part of the profiles, and (iii) the contribution from the outer wind does not change the position of the red and blue edges of the lines.

In analogy to models B and C, at $`i=85^{}`$, the scattered emission is seen as weaker by the observer for models C1 and C2 than for models B1 and B2 because the background continuum emission is stronger in the C models. At this high inclination, the background continuum of the C model is dominated by the WD contribution (see Table 1 in Paper 1).

## 4 Discussion and comparison with observations

Our main result is that the inclusion of a disk wind at larger radii changes qualitatively and quantitatively the line profiles predicted by the pure radiation-driven disk wind model. The models computed on a small grid – such as those in PSD 99, where the outer radius equals 10 $`r_{}`$ – are sufficient to calculate the gross properties of the disk wind. For example, PSD 99’s simulations are sufficient to calculate the wind mass-loss rate and the fast part of the wind, which are both associated with the outflow from the innermost disk. However, the previous simulations do not capture the entire region where lines are formed. As a result, they underpredict the line absorption and to a lesser extent the scattered emission. The previous simulations predicted a double-humped structure near the line center for intermediate inclinations. This structure is due to a nonnegligible red-shifted absorption that is formed in the slow wind where the rotational velocity dominates over expansion velocity. We showed here that by taking into account the downstream part of the same slow wind we were able to increase significantly the central absorption. As a result, the double-humped structure is reshaped to a more typical broad trough. We emphasize that all improvements in the shape as well as the strength of the absorption were achieved without changing the gross properties of the wind. In particular, our new models do not predict a higher mass-loss rate than the previous models. The changes in the line profiles are mainly caused by the fact that the ratio between the rotational and poloidal velocity decreases downstream.

In Paper 1, we showed that the double minimum apparent near line center can be replaced by a more plausible single, blueshifted minimum when we make the slow wind transparent. We concluded then, that the wind-formed line profiles seen at ultraviolet wavelengths cannot originate in a flow where rotation and poloidal expansion are comparable. The UV lines must trace gas that expands substantially faster than it rotates. Our new results support this reasoning.

Our new and old line profiles have many features that are qualitatively consistent with the features of the observed line profiles. These include the classical P-Cygni shape for a range of inclinations, the location of the maximum depth of the absorption component at velocities less than the terminal velocity (a requirement discussed in Drew 1987), and the transition from net absorption to net emission with increasing inclination. However, our new line profiles are also qualitatively consistent with the observations, in particular, the depth and width of the absorption. The main discrepancy between the predicted line profiles and the observed ones is in the line emission. Specifically, the model cannot produce the redshifted emission as strong as that seen, for example, in the Civ profile of many systems with intermediate inclinations (see below). However, this shortcoming is not a great surprise – this has been a problem for a while (see e.g. Mauche, Lee & Kallman 1997; Ko et al. 1998).

We are undertaking a systematic comparison between predicted line profiles and observations for many systems. Such a study is needed to support our finding that the pure radiation-driven disk wind model can work for CV winds in the sense that it can reproduce the observed line profiles for model parameters (e.g., system luminosities) suitable to CVs. Our preliminary results from limited survey of dynamical models and their predictions are promising. To illustrate this, we compare our new profiles with spectrally well-resolved UV observations for the brightest nova-like variable, IX Vel (Hartley ; in Paper 1, we made a model comparison with the same data). In Figure 2, we show a comparison between profiles derived from the model with $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$, $`x=0.25`$ (model E2) and observations of the Civ 1549Å and Siiv 1403Å transitions. To show how much line emission is required we plot only the absorption component of the line synthesized for $`i=50^{}`$ and $`60^{}`$. The observed $`i`$ for this system is $`60^{}`$ (Beuermann & Thomas 1990). The choice for the relative luminosity of WD (i.e., $`x=0.25`$) was motivated by the fact that the line profiles predicted by model B2 with $`x=0`$ and model C2 with $`x=1`$ bracket the observed profiles.

Figure 2 clearly shows that the model profiles well reproduce the blue-shifted absorption. Therefore, we conclude that we have significantly narrowed the gap between the kinematics of the radiation-driven models and reality. The gap is only narrowed but not bridged yet because the mass fluxes required to match the observed spectra at least of IX Vel are somewhat higher than those observed. In particular, it is believed that the luminosity of the system requires a mass accretion rate of at most $`10^8\mathrm{M}_{}\mathrm{yr}^1`$, whereas our line profiles require $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ and $`x=0.25`$, yielding a system luminosity higher than the observed one by a factor of $`4`$. However, our main point is that this discrepancy is much smaller than it used to be (i.e., it was more than 1 order of magnitude) and we believe that it can be reduced still further. We have computed many models, changing various model parameters such as $`\dot{M}_a`$ and the parameters of the force multiplier, $`\alpha `$ and $`M_{max}`$. In general, we find that there is a degeneracy in the model parameters as far as line profiles are concerned (PSD 98 found an analogous degeneracy for the wind properties). For example, models with slightly different parameters – such as $`\dot{M}_a`$ and $`\alpha `$ – produce similar line profiles for different $`i`$. Additionally, a model with the same parameters as model A2 but with $`\alpha =0.674`$ instead $`\alpha =0.6`$ predicts very similar line profiles to model B2 for the same $`i`$. Finally, we find that the product $`(1+x)L_DM_{max}`$, not its individual factors, appears to be a fundamental parameter determining the line profiles (i.e., their width and depth) for the parameter range applicable to CVs. This has an important implication for our models: to obtain a theoretical fit as good as shown in Fig. 2 for a fixed $`i`$, we need $`(1+x)L_DM_{max}1.3\times 10^5\mathrm{L}_{}`$ rather than specifically $`x=0.25`$, $`M_{max}=4400`$, and $`L_D=23.4\mathrm{L}_{}`$ as for model E2.

As in Paper 1, some aspects of our results depend on our assumptions about the disk and microphysics of the wind. Further work should include in particular, self-consistent calculations of the wind ionization state.

ACKNOWLEDGMENTS: We thank J.E. Drew and M.C. Begelman for their comments. We acknowledge support from NASA under LTSA grants NAG5-11736 and NAG5-12867. We also acknowledge support provided by NASA through grant AR-09532 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## REFERENCES

Beuermann, K., Thomas, H.-C., 1990, A&A, 230, 326

Drew, J.E. 1987, MNRAS, 224, 595

Hartley, L. E., Drew, J. E., Long, K. S., Knigge, C., & Proga, D. 2002, MNRAS, 332, 127

Knigge, C., Woods, J.A., Drew, J.E. 1995, MNRAS, 273, 225

Ko, Y., Lee, P. Y., Schlegel, E. M., & Kallman, T. R. 1996, ApJ, 457, 363

Long, K.S., & Knigge, C., 2002, ApJ, 579, 725

Mauche, C.W., Lee, Y.P., Kallman, T.R. 1997, ApJ, 477, 832

Mauche C.W., Raymond J.C. 1987, ApJ, 323, 690

Mauche C.W., & Raymond J.C. 2000, ApJ, 323, 690

Pereyra, N.A., Kallman, T.R., & Blondin, J.M. 1997, ApJ, 477, 368

Proga, D. 2003, ApJ, 585, 406

Proga, D., Kallman, T.R., Drew, J.E., & Hartley, L.E., 2002, ApJ, 572, 382 (Paper 1)

Proga, D., Stone J.M., & Drew J.E. 1998, MNRAS, 295, 595 (PSD 98)

Proga, D., Stone J.M., & Drew J.E. 1999, MNRAS, 310, 476 (PSD 99)

Shlosman I., & Vitello P.A.J. 1993, ApJ, 409, 372

Figure Captions

Figure 1 – Line profiles for hydrodynamical disk wind models as a function of inclination angle, $`i`$ (see top right corner of each panel for the value of $`i`$). The top panels, a-e, are for models with $`x=0`$ but with $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ (B models). The bottom panels, f-j, show results for models with $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ and $`x=1`$ (C models). The thin solid lines are the line profiles based on the models B and C from PSD 99 (these are the same profiles as in Paper 1, see figure 1 there). The dashed lines are line profiles based on our models B1 and C1 while the thick solid lines are profiles based on models B2 and C2. The zero velocity corresponding to the line center is indicated by the vertical line. Note the difference in the velocity and flux ranges in the planes for $`i=85^{}`$ (fifth column).

Figure 2 – Comparison between profiles derived from model E2 for $`i=50^{}`$ and $`60^{}`$ (thick solid and thick dashed line, respectively) and observations of the Civ 1549Å and Siiv 1403Å transitions in the spectrum of the brightest nova-like variable, IX Vel (Hartley et al. ). The synthesized lines show the absorption component without the contribution from the scattered emission (see the main text).


format: text

text:
# Interfacial mixing in heteroepitaxial growth

## I Introduction

Heterolayers, where e.g. ferromagnets are in contact with antiferromagnets, semiconductors or superconductors, give rise to new ordering and transport phenomena, which depend crucially on the interfacial structure. Examples for such ordering phenomena are the exchange bias , or the cryptomagnetism . Electronic transport through a ferromagnetic-nonmagnetic-ferromagnetic sandwich (“spin valve” geometry) gives rise to the giant magneto-resistance . Another example is the recently predicted possibility to enhance or reduce a Josephson current magnetically by replacing the tunnel barrier of a Josephson junction by a ferromagnetic-insulating-ferromagnetic sandwich . These phenomena belong to the growing field of spintronics , where the spin degree of freedom is used for electronic signal processing. Interfacial mixing affects all of them . For example, it leads to spin scattering disturbing the spin dependent transport properties.

Therefore it is important to be able to control the various physical processes that may spoil well defined interfaces. Some of them proceed after growth such as bulk interdiffusion or chemical interface reactions like silicide formation. However, there are also important processes taking place exclusively at the surface: For example the substrate may partially behave like a surfactant, when one grows a different material on it. It is this latter mechanism which we investigate in this paper. The questions we want to answer concern the asymptotic concentration profile, the width of the interdiffusion zone and possible correlations among the substrate impurities within the growing layer.

Specifically we consider growing some material B on a substrate A. Obviously the interfacial mixing requires that some substrate (A-) atoms get replaced by B-atoms and “float up” on the surface until they get incorporated into the growing film. There is ample experimental evidence that such a behaviour occurs in very different systems like Cr on Fe , AlAs on GaAs , Nb on Fe  or Au on Fe . This process depends on several important parameters including the lattice mismatch and the interaction between the different atoms including magnetic contributions. In particular the explanation of any ordering of A\- and B-atoms close to the interface would require a detailed investigation of these interactions .

The situation becomes considerably simpler, however, if one is interested in the physical properties further away from the interface. Then the concentration of A-atoms may be regarded as sufficiently low that their interaction as well as A-B-ordering become unimportant. The focus on this region justifies our simplified model, in which the interdiffusion zone depends only on the deposition rate $`F`$ of B-atoms, the diffusion constants $`D_\mathrm{A}`$ and $`D_\mathrm{B}`$ of the adatoms of type Aor B, respectively, and the rate $`E`$ for the exchange of B-adatoms with A-atoms. The limit $`E/D_\mathrm{B}\mathrm{}`$, where a B-adatom exchanges with the first A-atom it encounters, would be realized, if the B-adatoms diffuse by an exchange mechanism , while the A-adatoms diffuse by hopping. In the present paper we simplify the model even further by assuming that both kinds of atoms diffuse equally fast on the surface, $`D_\mathrm{A}=D_\mathrm{B}=D`$, with a diffusion constant independent of the surface composition.

Apart from the exchange there is a second crucial ingredient in the model: The A-atoms behave only partially as a surfactant in the sense that they can be overgrown by island edges. By contrast a perfect surfactant atom should “float up” also in front of an advancing island edge.

Naively one would expect an exponential decay of the density profile of A-atoms far from the interface. It is the main result of this investigation that this is not always the case: The incorporation of A-atoms is much slower, giving rise to a power law decay of the concentration profile in the limit of perfect layer-by-layer growth, $`D/F\mathrm{}`$. In this case the width of the interdiffusion zone diverges, provided there are no finite size effects. By contrast, we shall show that for finite $`D/F`$ the width of the interdiffusion zone is no longer infinite, but a power law of $`D/F`$.

This paper is organized as follows. In the next section we are going to define a simple solid-on-solid (SOS) model for epitaxial growth of a B-layer on an A-substrate, which allows for the irreversible exchange of B-atoms with A’s at the surface. In this model the A-atoms on the surface turn out to cluster in a time-periodic self-organized way, which is explained in Sec.III. The next three sections, Sec.IVVI, are devoted to the limit $`D/F\mathrm{}`$. First, in Sec.IV, we present a simple mean field argument leading to the prediction, that the concentration of A’s decays algebraically in the B-layer. For a finite system size this power-law is cut off leading to a finite width $`H`$ of the interdiffusion zone, which is discussed in Sec.V, where also a scaling ansatz for the surface concentration $`c_A`$ is proposed. This scaling ansatz is confirmed by simulation results for one- and two-dimensional surfaces in Sec.VI. The remaining sections deal with interdiffusion for finite $`D/F`$. Sec.VII contains simulation results and scaling arguments, and in Sec.VIII a rate equation theory is developed for the interdiffusion problem. In the Appendix we describe a very efficient implementation of the simulation model for one-dimensional surfaces in the limit $`D/F\mathrm{}`$.

## II The model

In order to model heteroepitaxial growth of a B-layer on an A-substrate including the interfacial mixing, we introduce a simple solid-on-solid (SOS) model, where the lattice mismatch and most interactions between the atoms are neglected. Moreover we describe the exchange mechanism of B-adatoms with A-atoms in the surface simply by a phenomenological constant exchange rate $`E`$, although it depends in reality on the system parameters as well as the local environment.

The model is defined on a simple cubic lattice by the following kinetic rules (cf. Fig. 1):
(1) Starting from an initially flat substrate consisting of A-atoms, B-atoms are deposited at randomly selected sites on the surface with deposition rate $`F`$.
(2) As long as they do not have a lateral neighbor, the B-atoms diffuse on the surface with diffusion constant $`D`$.
(3) When such a B-atoms happens to sit on top of a A-atom, it can exchange vertically with rate $`E`$ or continue to diffuse with rate $`D`$.
(4) After an exchange, the B-atom stays irreversibly bound, whereas the A-atom diffuses on the surface with diffusion constant $`D`$.
(5) There is no back exchange, when a A-atom sits on top of a B-atom.
(6) When two adatoms, regardless of their type, meet, they form a stable non-moving nucleation center of an island.
(7) When an adatom, regardless of its type, reaches a site adjacent to an island, it is irreversibly bound, increasing the size of the island.
(8) Both types of particles can diffuse down across terrace edges without being hindered by an Ehrlich-Schwoebel barrier.
(9) There are no overhangs, i.e., we assume SOS growth.
(10) The exchange of a B-atom with an A-atom underneath is forbidden, if the B-atom is already part of an island, i.e., if it has a nearest neighbor at the same height. Hence A-atoms can be overgrown by island edges.

Measuring time in units of monolayers (ML) and lengths in units of the lattice constant $`a`$, this model is controlled by the two dimensionless parameters $`D/Fa^4`$ and $`r_Ea^2E/D`$. (In the following we set $`a=1`$.)

This model will be investigated for one- and two-dimensional ($`d=1`$ and $`d=2`$) surfaces in the following. Note that it reduces to the usual model for MBE growth (for a recent review see e.g. ), if one does not distinguish the two particle types. Therefore, the $`D/F`$-dependence of quantities like island density, adatom density, surface width, etc. are the same as usual. For general values of $`D/F`$ and $`E/D`$ we use kinetic Monte Carlo simulations  in order to investigate the model. However, for one-dimensional surfaces in the limit $`D/F\mathrm{}`$ we implemented a much more efficient algorithm, which is described in the Appendix.

## III Correlations of the A-atoms

One of the most intriguing qualitative properties of this model is the time-periodic self organization of A-clusters on the growing surface with a period of one monolayer. Fig.2 shows that the A-atoms (black) are first clustered around the nucleation sites of a new layer, but migrate towards the holes remaining in that layer, when the islands coalesce. Thus the characteristic distance between these clusters agrees with the typical distance between the nucleation sites, the diffusion length

$$\mathrm{}_\mathrm{D}(D/F)^\gamma ,$$
(1)
as long as layer-by-layer growth persits. This can be verified by examining the lateral correlations of $`A`$-atoms on the surface after deposition of $`t`$ monolayers:

$$g(\stackrel{}{r},t)=\frac{1}{L^d}\underset{x=1}{\overset{L^d}{}}\rho _A(\stackrel{}{x},t)\rho _A(\stackrel{}{x}+\stackrel{}{r},t)c_A(t)^2,$$
(2)
where $`\rho _A(\stackrel{}{x})`$ denotes, whether there is an $`A`$-atom at the surface at site $`\stackrel{}{x}`$ \[$`\rho _A(\stackrel{}{x})`$ = 1\] or not \[$`\rho _A(\stackrel{}{x})`$ = 0\], and $`c_A(t)`$ denotes the surface density of $`A`$-atoms. $`d`$ is the dimension of the surface. In our simulations, $`d`$ was 1 or 2. A data collapse of these correlation functions for different values of $`D/F`$ is obtained, if the space coordinates are rescaled by $`\mathrm{}_\mathrm{D}`$ (see Fig. 3 for $`d=1`$), which shows that this is the characteristic distance between the A-clusters.

The mechanism of the periodic self organization can most clearly be seen in the limit $`r_\mathrm{E}=E/D\mathrm{}`$, where every $`B`$-adatom exchanges with the first exchange partner $`A`$ it encounters. When layer $`t`$ is completed and layer $`t+1`$ begins to grow, the first $`B`$-atoms deposited are likely to exchange with $`A`$-atoms from layer $`t`$. Hence the nuclei of islands in the new layer $`t+1`$ will consist predominantly of $`A`$-atoms. As growth proceeds, the lower terrace (layer $`t`$) gets depleted from exchange partners $`A`$, either because they are exchanged with freshly deposited $`B`$-adatoms or get overgrown by the islands. Then the core of the islands with a high concentration of $`A`$-atoms gets surrounded by mainly $`B`$-atoms (Fig. 2, left). However, as the island size increases, it becomes more and more likely that $`B`$-atoms are deposited on top of the islands, i.e. in layer $`t+2`$. These $`B`$-adatoms find many exchange partners at the core of the islands, which then become “washed out” and start decorating the island edges, because there are no Ehrlich-Schwoebel barriers in our model. Note that this edge decoration with A-atoms happens without lateral exchange of $`B`$\- with $`A`$-atoms at the edges of the islands, in contrast to the situation studied in . As a result, the interior of the islands gets cleared from $`A`$-atoms, which are collected in the holes of layer $`t+1`$ which get filled last (Fig. 2, right). Then the process starts again: Layer $`t+2`$ nucleates predominantly with $`A`$-atoms which were exchanged from layer $`t+1`$.

For finite $`r_\mathrm{E}`$ the mechanism is similar. However, for vicinal surfaces growing in step flow mode, the correlations among the $`A`$-atoms are different. Here, the terraces get cleared of $`A`$-atoms, which attach to the step edges. This decoration of advancing edges leads to a correlation pattern $`g(\stackrel{}{r},t)`$ with a spatial periodicity identical to the width of the terraces.

## IV The scale free limit

The limit of perfect layer-by-layer growth, $`D/F\mathrm{}`$, is particularly instructive. In this case, there is only one island on the surface, nucleating at a random site. Afterwards at most one adatom can be found on the surface at a time. If we assume maximal exchange in addition, i.e. an exchange rate $`E`$ much larger than both $`D`$ and $`F`$, $`r_E=E/D\mathrm{}`$, the model becomes parameter-free. A-atoms can be buried even in this limit, when they are overgrown by an island edge. However, the last A-atom will never get buried in this case. The nucleus of a new layer will always contain at least one of the remaining surface atoms of type Ain the limit we focus on, because the first B-atom deposited after completion of a layer exchanges with an A-atom before the next B-atom gets deposited: $`c_A(t\mathrm{})=L^d`$. As this case is scale free, we expect that the concentration of A-atoms falls off like a power law into the growing B-film.

In order to get a first idea about the distribution of A-atoms in the growing B-film, it suffices to consider the concentration $`c_A(n)`$ of $`A`$-atoms at the surface after the deposition of $`n`$ monolayers. During the growth of the next layer, a certain fraction $`(1q)[0,1)`$ of these $`A`$-atoms is transported to the next layer via vertical exchange:

$$c_A(n+1)=(1q)c_A(n).$$
(3)
If $`q`$ was constant, this would imply an exponential decay $`c_A(n)(1q)^n`$. In the present model, however, the probability $`q`$ that an $`A`$-atom gets overgrown decreases with decreasing $`c_A(n)`$, resulting in a decay which is slower than exponential.

Qualitatively this can be understood in the following way: In the limit $`D/F\mathrm{}`$ there is only one island on the surface. Moreover, any $`A`$-atom on the lower terrace gets transferred to the new layer as soon as it is reached by an adatom of type $`B`$, because $`r_\mathrm{E}\mathrm{}`$ is assumed, as well. Only $`A`$-atoms sufficiently close to the island have a chance to be overgrown by the island edge before being visited by a $`B`$-atom.

The island edge advances a characteristic distance $`\mathrm{}_{\mathrm{cover}}`$, while the lower terrace gets depleted from $`A`$-atoms by exchange with adatoms of type $`B`$. For a one-dimensional surface of length $`L`$ it is clear that $`\mathrm{}_{\mathrm{cover}}`$ is proportional to the number of $`A`$-atoms at the surface, $`c_A(n)L`$. Hence, the number of $`A`$-atoms with a chance to be overgrown is of the order of

$$\mathrm{}_{\mathrm{cover}}c_A(n)c_A(n)^2L.$$
(4)
This must be compared with $`qc_A(n)L`$, which shows that

$$qc_A(n).$$
(5)
Inserting Eq.(5) into Eq.(3) leads to the difference equation

$$c_A(n+1)c_A(n)c_A(n)^2$$
(6)
implying the asymptotic power law

$$c_A(n)1/n$$
(7)
for the concentration of $`A`$-atoms at the surface.

The concentration profile of $`A`$-atoms inside the grown film is given by

$$\rho (n)=c_A(n)c_A(n+1)1/n^2.$$
(8)
It is remarkable that for the mechanism discussed in this paper the width of the interdiffusion zone diverges logarithmically with the thickness $`T`$ of the film:

$$\underset{n=1}{\overset{T1}{}}\rho (n)n+c_A(T)T\mathrm{ln}T.$$
(9)
Nevertheless, the interface can be localized precisely, because $`B`$-atoms do not occur below layer $`n=0`$ due to the absence of bulk diffusion in this model. Below we show that these power laws are confirmed by simulations.

The argument leading to Eq.(5) ignores that the $`A`$-atoms at the surface are clustered, as shown in Fig.2, and was made plausible only for a one-dimensional surface. However, it can be refined such that it takes these spacial correlations into account and applies also for two-dimensional surfaces. For $`D/F\mathrm{}`$, we can imagine that there is only one cluster of size $`c_A(n)L^d`$ on the surface, when the new layer nucleates. The important point is that the nucleation happens anywhere on the surface with equal probability $`1/L^d`$ in this case. However, only if the nucleation site is within an area of about the size $`2^dc_A(n)L^d`$ centered at the middle of the cluster, there is a chance that some $`A`$-atoms get overgrown. In other words, only a fraction of nucleation sites $`c_A(n)`$ leads to overgrowth. The average number of $`A`$-atoms overgrown in such a case is proportional to the cluster size. Hence on average a fraction $`q`$ of $`A`$-atoms is overgrown which is proportional to the fraction of nucleation sites leading to overgrowth, i.e. this refined argument gives $`qc_A(n)`$ as in Eq.(5).

## V The width of the interdiffusion zone for $`D/F\mathrm{}`$

In the previous section we predicted that the width of the interdiffusion zone diverges in the scale free limit, where $`D/F\mathrm{}`$ and $`r_\mathrm{E}\mathrm{}`$. In this section we predict, that for finite $`r_\mathrm{E}`$ the power law Eq.(7) is only valid, if the system is infinitely large. For finite system size the power law is exponentially cut off at a characteristic width of the interdiffusion zone.

After the completion of several monolayers on a substrate of linear size $`L`$ the number of substrate atoms at the surface is $`c_AL^d`$. The $`A`$-atoms are concentrated in a cluster, which we assume to be compact, hence of diameter $`c_A^{1/d}L`$. This assumption is justified even for $`d=2`$, where the islands initially are fractal, because the $`A`$-cluster occupies the sites which were filled last in the uppermost monolayer. These sites do not form a fractal.

Now we imagine the surface to be coarse grained on the scale of the cluster diameter so that exactly one cell contains the $`A`$-cluster. The typical residence time of an adatom in such a cell is

$$\mathrm{\Delta }t=\frac{(c_A^{1/d}L)^2}{D}.$$
(10)
A $`B`$-adatom which enters the cell containing the $`A`$-cluster will almost certainly be replaced by an exchange partner $`A`$ within the residence time, if

$$E\mathrm{\Delta }t=r_E(c_A^{1/d}L)^21$$
(11)
(exchange dominates). For $`E\mathrm{\Delta }t1`$ the adatom changes from type $`B`$ into type $`A`$ only with probability $`E\mathrm{\Delta }t`$ (overgrowth dominates).

It is plausible to assume that the power law belongs to the exchange dominated slow decay of $`c_A`$ while the exponential cut off indicates the much faster decay when overgrowth dominates. Thus the width $`H`$ of the interdiffusion zone should be reached, when $`c_A`$ becomes so small that exchange is no longer guaranteed, i.e. when $`E\mathrm{\Delta }t`$ drops below 1. Inserting $`c_A1/t=1/H`$ into Eq.(11) one obtains

$$H(\sqrt{r_E}L)^d.$$
(12)
However, as $`c_A`$ cannot become smaller than $`L^d`$, this estimate is only valid for $`r_E<1`$, while $`HL^d`$ for $`r_E>1`$. This is our prediction for the width of the interdiffusion zone in the limit $`D/F\mathrm{}`$. Note that in this limit the cutoff of the power law is a finite size effect: For $`L\mathrm{}`$ the power law extends to infinity.

Based on the results of this and the previous paragraph we can conjecture the following scaling form for the surface concentration of A’s :

$$c_A(t,L;r_E)c_A(t\mathrm{})=\frac{1}{H}f\left(\frac{t}{H}\right)$$
(13)
where according to Eq.(7)

$$f(\tau )1/\tau \mathrm{for}\tau 1.$$
(14)
## VI Numerical Results for $`D/F\mathrm{}`$

In order to check the predictions of Secs. IV and V, we simulated the model described in Sec.II for $`D/F\mathrm{}`$ and varied the values of $`r_\mathrm{E}`$ and system size $`L`$ for one- and two-dimensional surfaces. For the case $`d=1`$ we used the algorithm described in the appendix, while kinetic Monte-Carlo simulations were done for $`d=2`$.

Fig. 4 (for $`d=1`$) and Figures 5, 6 (for $`d=2`$) show the concentration $`c_A`$ of A-atoms at the surface as a funcion of deposition time $`t`$ (in monolayers of B-atoms). All curves are averages over 200 - 400 independent runs. Both for $`d=1`$ and $`d=2`$ the exponent of the power law decay was found to be consistent with the value -1 derived in section IV.

As shown in Fig. 4, the predicted relations Eq.(12) and Eq.(13) lead to the expected data collapse for the one dimensional surface. The results in two dimensions are not as clear. In this case, we obtain the best data collapse with

$$HL^{1.93}r_\mathrm{E}^{1.2},$$
(15)
as shown in the Figures 5 and 6, whereas our predicted exponents (2 and 1, respectively, see Eq.(12)) were about 4% and 20% different. In fact the $`A`$-clusters are not as compact as assumed in the simple argument of Sec.V (see Fig.2).

## VII Scaling for finite $`D/F`$

For finite $`D/F`$ there are many $`A`$-clusters on the surface. Their typical distance is given by the diffusion length $`\mathrm{}_\mathrm{D}`$ as shown in Sec.III. The average size of the $`A`$-clusters is $`c_A\mathrm{}_D^d`$, provided this is much larger than 1. If the concentration $`c_A`$ becomes too small, less and less $`A`$-clusters will be found on the surface, and their typical distance will grow. Finally all $`A`$-atoms will be overgrown, in contrast to the situation of perfect layer-by-layer growth, where the last $`A`$-atom could never be overgrown. Apart from this, one might expect that the results of the previous three sections would essentially remain true, if one replaces $`L`$ by $`\mathrm{}_\mathrm{D}`$. Qualitatively, the surface concentration $`c_A`$ indeed decays first approximately as a power law of the deposition time, which is cut off at a characteristic width $`H`$ of the interdiffusion zone.

Quantitatively, however, the situation turns out to be more complex than this: All exponents are different, as the simulation results (Figs. 710) show. The power law decay of $`c_At^\beta `$ extends over at most two decades for the largest values of $`D/F`$ we simulated, so that the determination of the exponent $`\beta `$ from the slopes in the log-log-plots Fig.7 and 8 is not very accurate. We estimate

$$\beta =\{\begin{array}{ccc}0.78\pm 0.08\hfill & \mathrm{for}\hfill & d=1,\hfill \\ 0.53\pm 0.05\hfill & \mathrm{for}\hfill & d=2,\hfill \end{array}$$
(16)
which are indicated by the dashed lines in the two Figures. Both exponents are significantly smaller than $`\beta =1`$ obtained for infinite $`D/F`$, i.e. $`c_A`$ decays more slowly for finite than for infinite $`D/F`$.

This result is surprising on first sight, because there are more island edges on the surface for finite $`D/F`$, hence more possible places, where $`A`$-atoms may be overgrown. That $`c_A`$ decays more slowly nevertheless, may be explained by the fact, that the nucleation of islands does not happen anywhere with equal probability as for infinite $`D/F`$ , but preferentially far away from the holes in the previous layer, where the A-atoms are concentrated. Therefore overgrowth is less likely, and $`c_A`$ decays more slowly than for infinite $`D/F`$.

This raises the question, how big the parameter $`D/F`$ must be in order to see the exponent $`\beta =1`$ instead of the smaller one. The answer is, that the system size $`L`$ must be small compared to $`\mathrm{}_\mathrm{D}`$ in order to obtain the crossover to the faster decay of $`c_A`$. This was confirmed by simulation results in . In the four Figures belonging to this Section we carefully checked that the system sizes were big enough to exclude finite size effects.

In analogy to Eq.(11) we expect that the power law decay of $`c_At^\beta `$ stops, when

$$r_\mathrm{E}(c_A\mathrm{}_\mathrm{D}^d)^{2/d}1,$$
(17)
or

$$c_A\mathrm{}_\mathrm{D}^d1,$$
(18)
whichever happens first. Replacing $`c_A`$ by $`H^\beta `$, this implies that the width $`H`$ of the interdiffusion zone should be given by

$$H\mathrm{}_\mathrm{D}^{d/\beta }\mathrm{for}r_\mathrm{E}1$$
(19)
and

$$H(\sqrt{r_\mathrm{E}}\mathrm{}_\mathrm{D})^{d/\beta }\mathrm{for}r_\mathrm{E}1.$$
(20)
In analogy to Eq.(13) we postulate then that

$$c_A=\frac{1}{H^\beta }g\left(\frac{t}{H}\right)$$
(21)
with

$$g(\tau )1/\tau ^\beta \mathrm{for}\tau 1,$$
(22)
because $`c_A`$ is independent of $`H`$ for small $`t`$.

We first checked these conjectures for $`r_\mathrm{E}>1`$, where the surface concentration of A’s becomes independent of $`r_\mathrm{E}`$, as expected. If we insert the $`D/F`$-dependence Eq.(1) of the diffusion length in Eq.(19), Eq.(21) can be written in the form

$$c_A\left(\frac{D}{F}\right)^{\gamma d}=g_1\left(t\left(\frac{D}{F}\right)^{\gamma d/\beta }\right).$$
(23)
With $`\gamma =1/4`$  for $`d=1`$ and $`\gamma =1/(4+d_\mathrm{f})0.17`$ for $`d=2`$ ($`d_\mathrm{f}`$ is the fractal dimension of the islands) , and with the $`\beta `$-values determined above, the theory predicts

$$\gamma d=\{\begin{array}{ccc}0.25\hfill & \mathrm{for}\hfill & d=1\hfill \\ 0.35\pm 0.01\hfill & \mathrm{for}\hfill & d=2\hfill \end{array}$$
(24)
$$\gamma d/\beta =\{\begin{array}{ccc}0.32\pm 0.03\hfill & \mathrm{for}\hfill & d=1\hfill \\ 0.66\pm 0.07\hfill & \mathrm{for}\hfill & d=2\hfill \end{array}$$
(25)
The data collapses in Figs.7, 8 are in reasonable agreement with this prediction.

However, the $`r_\mathrm{E}`$-dependence Eq.(20) for $`r_\mathrm{E}<1`$ is not in agreement with the simulation results. For fixed $`D/F`$ the theory Eq.(21) predicts

$$c_Ar_\mathrm{E}^{d/2}=g_2(tr_\mathrm{E}^{d/2\beta }).$$
(26)
Inserting the value of $`\beta `$ determined above, the scaling exponent should be

$$d/2\beta =\{\begin{array}{ccc}0.64\pm 0.06\hfill & \mathrm{for}\hfill & d=1\hfill \\ 1.9\pm 0.2\hfill & \mathrm{for}\hfill & d=2\hfill \end{array}$$
(27)
For $`D/F=10^7`$ we could only check this for about one decade of $`r_\mathrm{E}`$-values: For $`d=1`$ we found that already for $`r_\mathrm{E}=0.3`$ the crossover into the regime, where $`H`$ becomes independent of $`r_\mathrm{E}`$, affects the data. For $`r_\mathrm{E}<10^3`$ the exchange was so weak that the surface concentration of A’s decayed very fast from the beginning, so that a convincing data collapse was not possible. Similar problems occurred for d=2. The best result of our attempts to get a data collapse in the available $`r_\mathrm{E}`$-interval are shown in Fig.9 for $`d=1`$ and Fig. 10 for $`d=2`$. The effective exponents turn out to be very different from the ones predicted in Eq.(26).

## VIII Rate Equation Approach

In this Section we extend the established rate equation approach for submonolayer homoepitaxial growth as described in , in order to apply it to our model for surface interdiffusion. Our approach describes the time evolution of four submonolayer quantities: The density of mobile adatoms $`\rho `$, the total island density $`I`$, the density of mobile adatoms of type $`B`$, $`\rho _B`$, and the density of potential exchange partners of type $`A`$ in the lower layer, $`\rho _A`$. With these quantities, the rate equations are as follows:

$`\dot{\rho }`$ $`=`$ $`FD\rho (I+2\rho ),`$ (28)
$`\dot{I}`$ $`=`$ $`D\rho ^2,`$ (29)
$`\dot{\rho }_B`$ $`=`$ $`FD\rho _B(I+\rho +\rho _B)Ea^2\rho _A\rho _B,`$ (30)
$`\dot{\rho }_A`$ $`=`$ $`Ea^2\rho _A\rho _B\rho _Aa^2D\rho (I+2\rho ).`$ (31)
While the first two equations are identical to the well known point-island model rate equations for homoepitaxial growth, the last two are specific to our heteroepitaxial model. The third one expresses the change in density of the $`B`$-type adatoms. Its positive contribution describes the depostition of new adatoms. The first negative term represents the loss of B-adatoms, when they get incorporated into islands or bind to another adatom to nucleate a new island. The extra term $`D\rho _B^2`$ accounts for the fact that nucleation events involving two B-adatoms count twice as much as those between a B\- and an A-adatom, because they remove two B-adatoms simultaneously. The last term of Eq.(30) describes the exchange of mobile $`B`$’s with $`A`$’s. Equation (31) expresses the annihilation of possible $`A`$-type exchange partners in the lower layer. Since this value is monotoneously decreasing, there is no positive contribution. The negative terms describe the exchange of $`A`$’s with mobile $`B`$’s, and the overgrowth of $`A`$’s due to propagating island edges and nucleation events.

Rescaling the variables (see ) according to $`\widehat{t}=tF\mathrm{}_0^2`$, $`\widehat{\rho }=\rho \mathrm{}_0^2`$, $`\widehat{I}=I\mathrm{}_0^2`$, $`\widehat{\rho }_A=\rho _A\mathrm{}_0^2`$, $`\widehat{\rho }_B=\rho _B\mathrm{}_0^2`$, where $`\mathrm{}_0=(D/F)^{1/4}`$, leads to the dimensionless equations

$`\dot{\widehat{\rho }}`$ $`=`$ $`1\widehat{\rho }(\widehat{I}+2\widehat{\rho }),`$ (32)
$`\dot{\widehat{I}}`$ $`=`$ $`\widehat{\rho }^2,`$ (33)
$`\dot{\widehat{\rho }}_B`$ $`=`$ $`1\widehat{\rho }_B(\widehat{I}+\widehat{\rho }+\widehat{\rho }_B)r_E\widehat{\rho }_A\widehat{\rho }_B,`$ (34)
$`\dot{\widehat{\rho }}_A`$ $`=`$ $`r_E\widehat{\rho }_A\widehat{\rho }_B\left(a/l_0\right)^2\widehat{\rho }_A\widehat{\rho }(\widehat{I}+2\widehat{\rho }).`$ (35)
If we consider systems in perfect layer-by-layer growth mode, this approach not only holds for the submonolayer regime starting from the substrate, but also starting after integer numbers of deposited monolayers. The initial conditions of these equations for a flat surface after $`n`$ deposited monolayers are

$$\widehat{\rho }(0)=\widehat{I}(0)=\widehat{\rho }_B(0)=0,\widehat{\rho }_A(0)=\widehat{c}_A(n).$$
(36)
Disregarding the point island model nature of this approach, which only holds for early stages of the submonolayer regime, we can establish the surface concentration of A’s after the deposition of one additional monolayer, $`c_A(n+1)`$, as the integral over the density of all exchanged atoms:

$$\widehat{c}_A(n+1)=_0^{(\mathrm{}_0/a)^2}𝑑\widehat{t}r_E\widehat{\rho }_A\widehat{\rho }_B$$
(37)
The upper integration boundary is the dimensionless time for depositing one monolayer. This approximation can be justified by taking into account that the transport of A-atoms from the $`n`$th to the $`(n+1)`$th layer mainly takes place at early times, that is, the nucleation regime and early stages of the intermediate coverage regime, as explained in section III. The chosen approach describes these regimes with sufficient accuracy.

The solution of the first two equations can be taken directly from the literature : For early times, $`\widehat{t}1`$, $`\widehat{\rho }`$ is linear in $`\widehat{t}`$, and $`\widehat{I}`$ increases with $`\widehat{t}^3`$. At late times, $`\widehat{t}1`$, one gets $`\widehat{\rho }\widehat{t}^{1/3}`$ and $`\widehat{I}\widehat{t}^{1/3}`$.

With these results, the last two equations can be solved analytically in a similar way for the early-time regime, $`\widehat{t}1`$. For equation (34), the second term on the right hand side can be neglected in this limit. If we also neglect the time dependence of $`\widehat{\rho }_A(t)\widehat{\rho }_A(0)=\widehat{c}_A(n)`$, we get

$$\dot{\widehat{\rho }}_B1r_E\widehat{c}_A\widehat{\rho }_B$$
(38)
This equation relaxes into a steady state with $`\widehat{\rho }_{B,\mathrm{}}1/(r_E\widehat{c}_A)`$ after a characteristic time $`\widehat{t}^{}1/(r_E\widehat{c}_A(n))`$. For even smaller times, $`\widehat{t}\widehat{t}^{}`$, we can also neglect the other right hand side term, and we get $`\widehat{\rho }_B\widehat{t}`$.

Plugging these results into Eq.(35), and realizing that the second term on the right hand side can be neglected compared to the first one, we get

$$\widehat{\rho }_Ae^{\widehat{t}/\widehat{c}_A}(1\frac{\widehat{t}}{\widehat{c}_A})\text{for}\widehat{t}^{}\widehat{t}1$$
(39)
and

$$\widehat{\rho }_Ae^{r_E\widehat{t}^2}(1r_E\widehat{t}^2)\text{for}\widehat{t}\widehat{t}^{}.$$
(40)
To relate these findings to our results in the other sections, we employed an iteration scheme for the rate equation system to obtain the surface concentration of $`A`$’s for every deposited integer monolayer. Starting from the substrate ($`\rho _A(0)=c_A(0)=1`$) we can obtain $`\rho _A(1)=c_A(1)`$ from Eq.(37) by solving Eqs.(32) – (35) numerically. Plugging $`c_A(1)`$ back into our rate equations as the initial surface concentration, that is $`\rho _A(0)=c_A(1)`$, we get $`\rho _A(1)=c_A(2)`$ by using Eq.(37) again. Repeating this iteration scheme leads to $`c_A(t)`$.

The rescaled results of this approach for different $`D/F`$ are shown in Fig. 11. One can clearly observe power law decay for high $`D/F`$ values at intermediate times, and a similar scaling behaviour as obtained from the simulations. The exponent of the power law behaviour is approximately $`1`$, which is identical to the result from the simulations for perfect layer-by-layer growth mode, $`D/F\mathrm{}`$. This fact supports our argumentation concerning the different exponents for $`D/F\mathrm{}`$ and finite $`D/F`$ in section VII: Since the rate equations cannot describe the clustering of the A-atoms, they also do not reflect the preference of the nucleation sites to be far from the A-clusters. The observed scaling exponent of $`l_0^2=\left(D/F\right)^{1/2}`$ doesn’t match any of the exponents resulting from the simulations. This is in general agreement with the analytical calculation of the exponents in section V: There we derived the scaling exponents from characteristic properties of the A-clusters, which are totally neglected in the presented rate equation approach.

## IX Conclusion

In the present work we have investigated heteroepitaxial growth of B-particles on an A-substrate. Introducing an exchange mechanism for B-adatoms, when they encounter an A-atom in the uppermost layer, we observed that in the limit of layer-by-layer growth the top layer concentration of A-atoms decays algebraically. Therefore, the resulting interdiffusion zone has a broad profile with a diverging width. Varying the rate $`E`$ at which A-atoms and B-atoms are exchanged did not change the exponent of this power law. A different situation has been observed, as we varied the diffusion constant $`D`$. For finite values of $`D`$ a crossover from power law to exponential decay has been found. The crossover time $`H`$ is given by $`(D/F)^{0.3}`$ for $`d=1`$ and by $`(D/F)^{0.6}`$ for $`d=2`$.

It would be very interesting to have some experimental results for the concentration of substrate atoms in the growing layer in the case of heteroepitaxy. The results we found suggest that in order to provide a sharp interface between substrate and deposited material one should use small $`D/F`$.

###### Acknowledgements.

We are very grateful to H. Hinrichsen, L. Brendel, V. Uzdin, U. Koehler, C. Wolf and A. Lorke, with whom we had many fruitful discussions. This work was supported by the DFG within Sonderforschungsbereich 491 (Magnetic Heterolayers) and GK 277 (Strukture and Dynamics of Heterogeneous Systems), as well as by the INTAS-project 2001-0386. Work done by SBL was supported in part by Korea Research Foundation Grant (KRF-2001-015-DP0120).

## Appendix

Here we describe, how we implemented the model introduced in Sec. II for one dimensional surfaces in the limit $`D/F\mathrm{}`$. We first describe the idea for the scale free limit, where also $`E/D\mathrm{}`$.

For $`D/F\mathrm{}`$ one has perfect layer-by-layer growth. The nucleation of a new layer happens at an arbitrary position. Afterwards there is at most one adatom on the surface. The idea is to calculate the probabilities exactly, with which the adatom reaches the nearest sinks to its left and to its right. For an A-adatom these are the island edges, while for a B-adatom it might also be an A-atom, with which it could exchange. Let $`d_L`$ ($`d_R`$) denote the distance to the nearest sink to the left (right).

As shown in , the probability $`p_\mathrm{L}`$ to reach the left position prior to the right one with unbiased diffusion is given by

$$p_\mathrm{L}=\frac{d_\mathrm{R}}{d_\mathrm{R}+d_\mathrm{L}}$$
(41)
Correspondingly, $`p_\mathrm{R}=1p_\mathrm{L}`$. Therefore, it is not necessary to simulate the whole random walk of an adatom, but it suffices to select the final position according to (41).

Thus in the scale free limit the model (after nucleation of a new layer) may be simulated as follows:
(1) Deposition of a Bat a randomly chosen site $`i`$.
(2) Determination of the distances $`d_\mathrm{L}`$ and $`d_\mathrm{R}`$ followed by a decision for a side according to the probabilities given in (41).
(3) If the final position of the B-adatom is an A-site, the atoms exchange (as $`E/D\mathrm{}`$). In this case the A-adatom goes to the left or right island edge according to (41). If the final postition of the B-adatom is an island edge, it is bound there irreversibly possibly overgrowing an A-atom. Then one returns to step (1) and deposits the next B-atom.

This algorithm can be generalized for finite $`r_\mathrm{E}=E/D`$: Not always, when a B-adatom encounters an exchange partner A, they exchange immediately. This happens only with probability $`p_\mathrm{E}=E/(E+2D)`$, where the denominator is the sum of the rates for the three possible actions of the adatom – exchange with the A-atom underneath, a hop to the right neighbour and a hop to the left neighbour. With probability $`p_\mathrm{E}`$ the B-adatom is replaced by an A-adatom, which attaches to the island edges to its left with probability Eq.(41), and otherwise to the island edge to its right; $`1p_\mathrm{E}`$ is the probability that the B-adatom continues to diffuse until it encounters the next A-atom or attaches to the island edge.

In order to avoid simulating the random walk explicitely, one has to calculate the probabilities analytically, that the B-adatom exchanges with any particular of the A-atoms or attaches to the island edges. Technically speaking, the B-atom is a random walker on a one-dimensional lattice with fixed partial absorbers (the A-atoms) and two full absorbers (the island edges) (Rosenstock trapping model with partial absorbers). In order to calculate the absorption probabilities at the different absorbers, which depend on the deposition site, we consider an incoming flux (normalized to 1) of independent random walkers at the deposition site $`x_\mathrm{S}`$ (source) and determine the outgoing fluxes at the absorption sites (sinks). The absorption probability is then the steady state fraction of the incoming flux that leaves the system at the respective absorption site.

The density of random walkers at a site $`x`$ evolves according to

$`\dot{\rho }(x,t)`$ $`=`$ $`D[\rho (x1,t)2\rho (x,t)+\rho (x+1,t)]`$ (43)
$`E\rho (x,t)\rho _\mathrm{A}(x)+\delta _{x,x_\mathrm{S}},`$
where the density of partial absorbers, $`\rho _\mathrm{A}(x)`$, is 1 at all the sites $`x_\mathrm{A}`$, where an A-atom sits, and 0 otherwise:

$$\rho _\mathrm{A}(x)=\underset{\nu =1}{\overset{n}{}}\delta _{x,x_{\mathrm{A}\nu }}.$$
(44)
The terms on the right of Eq.(43) which are proportional to $`D`$ are the gain and loss terms due to hopping from a neighbor site to $`x`$, respectively away from $`x`$. The term proportional to the exchange rate $`E`$ describes the loss of walkers at the partial absorption sites. The last term is the gain term due to the normalized influx of walkers at site $`x_\mathrm{S}`$. The perfect sinks corresponding to the island edges are represented by the boundary conditions $`\rho (1)=\rho (L)=0`$, where $`L`$ is the size of the terrace, on which the source is located.

The probability of absorption at site $`x_\mathrm{A}`$ is then obtained from the steady state solution of Eq.(43) by

$$p(x_\mathrm{A})=E\rho (x_\mathrm{A}),$$
(45)
and the ones at the island edges by

$$p(1)=D\rho (2),p(L)=D\rho (L1).$$
(46)
Introducing the diffusion current between $`x`$ and $`x+1`$ (i.e. the current to the right of $`x`$ and to the left of $`x+1`$),

$$j_\mathrm{R}(x)=j_\mathrm{L}(x+1)=D(\rho (x+1)\rho (x)),$$
(47)
equation (43) can be rewritten in the steady state as

$$j_\mathrm{R}(x)j_\mathrm{L}(x)=E\rho (x)\underset{\nu =1}{\overset{n}{}}\delta _{x,x_{\mathrm{A}\nu }}+\delta _{x,x_\mathrm{S}}.$$
(48)
This shows that $`\rho (x)`$ is a piecewise linear function with slope discontinuities at the source and the sinks. Hence Eq.(43) reduces to a set of $`2n+2`$ coupled linear equations for the $`2n+2`$ unknowns $`j_\mathrm{R}(x_{\mathrm{A}\nu }),\rho (x_{\mathrm{A}\nu })`$ and the boundary values $`j_\mathrm{R}(1)`$ and $`j_\mathrm{L}(L)`$.

The solution determines the probabilities Eqs.(45), (46) with which a freshly deposited B-atom is exchanged at the different A-sites or absorbed by the island edges. By choosing a random number we decide which site to pick. If it is an island edge, the B-atom is moved there, and the next B-atom is deposited at a random position. Otherwise we move the B-atom to the chosen site, exchange it with the A-atom there, let another random number determine, whether to attach the A-atom to the left or right island boundary, and deposit the next B-atom at a random position.

The complexity of this algorithm is linear in the number of A-atoms left on the surface, while a brute force simulation of the diffusion would cost much more computing time proportional to $`L^2`$.


text:
# Ferromagnetic transition temperature enhancement in (Ga,Mn)As semiconductor by carbon co-doping

## Abstract

We present a theoretical study of (Ga,Mn)(As,C) diluted magnetic semiconductors with high C acceptor density that combines insights from phenomenological model and microscopic approaches. A tight-binding coherent potential approximation is used to describe the electronic structure in the presence of MnGa and CAs impurities. We find only a small effect of C on the distribution and coherence of electronic states close to the top of the valence band and on the coupling between Mn moments, even at doping levels of several per cent. These results justify applying the model of ferromagnetic Mn-Mn coupling mediated by itinerant holes in the valence band also to C doped samples. The increase of ferromagnetic transition temperature due to the presence of C acceptors is illustrated by calculations that use the $`𝐤𝐩`$ Kohn-Luttinger description of the GaAs valence band and assume systems where Mn local moment and itinerant hole densities can be varied independently.

Prospects for new device functionalities in all-semiconductor spin-electronic structures rely on the realization of a ferromagnetic semiconductor operating at room temperature. An important milestone in this material research was the discovery five years ago ohno\_science98  of ferromagnetism in Mn-doped GaAs firstGaMnAs  with the Curie temperature $`T_c=110`$ K ohno\_science98 . In the Ga1-xMnxAs diluted magnetic semiconductor (DMS) with $`x110\%`$, Mn substituting for Ga provides a local moment $`S=5/2`$ and a delocalized hole bookchapter ; mm03 . Ferromagnetic coupling between Mn moments is mediated by the itinerant holes via a kinetic-exchange interaction dietl . Recent progress in low-temperature MBE growth and post-growth annealing techniques has led to the increase of the transition temperature in (Ga,Mn)As DMS’s by nearly 50 K anneal . This success is attributed to a smaller concentration of carrier and moment compensating defects, especially interstitial Mn, in the optimally annealed samples Mn-inter . The latest observations are consistent with an approximately linear dependence of $`T_c`$ on $`x`$ and hole Fermi wavevector predicted by theory jungwirth\_prb99 ; dietl\_science00 ; jungwirth\_prb02 .

In this paper we address theoretically the possibility of increasing $`T_c`$ in (Ga,Mn)As DMS’s by a non-magnetic acceptor co-doping, namely by introducing substitutional CAs impurities. Our work is partly motivated by a recent experimental observation of a marked enhancement of the Curie temperature in Mn-implanted GaAs:C samples compared to the Mn-implanted undoped GaAs layers park\_tbp . After briefly discussing the different effects of CAs and BeGa Be-Furdyna  acceptors on substitutional Mn incorporation in (Ga,Mn)As DMS’s, we focus on two key issues related to high-density carbon doping: (i) The effect of C impurities on the density of states in the semiconductor valence band is assessed using the tight-binding/coherent-potential-approximation (TB/CPA) description of the disordered semiconductor. (ii) The strength of the hole-mediated ferromagnetic Mn-Mn coupling is compared for systems with and without C co-doping using the TB/CPA results and $`T_c`$ is then estimated from a model that combines the $`𝐤𝐩`$ Kohn-Luttinger description of the GaAs valence bands and a mean-field treatment of the kinetic-exchange coupling between Mn local moments and the band-holes.

An important issue for hole co-doping in (Ga,Mn)As DMS’s is the change of the substitutional MnGa formation energy caused by the presence of additional non-magnetic acceptors. Although the formation energy is, strictly speaking, an equilibrium characteristic we assume that the basic features of its compositional dependence are reflected also in the dynamics of the non-equilibrium growth of the co-doped materials. A systematic ab-inito study of the formation energies is beyond the scope of this short paper and will be discussed elsewhere masek\_tbp . Here we present a qualitative analysis based on the TB/CPA calculations, described in more detail below. Our results suggest that Be acceptors substituting for the same element as Mn, i.e. for Ga, lead to a strong enhancement of the MnGa formation energy. At the same time, the formation energy of interstitial Mn is suppressed and, hence, an increasing fraction of Mn is incorporated in the form of interstitial donors or electrically neutral MnAs or Mn clusters. This scenario, which leads to a decrease rather than an increase of $`T_c`$, has recently been established by extensive experimental studies of Be co-doped (Ga,Mn)As samples Be-Furdyna .

The TB/CPA calculations indicate that C acceptors lead to a much weaker increase of the MnGa formation energy. We surmise that this property stems from the anomalous nature of CAs which acts as an acceptor, yet has a larger Pauling’s electro-negativity than As. The difference between the effects of BeGa and CAs impurities on MnGa incorporation is further enhanced at high doping levels where the competition of Be and Mn for the same lattice site starts to play a role.

The above qualitative analysis of formation energies suggests that CAs co-doping is favorable for achieving high Curie temperatures. Unlike typical acceptors, however, C has a very different atomic size and, as already mentioned, a large electro-negativity compared to the atom it substitutes for. The crucial question, addressed in the following paragraphs, is then to what extent doping by several per cent of CAs changes the semiconductor band structure and whether the model of shallow acceptor carrier-induced ferromagnetism still applies in (Ga,Mn)(As,C).

Our calculations of the electronic structure of GaAs in the presence of MnGa and CAs impurities are done using the tight-binding version of the coherent potential approximation (see e.g. Ref. Masek87 ). The CPA, in contrast to supercell calculations, is well-suited for the mixed crystals with low concentrations of randomly distributed impurities. It provides estimates of a configurationally averaged density of states (DOS) and related quantities. The configurational averaging restores full translational symmetry of the lattice and makes it possible to decompose the DOS into contributions from specific points in the Brillouin zone. The spectral density $`A(k,E)`$ then includes a detailed information about the dispersion of the electronic states in the reciprocal space. In addition, the width of the peaks of $`A(k,E)`$ defines the scattering rate of the band electrons in various parts of the Brillouin zone due to the impurities.

The parameterization of the TB Hamiltonian provides a correct band gap for the pure GaAs crystal Talwar82  and an appropriate exchange splitting of the Mn d–states. Local changes of the crystal potential at both Mn and C impurities, represented by shifted atomic levels, are estimated using Ref. Harrison\_book . Long-range tails of the impurity potentials, which become less important with increasing level of doping, are neglected. (Note, that the Thomas-Fermi screening length is only 3-5 $`\AA `$ for typical carrier densities jungwirth\_apl02 , i.e., comparable to the lattice constant.) Also lattice relaxation effects are neglected within the CPA. This is well justified for MnGa impurities mirbt\_jpcm02 ; masek\_prb03  but becomes a more important issue in the case of CAs impurities. Previous density-functional studies found a -0.6% relative change of the GaAs lattice constant in a 64-atom super-cell with a single CAs impurity latham\_prb01 . Although the lattice relaxation around CAs may change our results quantitatively, we expect that the dominant effects on the band structure and on the spectral broadening arise from the different atomic levels of C compared to As which is readily accounted for in the TB scheme.

In Fig. 1 we plot the spin-polarized DOS in Ga1-xMnxAs1-yCy with $`x=y=4`$% together with local DOS on host (As) and impurity (CAs) atoms. Despite the remarkable difference between the atomic levels of As and C, $`\mathrm{\Delta }\epsilon _p`$ -1 eV, the local DOS on CAs sites (thin full line) near the valence band edge does not differ much from the local DOS on the As sites (dashed line). (A still larger full width of the valence band than $`\mathrm{\Delta }\epsilon _p`$ may partly explain this result.) Similarly, the total DOS (thick full line) is only weakly affected by the presence of C in this spectral region. In particular, the spin-splitting 0.26 eV of the valence band edge, directly related to the kinetic-exchange parameter $`J_{pd}`$, is nearly the same in systems with and without CAs impurities. The shaded region in Fig. 1 shows the DOS of Mn d-states peaked near the energy of -4eV which is consistent with results in C free (Ga,Mn)As.

It is important to point out that the substitution of C also has a very small effect on the line-width of $`A(k,E)`$ for the states close to top of the valence band, i.e., electron scattering on C impurities does not disturb substantially the coherence of the Bloch states. This, together with an unchanged value of $`J_{pd}`$, implies that the additional disorder due to the co-doping with C should not have any marked effect on the carrier mediated coupling between Mn moments. We check this more explicitly using the compatibility of the CPA with the Weiss mean-field theory. The strength of the Mn-Mn coupling is characterized by the energy cost of flipping one MnGa moment, which can be calculated for a given chemical composition masek\_ssc91 . In Fig. 2, we plot this quantity as a function of the hole density, $`p`$, for $`x=y=4`$% (solid line), and $`x=4`$% and $`y=0`$ (dashed line) samples. The hole density is varied in the calculations indepenedently of $`x`$ and $`y`$, i.e., we assume implicitly a compensation whenever $`p`$ is smaller than the total density of MnGa and CAs acceptors. The curves are nearly identical for systems with and without C, as anticipated above. This means that, for a given hole concentration, the exchange coupling is quite insensitive to the large amount of the additional CAs defects. We conclude our TB/CPA considerations by stating that, as in the MnGa acceptor case, the main effect of CAs doping is a downward shift of Fermi energy with respect to the band edge, i.e., an increase of the number of holes in the valence band. In the following paragraphs we discuss prospects for ferromagnetic transition temperature enhancement by adding extra holes into the DMS valence band.

We use a model in which the itinerant holes are described by the GaAs host bands and the coupling to the local moments by a phenomenological constant $`J_{pd}`$ dietl ; ohno\_jmmm99 . This theory has been successful in describing semiquantitatively many non-trivial thermodynamic and transport properties of (Ga,Mn)As DMS’s dietl\_science00 ; jungwirth\_prb02 ; aniso ; domains ; jungwirth\_prl02 ; jungwirth\_0302060 ; sinova\_prb02 . In the simplest, virtual-crystal mean-field version of the model, that assumes ferromagnetic indirect coupling between Mn ions, the energy to flip the MnGa moment is proportional to the effective field $`H_{eff}=J_{pd}s`$. Here $`s`$ is the mean spin-polarization density of the itinerant holes which increases with the hole density. This result is consistent with our TB/CPA calculations, shown in Fig. 2, for hole densities not too much larger than the MnGa density. Recall that the concentration of the substitutional Mn ions is given by $`N_{Mn}=4x/a_{lc}^3`$ where $`a_{lc}`$ is the GaAs lattice constant ($`N_{Mn}=0.88`$ nm-3 for $`x=4`$%, e.g.). For $`p>N_{Mn}`$, RKKY oscillations of the Mn-Mn coupling dietl  start to play a role and the increasing number of antiferromagnetically coupled Mn moments leads to a saturation or even to a suppression of $`T_c`$, as seen in Fig. 2 for $`p>1.3`$ nm-3. The mean-field theory that allows only for collinear ferromagnetic states is therefore likely to break down in this high hole-density region. Note also that for high $`N_{Mn}`$ and low $`p`$, the direct antiferromagnetic Mn-Mn interaction takes over, as suggested by negative Weiss exchange-field values in Fig.2, which sets another limit on the validity of the simple mean-field model. This low hole-density region is, however, not important for the high-$`T_c`$ (Ga,Mn)(As,C) DMS’s we focus on in this paper.

We will use now the mean-field model to estimate the Curie temperature as a function of the MnGa local moment density and of the density of itinerant holes that can be varied independently. We emphasize that hole doping due to C is relevant also for systems with $`p<N_{Mn}`$ since most of the experimental (Ga,Mn)As samples show some level of compensation, usually caused by interstitial Mn defects. In the calculations, the GaAs host band-structure is obtained from the $`𝐤𝐩`$, Kohn-Luttinger model aniso . We neglect the $`1020`$% suppression of $`T_c`$ due to spin-wave fluctuations since the effect is compensated, to a large extent, by hole-hole exchange enhancement of $`T_c`$ jungwirth\_prb02 , also neglected in the present calculations. Details of the model are described elsewhere jungwirth\_prb99 ; bookchapter ; aniso ; jungwirth\_prb02 ; here we recall only the general $`T_c`$ expression we use dietl ; jungwirth\_prb99 ; jungwirth\_prb02 :

$$k_BT_c=\frac{N_{Mn}S(S+1)}{3}\frac{J_{pd}^2\chi }{(g\mu _B)^2},$$
(1)
where $`\chi `$ is the band-hole magnetic susceptibility which is roughly proportional to $`p^{1/3}`$.

Fig. 3 shows constant-$`T_c`$ curves calculated for critical temperatures ranging from 50 K to room-temperature. The dashed line corresponds to $`p=N_{Mn}`$ and also indicates approximately the onset of RKKY oscillation effects on $`T_c`$. Note that the simple mean-field theory used to calculate the constant-$`T_c`$ curves is expected to break down above the dotted-dashed line where the large number of antiferromagnetically RKKY-coupled Mn moments might lead to a decrease rather than increase of $`T_c`$ with increasing hole density.

Fig.3 suggests that a substantial enhancement of the ferromagnetic transition temperature may be expected in hole co-doped samples with high Mn moment concentration. Assuming, e.g., $`x=10`$%, an increase of the hole density from $`p0.1`$ nm-3 to $`p0.7`$ nm-3 leads to an increase of the theoretical $`T_c`$ from 50 K to 300 K. For smaller $`x`$, a larger $`\mathrm{\Delta }p`$ is needed to enhance $`T_c`$ by the same amount. Recent experiments have demonstrated that good quality (Ga,Mn)As DMS’s can be grown with $`x`$ reaching 8% anneal . Based on the curves in Fig. 3 and our TB/CPA results we conclude, that (Ga,Mn)As DMS’s co-doped with several per cent of C should not be overlooked among potential candidates for a room-temperature ferromagnetic semiconductor.

###### Acknowledgements.

We acknowledge fruitful discussions with Cammy Abernathy, Miroslav Cukr, Vít Novák, and Dan Park. This work was supported in part by the Welch Foundation, DOE under grant DE-FG03-02ER45958, the Grant Agency of ASCR under grant A1010214, and the Grant Agency of the Czech Republic under grant 202/02/0912.


text:
# Constraining superheavy dark matter model of UHECR with SUGAR data

## I Introduction

The origin of ultra-high energy cosmic rays (UHECR) is a long-standing puzzle. Observations of cosmic rays with energy above $`10^{19}`$ eV reveal that the GZK cutoff gzk  in the energy spectrum is not observed by AGASA agasa  (though it seems to be consistent with HiRes hires ), and that the arrival directions of cosmic rays are quite isotropic at large scales but cluster at small anglescluster . This is a perplexing situation for seeking UHECR sources. The isotropic distribution of UHECR arrival directions, in the absence of strong magnetic field which can produce it, either Galactic or extraGalactic, implies that the sources are distributed over a few hundred Mpc for which the matter distribution looks homogeneous, while the absence of the GZK cutoff implies that the UHECR sources are within the GZK radius which is about 50 Mpc. The observed correlations of UHECR with BL Lacertae objects bllacs  do not make the situation easier: although perfectly consistent with clustering, BL Lacs as sources of UHECR cannot explain super-GZK events under assumption of standard physics Kachelriess:2003yy .

A possible way to reconcile these contradictory facts is to assume that, apart from UHECR accelerated in BL Lacs, there exists a second, GZK-violating component Berezinsky:2002vt  originating from decays or annihilation of superheavy particles within the GZK radius shdm ; bdk02 . These superheavy particles are produced in the early universe kt99  and constitute a fraction of cold dark matter. The latter is concentrated in galactic halos. Since we are 8.5 kpc away from the Galactic center, the clear signature of superheavy particle decays is the significant excess of the UHECR flux from the direction of the Galactic center dt98 . This kind of anisotropy toward the Galactic center (or toward any other strong source) has not been observed by the Northern hemisphere arrays Hayashida:1999ab ; bsw99 . However, Northern arrays do not see the Galactic center directly, and the small statistics of cosmic rays with energy above the GZK cutoff makes it difficult to draw strong constraints.

In this paper we use the SUGAR data sugar  to examine the viability of the superheavy dark matter (SHDM) scenario of UHECR, for SUGAR is the unique experiment which operated in the Southern hemisphere before Pierre Auger Observatory and covered the Galactic center. We use the Kolmogorov-Smirnov (K-S) test to check the consistency of the SUGAR data with the arrival direction distribution expected from superheavy particle decays in the Galactic halo. Our purpose is to limit a fraction of UHECR which are produced by this mechanism.

We perform the analysis in two steps. In Sects. II and III we analyze the SUGAR data and derive bounds on the relative weight $`\eta `$ of the halo component in the total UHECR flux. These bounds are presented in Figures 3 and 4. They depend exclusively on the halo profile and statistics of SUGAR data. We then move in Sect. IV to the implications of these constraints for the superheavy dark matter models. The difference with the previous case is that SHDM contribution contains both Galactic and extragalactic parts. Their relative weight depends on a particular SHDM model (spectrum and maximum energy) and on absolute dark matter density in the halo of our Galaxy. Because of the uncertainties in these parameters we separate this question from the more robust analysis of sects.II and III. We will see, however, that in SHDM models the extragalactic contribution is generically subdominant in the energy range of interest, and therefore boudns of Secs.II and III apply to SHDM part of cosmic rays as a whole. We present our conclusions in Sect.V.

## II Galactic halo models and SUGAR data

We assume that superheavy particles constitute a part of cold dark matter and their spatial distribution follows the dark matter density. We examine two different Galactic halo models. The first one is the isothermal (ISO) halo model iso , which is characterized by the density profile

$$\rho _{\mathrm{ISO}}(R)\frac{1}{R^2+R_c^2}.$$
(1)
The second one is the Navarro-Frenk-White (NFW) model nfw96 , which is obtained from the $`N`$-body simulation and thus can be more realistic

$$\rho _{\mathrm{NFW}}(R)\frac{1}{R(R+R_s)^2}.$$
(2)
Here $`R`$ is the distance from the Galactic center, $`R_c`$ is the scale of the halo core, and $`R_s`$ is the characteristic scale of the NFW model. For our galaxy, the size of $`R_c`$ is supposed to be a few kpc. The value of $`R_s`$ is expected to be $`10`$$`20`$ kpc on the basis of correlation between the predicted total mass and the value of $`R_s`$ nfw96 . In our analysis, we leave $`R_c`$ and $`R_s`$ as free parameters and use $`R_c=4`$ kpc and $`R_s=15`$ kpc as prototype values.

For cosmic rays coming from the Galactic halo, the expansion of the universe and the evolution of halo can be neglected. We also neglect the effects of the Galactic magnetic field on the propagation of UHECR and will make a comment on it afterward. Since we consider the halo model which is spherically symmetric about the Galactic center and the distribution of sources follows the halo profile, the cosmic ray flux can be calculated as

$$f(\theta )=\frac{1}{4\pi }_0^{r_{\mathrm{max}}(\theta )}L(\sqrt{r^22rR_0\mathrm{cos}\theta +R_0^2})𝑑r,$$
(3)
where $`L(R)`$ is the luminosity of cosmic ray sources, $`\theta `$ is the angle measured from the Galactic center direction, which is related to the Galactic coordinates $`l`$ and $`b`$ by $`\mathrm{cos}\theta =\mathrm{cos}b\mathrm{cos}l`$, $`r_{\mathrm{max}}(\theta )=R_0\mathrm{cos}\theta +\sqrt{R_{\mathrm{max}}^2R_0^2\mathrm{sin}^2\theta }`$, $`R_{\mathrm{max}}`$ is the extension of the Galactic halo, and $`R_08.5`$ kpc is the distance to the Earth from the Galactic center. In the case of decay, the luminosity $`L(R)`$ is proportional to the number density of superheavy particles, that is $`L(R)n_{\mathrm{SH}}(R)\rho _{\mathrm{Halo}}(R)`$. In the case of annihilation, $`L`$ is proportional to $`n_{\mathrm{SH}}^2`$. Since the density profile has a peak at the Galactic center, the annihilation case predicts stronger anisotropy than the decay case and leads to stronger constraints on the halo contribution to UHECRs. In the following, we will consider the decay case only.

To obtain the arrival direction distribution predicted from the Galactic halo model, we use the Monte-Carlo simulation based on the flux (3). In simulating cosmic ray data detected by an array of detectors, we also need to consider the exposure function of the detector array. Here, we use the simple geometrical exposure function which is a function of declination $`\delta `$ only

$$h(\delta )=\frac{1}{\pi }\left[\mathrm{sin}\alpha _m\mathrm{cos}\lambda \mathrm{cos}\delta +\alpha _m\mathrm{sin}\lambda \mathrm{sin}\delta \right]$$
(4)
where $`\lambda `$ is the latitude of the detector array, $`\theta _m`$ is the zenith angle cut,

$$\xi =\frac{\mathrm{cos}\theta _m\mathrm{sin}\lambda \mathrm{sin}\delta }{\mathrm{cos}\lambda \mathrm{cos}\delta },$$
(5)
and $`\alpha _m=0`$ for $`\xi >1`$, $`\alpha _m=\pi `$ for $`\xi <1`$ and $`\alpha _m=\mathrm{cos}^1\xi `$ otherwise.

The SUGAR array, located in the Southern hemisphere, has the latitude $`\lambda =30.53^{}`$. It covered the Galactic center direction best so far, because other detectors are all located in Northern hemisphere and cover the Galactic center marginally at best. Therefore, the SUGAR data is most suited at the moment for the study of cosmic rays coming from the Galactic halo, though in the near future Pierre Auger array will provide much better data. The published SUGAR data contain events with zenith angle up to 72 sugar . We use the data with zenith angle smaller than 60 considering increasing uncertainty at large zenith angles. Figure 1 shows arrival directions of 316 SUGAR events with zenith angles $`z60^{}`$ and energy $`E1.0\times 10^{19}`$ eV (upper panel), and simulated data obtained using the flux (3) for two halo models (lower panels). Note that the NFW model with $`R_s=15`$ kpc and the isothermal model with $`R_c=4`$ kpc predict similar concentration of events around the Galactic center.

## III Kolmogorov-Smirnov test

The most direct way to compare the arrival direction distribution of the real data and that predicted by the superheavy dark matter model is to use the two-dimensional K-S test over the pair $`(\alpha ,\delta )`$ or $`(l,b)`$. Here, we use a different method which is simpler and gives better results. Since the cosmic rays in the halo model are distributed symmetrically around the Galactic center, we analyze the distribution of cosmic rays over the angle $`\theta `$ between their arrival directions and the direction to the Galactic center. We use the one-dimensional K-S test to compare the SUGAR data and model predictions. Since one-dimensional K-S test is reparametrization-invariant, we use $`x=\mathrm{cos}\theta `$ which is more convenient in computation.

The K-S test makes use of the cumulative probability distribution function, which is defined by

$$S_N(x)=\frac{1}{N}\underset{i}{}\theta (x_ix),$$
(6)
where $`N`$ is the number of data points. For the distribution of arrival directions expected in the Galactic halo model, we generate large number $`M`$ of simulated events and obtain corresponding cumulative distribution $`S_M(x)`$. We then calculate the K-S statistic

$$D=\mathrm{max}\left|S_N(x)S_M(x)\right|,$$
(7)
which determines the probability that the two sets of events are drawn from the same distribution nrc . Small probabilities indicate that the two sets significantly differ. Figure 2 shows the cumulative distribution functions used in the K-S test. The difference between the SUGAR data and halo models is clearly seen (it can be expected already from Figure 1).

As mentioned in the Introduction, the SHDM model does not explain clustering and correlations of UHECR with BL Lacs, so another (extragalactic) component is needed which is isotropic at large angles. The SHDM model itself contains the isotropic extragalactic contribution coming from SHDM decays in other galaxies. Thus, let us consider a generic model where UHECR flux consists of the isotropic and Galactic halo components (their relative strength in the SHDM models is discussed in the next section). An important quantity is the fraction $`\eta `$ of cosmic rays coming from the Galactic halo,

$$\eta (E_{\mathrm{min}})=\frac{_{E_{\mathrm{min}}}f_\mathrm{G}(E)𝑑E}{_{E_{\mathrm{min}}}\left(f_\mathrm{G}(E)+f_\mathrm{U}(E)\right)𝑑E},$$
(8)
where $`f_\mathrm{G}(E)`$ and $`f_\mathrm{U}(E)`$ are the energy spectra of cosmic rays from the Galactic halo and from the uniform background, respectively. We consider $`\eta `$ as a function of the lower energy cut $`E_{\mathrm{min}}`$. When two spectra $`f_\mathrm{G}(E)`$ and $`f_\mathrm{U}(E)`$ have quite different shape, $`\eta (E_{\mathrm{min}})`$ depends sensitively on the value of $`E_{\mathrm{min}}`$. Instead of following a specific model to evaluate $`\eta (E_{\mathrm{min}})`$, we use $`\eta `$ as a parameter and calculate the probability to get the SUGAR data distribution. In this way bounds on $`\eta `$ at a given energy cut can be obtained.

Figure 3 shows the dependence of the probability that all cosmic rays come from the Galactic halo ($`\eta =1`$) as a function of the energy cut $`E_{\mathrm{min}}`$ and the halo parameters $`R_c`$ and $`R_s`$. The dependence on $`R_s`$ in the case of NFW halo profile is rather weak. At realistic values of parameters $`R_c4`$ kpc and $`R_s15`$ kpc the two halo models give similar constraints. Dominant fraction of halo component is excluded up to energies as high as $`5\times 10^{19}`$ eV. For higher energies, the constraints become progressively weaker; the reason may be much smaller statistics.

Figure 4 shows the lines of constant probability (corresponding to $`1\sigma `$$`4\sigma `$ levels) in the plane of the model parameters $`R_c,R_s`$ and the fraction of halo component $`\eta `$. Two energy cuts are represented: $`E_{\mathrm{min}}=1\times 10^{19}`$ eV (upper panels) and $`E_{\mathrm{min}}=4\times 10^{19}`$ eV (lower panels). In the case of NFW halo profile, isotropy of SUGAR data at $`E>1\times 10^{19}`$ eV excludes fractions of halo component larger than 15% at $`2\sigma `$ level, while for $`E>4\times 10^{19}`$ eV the halo fraction has to be less than 45%. Constraints are similar for isothermal halo profile eq.(1) at $`R_c4`$ kpc.

Finally, let us comment on how robust are the results presented here. One thing one should check is the influence of the Galactic magnetic field during the propagation. This can be important because significant portion of cosmic rays from the Galactic halo are produced near the Galactic center and experience strong magnetic fields. To see the effect of this spread on our results, we added a Gaussian random deflections with mean deviation up to $`10^{}`$ in the simulated data set. This effect also mimics possible poor angular resolution of SUGAR events. We found that results presented above are stable against this modification. The maximum change occurs at $`E_{\mathrm{min}}1\times 10^{19}`$ eV and $`\eta 1`$. For the moderate energy cut or with the uniform background component added, its effect becomes less than 10% and can be completely neglected.

## IV Application to SHDM model

In order to see what the above constraints imply for a given SHDM model, one has to estimate the strength of the halo (anisotropic) component relative to the extragalactic (isotropic) one. The latter consists of a part related to SHDM decays outside of the Galaxy and a part due to astrophysical sources. Let us estimate the former. The Galactic and extragalactic contributions in the SHDM model can be written as

$$F_G(E)=Cd^3x\frac{n(x)}{(\stackrel{}{R}\stackrel{}{x})^2},$$

$$F_{\mathrm{ext}}(E)=Cd^3x\frac{\overline{n}}{r^2}w(r).$$

Here $`\overline{n}`$ and $`n(x)`$ are the average and halo dark matter densities, $`\stackrel{}{R}`$ is the vector from Galactic center to the Earth, while $`w(r)`$ is the weight function which accounts for the flux attenuation at large distances. The proportionality coefficient $`C`$ which accounts for the efficiency of UHECR production from SHDM is the same in both cases. It therefore cancels in the ratio of the two contributions.

The extragalactic part can be written in the form

$$F_{\mathrm{ext}}=2\pi C\overline{n}R_{\mathrm{eff}},$$

where by definition $`R_{\mathrm{eff}}=𝑑rw(r)`$. This quantity represents the size of the region from which the extragalactic flux is collected (for instance, for very high-energy protons $`R_{\mathrm{eff}}`$ is the GZK distance). It can be estimated as the attenuation length divided by the spectral index. The factor $`2\pi `$ in this equation takes into account the fact that SUGAR sees roughly half of the sky.

To estimate the halo contribution we take the NFW profile at $`R_s=15`$ kpc and calculate the flux coming from half of a sphere overlooking the Galactic center. The result can be written as follows,

$$F_G2\pi C\overline{n}\times 5000\text{ Mpc}.$$

When obtaining this number the ratio of the local to average dark matter density was taken to be $`2\times 10^5`$. This number is uncertain by a factor 2 or so.

The ratio of the two components is

$$\frac{F_{\mathrm{ext}}}{F_G}\frac{R_{\mathrm{eff}}}{5000\text{ Mpc}}.$$

The UHECR flux in the SHDM models is dominated by photons which have attenuation length less than 100 Mpc in the range of energies $`10^{19}10^{20}`$ eV even at very weak assumptions about radio photon background BhattSigl . We therefore arrive at the conclusion that in the SHDM models the halo component dominates already at energies around $`10^{19}`$ eV, i.e., quite below the GZK cutoff. Since the halo component dominates, the bounds on the fraction $`\eta `$ in Figures 3 and 4 essentially apply to the total fraction of UHECR which originate from decays of SHDM particles.

The situation becomes even better at high energies where the extragalactic component dies away and the halo contribution reaches 100%. In this energy range $`60`$ uniformly distributed events would be enough to rule out the SHDM model at $`45\sigma `$ level following the technique developed here.

## V Conclusions

We have used the SUGAR data to check the viability of the hypothesis that a part of the highest-energy cosmic rays originate from decays of relic SHDM particles. We found that the absence of the anisotropy toward the Galactic center in the SUGAR data imposes strong constraints on corresponding fraction of UHECR. At 95% confidence level, the fraction of SHDM-related cosmic rays should be less than 20% and 50% at $`E>1\times 10^{19}`$ eV and $`E>4\times 10^{19}`$ eV, respectively. This implies that the SHDM component of UHECR, if exists, must have harder spectrum than the uniform one and should take over at energies higher than $`4\times 10^{19}`$ eV.

###### Acknowledgements.

The authors are grateful to S. Dubovsky for helpful discussions. This work is supported by the Swiss Science Foundation, grant 20-67958.02.
Note Added: While completing this work, we received a preprint by M. Kachelrieß and D. V. Semikoz ks03  where similar results were obtained. We thank M. Kachelrieß and D. V. Semikoz for sending us the preprint before publication.