text:
# The Cosmic Infrared Background at 1.25 𝜇m and 2.2 𝜇m using DIRBE and 2MASS: a contribution not due to galaxies ?

## 1 Introduction

The cosmic infrared background (CIB) is an important cosmological constraint on the star formation history in the Universe (Dwek et al., 1998; Gispert et al., 2000). The first detection of the CIB was reported by Puget et al. (1996) for sub-millimetric wavelengths using data from the Far Infrared Absolute Spectrometer (FIRAS) on-board the Cosmic Background Explorer (COBE). Shorter wavelengths have been investigated with the help of the COBE Diffuse Infrared Background Experiment (DIRBE) from 240 $`\mu `$m to 1.25 $`\mu `$m. Hauser et al. (1998a) measured the CIB at 240 and 140 $`\mu `$m and set upper limit for other wavelengths. Lagache et al. (2000) proposed an estimation of the CIB at 100 $`\mu `$m.

The CIB is the signal that remains after subtracting from the total celestial brightness the emission from the interstellar medium, stars and interplanetary dust (scattering and emission). However, estimating the contamination due to these three foreground components is challenging. The interstellar medium is the main contaminant for long wavelengths whereas starlight and zodiacal light dominate at wavelengths studied in the present paper ($`J=1.25`$ $`\mu `$m and $`K=2.2`$ $`\mu `$m).

Stellar population models were used to estimate the brightness of starlight in the large beam DIRBE data (Arendt et al., 1998) but this technique led only to upper limits. Detections of the CIB at 3.5 $`\mu `$m ($`L`$) and 2.2 $`\mu `$m (K) were reported by Gorjian et al. (2000) using near-infrared observations of a $`2^{}\times 2^{}`$ dark spot to measure the brightness of stars in a few DIRBE beams. Wright (2001) confirmed the result at 2.2 $`\mu `$m and proposed a weak limit for 1.25 $`\mu `$m ($`J`$) using data from the 2 Micron All-Sky survey (2MASS) for 4 dark spots. Matsumoto et al. (2000) estimated the CIB from 1.4 $`\mu `$m to 4 $`\mu `$m with the Near Infrared Spectrometer (NIRS) on-board the InfraRed Telescope in Space (IRTS). Although NIRS provided a point source catalog limited to stars brighter than 7.5 mag at 2.24 $`\mu `$m, significant stellar contamination remains and the CIB estimate is strongly dependent on the stellar population model used to remove contribution from fainter stars.

This work presents a large scale study of the CIB using 2MASS data for 1400 deg2 of the sky in order to estimate accurately the stellar contribution to the surface brightness observed by DIRBE. Section 2 describes the 2MASS and DIRBE data. The method used to compare these two data sets is presented is Section 3, and the following section presents the determination of the CIB for $`J`$ and $`K`$ with the associated uncertainties. Section 5 is dedicated to the comparison with previous results and with galaxy counts.

## 2 Data

### 2.1 2MASS

The 2MASS second Incremental Data Release from the 2MASS survey covers 48% of the sky and contains $`1.6\times 10^8`$ stars (Cutri et al., 2000). The Point Source Catalog completeness limits are 15.8 and 14.3 mag with a signal-to-noise ratio greater than 10:1 for $`J`$ and $`K_s`$, respectively. No photometry is provided for stars brighter than 4-5 mag (see Sect. 3.2).

Since the catalog of the whole release is $`7`$ GB in binary format, we choose to work on integrated maps in which each pixel corresponds to the integrated flux in a $`5^{}\times 5^{}`$ box. These maps are constructed from the point source catalog clipped to remove sources fainter than the completeness limit. Magnitudes are converted to flux density using the flux for zero magnitude $`F_0(J)=1603`$ Jy defined by Campins et al. (1985) and $`F_0(K_s)=698`$ Jy which is an extrapolation of the $`J`$ value to the $`K_s`$ filter assuming a black body at 9700 K (Vega photometric system). Uncertainties on these photometric zero points are about 2% and do not affect our results. Each $`5^{}\times 5^{}`$ box contains the sum of the point source fluxes in the area.

The resulting image in an Aitoff (equal area) projection consists of $`4321\times 2161`$ pixels. This size is reasonable and the image can be handled globally. Errors resulting from this pixelization are negligible since $`5^{}`$ is small compared to the DIRBE beam of $`1^{}`$ (see below).

We should point out that the 2MASS $`K_s`$ filter (2.17 $`\mu `$m) and the DIRBE $`K`$ filter (2.2 $`\mu `$m) are slightly different. Persson et al. (1998) discuss the $`K`$ to $`K_s`$ transformation and present a list of photometric standards in both colors. Moreover, as Pahre & Mould (1994) point out, the $`K_s`$ bandpass is intermediate between $`K`$ and $`K^{}`$ and rough corrections of $`K_s`$ to $`K`$ can be obtained by averaging between $`K^{}`$ and $`K`$. This transformation is well defined by Wainscoat & Cowie (1992): $`K^{}K=(0.22\pm 0.03)(HK)`$ for $`0<HK<0.4`$ ($`K^{}K0.07`$ otherwise). Hence, for a K5 spectral type star we would obtain $`K_sK0.01`$ for example. This correction is negligible compared to uncertainties described further in the paper and differences between 2MASS and DIRBE filters will be ignored in the following analysis.

### 2.2 DIRBE

Cryogenic DIRBE operation took place from 1989 December to 1990 September. During these 10 months, the sky was observed at the rate of half of the sky per week in 10 bands from 1.25 $`\mu `$m to 240 $`\mu `$m (Hauser et al., 1998a). The DIRBE instrument was designed to make accurate absolute sky brightness measurements, with a straylight rejection of $`<1`$ nW m-2 sr-1 and an absolute brightness calibration uncertainty of 0.05 and 0.03 nW m-2 sr-1, at 1.25 and 2.2 $`\mu `$m, respectively.

The zodiacal light is the first component along the light of sight and must be removed using a dust interplanetary model (e.g. Kelsall et al., 1998). However these models are not unique (see Wright, 1998; Gorjian et al., 2000 for another model). Artifacts still remain. In the following analysis, we use the Zodi-Subtracted Mission Average (ZSMA) maps produced by the DIRBE team (Hauser et al., 1998b). The zodiacal light intensities were subtracted week by week and the residual intensity values were averaged to create the ZSMA maps.

To compare DIRBE and 2MASS, we need a precise knowledge of the DIRBE beam in both $`J`$ and $`K`$ bands. The $`42^{}`$ beam size commonly mentioned in the literature is only valid for daily maps, not the annual maps used here. The beam for the $`J`$ band is presented in the Figure 1. This is an effective beam profile (provided with DIRBE maps), which measures the relative response of DIRBE to a point source, including the effects of sky-scanning and data-sampling rates. The beam for the annual maps corresponds to an average of the effective beam for all orientations (Fig. 1). Similar results are obtained for $`J`$ and $`K`$ and the full width half maximum of both beams is $`1^{}`$ (solid angle $`0.78`$ deg2).

## 3 Method

### 3.1 2MASS integration on DIRBE pixels

First, we exclude from this study regions of the sky for which we know that the CIB cannot be straightforwardly extracted, for example due to high stellar density fields, interstellar cirrus (which is responsible for emission and scattering), and regions of residual zodiacal light structures in the DIRBE maps.

We keep only regions that satisfy the following criteria:

* high galactic latitude: $`|b|>40^{}`$
* low DIRBE 240 $`\mu `$m flux in order to eliminate diffuse emission or scattered light from cirrus clouds: $`I_{240}<3`$ MJy sr-1
* high ecliptic latitude: $`|\beta |>30^{}`$
* exclusion of the Magellanic clouds

The most straightforward way to remove the stellar component (2MASS) from the DIRBE flux is to work directly in the COBE coordinate system (projection on a cube) and to integrate the 2MASS maps on each DIRBE pixel. However, since the 2MASS data still contain coverage holes, we kept only those that are completely covered by 2MASS. The effective useful region covers finally about 1400 deg2 (Fig. 2).

### 3.2 Bright stars

In the remaining 1400 deg2 piece of sky, $`25`$% of pixels are contaminated by bright stars. Stars bright enough to be detected as point sources in DIRBE are saturated in 2MASS. Consequently their photometry cannot be derived from 2MASS to be subtracted from DIRBE.

Moreover, the DIRBE beam has a size of $`1^{}`$, but the sampling of the data is $`22^{}`$. Several pixels are therefore contaminated by a bright star, and this contamination depends on the star position in the beam. To identify pixels affected by these stars we use a median filter which deals with the 2MASS holes. A pixel is considered as contaminated if it is brighter than the threshold, $`𝒯`$, defined by the median plus $`2\sigma `$, $`M+2\sigma `$, of the surrounding pixels in a $`3^{}`$ radius circle of the DIRBE$``$2MASS map:

$$𝒯^i=M^i\left(\underset{i,r=3^{}}{}\mathrm{DIRBE}2\mathrm{M}\mathrm{A}\mathrm{S}\mathrm{S}\right)+2\sigma ^i\left(\underset{i,r=3^{}}{}\mathrm{DIRBE}2\mathrm{M}\mathrm{A}\mathrm{S}\mathrm{S}\right)$$
(1)
Pixels fainter than $`M2\sigma `$ are also removed to avoid bias. In order to validate this operation, we use the 2MASS Bright Infrared Star Compilation (Tam et al. 1999111http://spider.ipac.caltech.edu/staff/raymond/2mass/birsc/birsc.html) which contains photometry for stars as bright as $`4`$ mag in $`JHK_s`$. This compilation consists of data taken from the literature and photometry extrapolations from IRAS, MSX or 2MASS. Among the 343 stars brighter than $`J=5`$ mag in our regions of interest, 11 are not identified with the median filter (Eq.1). However, these remaining contaminated pixels are rejected by the robust linear fitting of DIRBE/2MASS brightness (see Fig. 4).

After removing bright stars, 1040 deg2 of the sky remain.

### 3.3 Faint stars

The 2MASS integrated brightness maps were limited to the formal completeness limits of 2MASS ($`J^{\mathrm{lim}}=15.8`$ mag, $`K_s^{\mathrm{lim}}=14.3`$ mag). A model is required to estimate the contribution of the fainter stellar population with respect to the position in the sky. We use a model in which the number and distribution of Milky Way stars as seen in the near to mid-infrared is adapted from the Bahcall & Soneira (1980) optical star count model. We have employed the discrete population formalism of Elias (1978), Jones et al. (1981), Garwood & Jones (1987) and Jarrett (1992). The model includes the class III (evolved giant), class IV (subdwarfs) and class V (main sequence) stellar populations. These are further divided into disk and spheroid spatial distributions (ala Bahcall & Soneira). The stars are discretely binned according to their spectral types, ranging from the hottest O stars to the coolest M dwarfs, giving a total of 22 separate spectral bins for the main sequence, and 12 for the luminous giants (ranging from G2 to M7 giant). The optical/infrared colors, fluxes and luminosity function per spectral type are based on empirical data (e.g., Koornneef 1983; 2MASS, ISO and IRAS). Interstellar extinction is applied as a smooth exponential function of Galactic position, characterized by a scale height and disk length. The model parameters were tuned using deep optical and infrared star counts. The resultant model was validated using the 2MASS survey down to $`J=15`$ and $`K_s=14`$ mag ($`2`$ mJy).

The result of the model for a typical high galactic latitude field ($`l=210^{}`$, $`b=50^{}`$) is presented in Figure 3. As expected, 2MASS counts are complete at least to $`K_s=14.3`$ mag. According to the model, less than 3% (0.1%) of the total stellar energy is contained in stars fainter than $`K_s=14.3`$ (19). The model is useful for magnitudes ranging from 14.3 to $`20`$ in $`K_s`$. Assuming a conservative number density uncertainty of a factor of 2 for these faint magnitudes, the resulting error will not exceed 3% of the total stellar flux.

## 4 Result

### 4.1 DIRBE / 2MASS correlation

We choose 4 areas with a radius of $`5^{}`$ in order to illustrate the correlation between DIRBE and 2MASS. Pixels contaminated by bright stars have been filtered. Figure 4 show the correlations for $`J`$ and $`K`$ bands. The main characteristics for each region are summarized in Table 1. We note that slopes are not exactly 1 ($`1.07`$ for $`J`$, $`0.94`$ for $`K`$). This arises from the different calibration strategy used by DIRBE and 2MASS: DIRBE is calibrated with Sirius, 2MASS with a list of calibration stars of many different spectral types (Cutri et al., 2000). These calibration differences affect only the slope of the DIRBE/2MASS relation, not the value of the intercept which is the measure of the CIB.

Correlation coefficients show that the correlation is good for both colors (see Table 1). The CIB is obtained by removing the faint star contribution to the intercept value.

### 4.2 Zodiacal contamination

Since the main error is suspected to come from the zodiacal subtraction in the DIRBE maps, a representation of the CIB versus the ecliptic latitude $`\beta `$ is useful. For each pixel, we estimate the intercept of the DIRBE/2MASS linear correlation in a $`5^{}`$ radius area, and we remove the faint star model from the DIRBE data. After this operation we select regions with the best DIRBE/2MASS correlation as defined in Table 2. These values were obtained by examination of the histogram of each quantity, and correspond to the peak value divided by 4. This selection leads to $`550`$ deg2 of reliable regions in the sky.

Figure 5 shows the resulting CIB versus ecliptic latitude for both colors with an additional plot of the DIRBE 25 $`\mu `$m Zodi-subtracted surface brightness. For this figure, we include both low and ecliptic latitude points. The 25 $`\mu `$m band is the most sensitive to the zodiacal light and contains subtraction residuals from the interplanetary dust model. The comparison of the 25 $`\mu `$m plot with the CIB plots confirms that the scatter at low ecliptic latitude comes mainly from residual zodiacal effects.

Systematic uncertainties associated with the zero point of the interplanetary dust models are difficult to estimate. Kelsall et al. (1998) estimate them by choosing the largest among three models at high ecliptic latitude. Accordingly, a conservative estimate of the systematic uncertainties in $`J`$ and $`K`$ is 6.25 kJy sr-1 and 4.4 kJy sr-1, respectively. If the three models used to derive these numbers were equally distributed in the space of possible models, the given uncertainties would correspond to $`1.6\sigma `$. Unfortunately, the models are probably not equally probable and it is hard to interpret these numbers in terms of a confidence level. To be conservative, we assume they correspond to $`1\sigma `$. Theoretical uncertainties in the zodiacal light model are further discussed by Dwek et al. (1998).

Statistical uncertainties that comes from the DIRBE/2MASS correlation and from the zodiacal light can be estimated directly from Figure 5 by measuring the root mean square deviation of the distribution. The resulting uncertainties are 2.7 kJy sr-1 for $`J`$ and 2.1 kJy sr-1 for $`K`$. A summary of the different uncertainty contributions is presented in Table 3. The presence of fluctuations in the residual sky brightness, especially in $`J`$, suggests a small contribution to the sky brightness by starlight scattered off of interstellar dust grains.

### 4.3 Local interstellar medium contamination

Interstellar dust passing through the Solar System will scatter and reemit absorbed sunlight. The volume within which sunlight is significant (compared to the interstellar radiation field \[ISRF\]) is smaller than the size of the Solar System itself: at visible wavelengths, sunlight exceeds the ISRF out to the Oort cloud ($`10^4`$ AU), while at far-ultraviolet wavelengths, sunlight exceeds the ISRF only out to 800 AU (at 0.2 $`\mu `$m) or less. The amount of interstellar dust close to the Sun was recently measured by the Ulysses and Galileo spaceprobes (Grün et al., 1994); the volumetric cross-section is $`n_d\sigma _d2.9\times 10^{23}`$ cm2 cm-3. It was found that the size distribution is deficient in small particles relative to the interstellar size distribution, such that the surface area is dominated by particles with size $`0.4`$ $`\mu `$m.

To estimate the brightness of sunlight scattered by local interstellar dust, we assume the particles detected by Ulysses are spread uniformly throughout the Solar System. The emission from dust in the inner Solar System will have a detectable dependence on solar elongation angle, and would have been included in the DIRBE zodiacal emission model (Kelsall et al., 1998). Therefore, we integrate from 3 AU outward in calculating the brightness of sunlight scattered by interstellar dust and potentially contributing to an isotropic background. The brightness is

$$I_{\mathrm{lism}}=_3^{\mathrm{}}\omega n_d\sigma _dI_{}\mathrm{\Phi }𝑑L,$$

where $`\omega `$ is the albedo, $`I_{}`$ is the solar intensity, and $`\mathrm{\Phi }`$ is the scattering phase function. Using the properties of a mixture of astronomical silicates and graphite (Draine & Lee, 1984), the albedo in the $`J`$ and $`K`$ bands is 0.42 and 0.21, respectively. For the phase function, we used a Henyey-Greenstein function with the asymmetry factor appropriate for the same mixture of astronomical silicates and graphite ($`g_J=0.15`$, $`g_K=0.02`$). The resulting brightness is $`I_{\mathrm{lism}}=0.8`$ and $`0.3`$ kJy sr-1 in the $`J`$ and $`K`$ bands, respectively.

Uncertainty in the estimate of the brightness of the local interstellar dust is due to three sources. First, our lack of accurate knowledge of the composition of the interstellar grains leads to a factor of approximately 2 uncertainty in their albedo; this estimate is based on the difference between the albedos of silicate and graphite grains (Draine & Lee, 1984). Second, the lack of precise distinction between the emission that would be effectively isotropic and that which would be incorporated into the zodiacal light model leads to a $`20`$ % uncertainty; this is based on changing the minimum integration distance from 3 to 2.5 AU. Third, smaller particles that are not detected by Ulysses may be present at larger distances from the Sun; these particles are likely to be produce an insignificant scattering because of the weakness of sunlight at large distances from the Sun.

To summarize, we estimate the CIB as follows: $`\mathrm{CIB}(J)=22.9\pm 7.0`$ kJy sr-1 ($`54.0\pm 16.8`$ nW m-2 sr-1) and $`\mathrm{CIB}(K)=20.4\pm 4.9`$ kJy sr-1 ($`27.8\pm 6.7`$ nW m-2 sr-1).

## 5 Discussion

### 5.1 Comparison with previous studies

Previous estimations of the CIB in $`K`$ are $`16.4\pm 4.4`$ kJy sr-1 and $`14.8\pm 4.1`$ kJy sr-1 by Gorjian et al. (2000) and Wright (2001), respectively. The difference between these values and the estimation given in this paper are within the measurement uncertainties ($`1\sigma `$) and thus are in good agreement. Most of the difference can be attributed to the zodiacal light subtraction because we used the Kelsall et al. (1998) model and they used the Wright (1998) model.

The estimation of the $`J`$ CIB can be compared to the weak limit of $`12.0\pm 6.3`$ kJy sr-1 proposed by Wright (2001). The discrepancy with our value ($`22.9\pm 7.0`$ kJy sr-1) can be explained by the zodiacal light subtraction, since Wright would find 25.8 kJy sr-1 with the Kelsall et al. (1998) dust interplanetary model.

Kashlinsky & Odenwald (2000) have investigated the CIB fluctuations in the DIRBE maps from 1.25 $`\mu `$m to 4.9 $`\mu `$m. Fluctuations on the DIRBE beam scale are expected to be 10-20% of the total CIB flux (Kashlinsky, private communication). Their values of $`6.5_{2.9}^{+1.5}`$ and $`4.3_{2.7}^{+1.2}`$ kJy sr-1 (92% confidence level) for $`J`$ and $`K`$ respectively are therefore in agreement with our results in both bands.

### 5.2 Galaxy counts

Galaxy counts and the CIB are closely related: at least part of the CIB is due to the integrated surface brightness of galaxies for all redshifts. In the past decade, many efforts have been made to obtain galaxy counts in a wide range of magnitudes. Jarrett & et al. (2001) performed large scale studies with 2MASS (2350 deg2) which provide galaxy counts for $`8J15`$ and $`8K_s14`$. Bershady et al. (1998) obtained ground based observations up to $`24^{\mathrm{th}}`$ magnitude in $`J`$ and $`K`$ with the Keck 10m telescope. Saracco et al. (1999) used ESO NTT telescope and obtained galaxy counts for $`17J24`$ and $`16K22`$. Deeper observations are possible with the Hubble Deep Field (HDF) of the Hubble Space Telescope (HST) and Thompson et al. (1999) reached the $`30^{\mathrm{th}}`$ magnitude in the HST filters F160W and F110W which correspond roughly to $`I+J`$ and $`H`$, respectively.

Figure 6 shows the cumulative brightness of galaxies up to $`J=30`$ and $`K=30`$ mag. For the faintest $`J`$ and $`K`$ region (magnitudes 24 to 30) we used the HST measurements assuming $`m_{F160W}m_J=0`$ and $`m_{F110W}m_K=0`$. The resulting error can be neglected, since the contribution of these very faint galaxies to the total brightness is small (9% of the total energy in $`J`$, only 3% in $`K`$). But note that 90% of the energy is in the range $`12J24`$ and $`14K22`$.

Uncertainties on galaxy counts come from the Poisson statistics ($`\sigma =\sqrt{N}`$), and from completeness limitations and corrections. By comparing galaxy counts in the literature, we find that the number of galaxies by range of magnitude $`N`$ is consistent to within a factor of 1.75, which we adopt as a conservative value for the uncertainty. We find a $`J`$ total brightness for galaxies of $`3.6\pm 0.8`$ kJy sr-1, and a $`K`$ total brightness of $`5.3\pm 1.2`$ kJy sr-1, in good agreement with Pozzetti et al. (1998) who obtained 5.8 kJy sr-1 in the $`K`$ band.

The integrated galaxy counts are smaller than our DIRBE-2MASS CIB estimations, even compared with the 95% confidence level. We can address the possibility of large diffuse galaxy halos that would not be included in the galaxy luminosities (because of their faintness) and we find that to reduce the discrepancy between galaxy surface brightness and our CIB value to less than $`1\sigma `$, an increase of each galaxy luminosity by a factor 3.6 and 2.4 is required in $`J`$ and $`K`$, respectively. If we assume that galaxy photometries are correct, a significant part of the energy in the CIB is not detected in sources. The slope of the galaxy luminosity function is observed to flatten for $`K>17`$ (Saracco et al., 1997) and the resulting slope for the cumulative brightness is almost 0 for $`K>25`$ mag. For $`J`$, the slope change in the luminosity function is not as pronounced as for $`K`$, but the slope of the cumulative brightness is also close to 0 for $`J>25`$. Consequently, the slope of galaxy luminosity function must increase for fainter magnitudes, or other contributions to the near-infrared background must be invoked. To summarize, the brightness of the cosmic background found in this paper (Sect. 4) is greater than the integrated brightness of galaxy counts by $`19.3\pm 7.8`$ kJy sr-1 ($`46.3\pm 18.7`$ nW m-2 sr-1) at $`J`$ and $`15.1\pm 6.1`$ kJy sr-1 ($`20.6\pm 8.3`$ nW m-2 sr-1) at $`K`$.

### 5.3 Other possible contributions to the CIB

Although galaxies appear to be unable to explain the brightness of the cosmic background, there are other possible sources of energy production that could lead to an isotropic cosmic background in the near-infrared. An early burst of star formation inside or outside of galaxies, (that is, either in primeval galaxies or in Population III stars), would result in a nearly uniform and isotropic background (Bond et al., 1986). The photospheric emission from these stars would be redshifted into the infrared. For $`z_{}10(100)`$, the background due to high-mass stars will peak at wavelengths 1 (10) $`\mu `$m, and the background due to intermediate-mass stars will be in the range 0.5–2 (6) $`\mu `$m. The gap between our background estimation and galaxy count contributions could be due to Population III stars. The present abundance of heavy elements limits the amount of early star formation such that Population III stars are not likely to be a very significant mass contribution to the universe as a whole (Bond et al., 1986), but the amount of star formation required to explain the near-infrared cosmic background does not violate metallicity constraints. The implications of such early star formation, for the small-scale density fluctuations of primordial matter and for the nature of the intergalactic medium, would be significant.

Depending on the power spectrum of small-scale density fluctuations, a population of black holes could form before galaxies. The accretion of surrounding material by these holes would generate energy that would make an isotropic background today. The temperature of the accretion disk would be $`10^5`$ K, and the accretion luminosity would be redshifted into the near-infrared if the accretion is active at $`z_{bh}30`$–50 (Carr, 1994). Our observed $`J`$ band background brightness, after subtracting the integrated light from galaxies, could be explained by the accretion luminosity of a population of black holes. Such black holes would not comprise enough dark matter to close the universe, and their total mass would not violate constraints on the total baryon density inferred from models of primordial nucleosynthesis. But such a large population of objects (comparable to the total mass of galaxies) would contribute significantly to the total inventory of baryons.

Massive particles that survive until redshift $`z_\nu `$ and subsequently decay into photons will produce a background of decay photons that is subsequently redshifted (Bond et al., 1986). For example, part of the $`J`$ and $`K`$ band background could be produced by the decay of relic 100 eV neutrinos at redshift $`z_\nu 300`$. The number of plausible models is very large, and we can only speculate whether the background is due to such decaying particles.

## 6 Conclusion

The integration of the 2MASS catalog over 550 deg2 (i.e. $`5500`$ DIRBE pixels) leads to an estimate of the CIB for $`J`$ and $`K`$ of $`22.9\pm 7.0`$ kJy sr-1 ($`54.0\pm 16.8`$ nW m-2 sr-1) and $`20.4\pm 4.9`$ kJy sr-1 ($`27.8\pm 6.7`$ nW m-2 sr-1), respectively.

We have selected the most reliable areas of the sky by eliminating pixels contaminated by bright stars, low galactic latitude regions, low ecliptic latitude regions and possible cirrus clouds contribution. A study of DIRBE$``$2MASS versus ecliptic latitude clearly shows that the zodiacal light is responsible of most of the statistical uncertainty ($`2`$ kJy sr-1). The interplanetary dust model contains a systematic error which is actually the dominant uncertainty in our results ($`5`$ kJy sr-1).

The integrated brightness of known galaxies (down to the 30th magnitude) represents only a fraction of the CIB. Although the uncertainties due to residual contribution of the zodiacal light are significant, the discrepancy between DIRBE$``$2MASS and galaxy counts would suggest that a significant part of the energy is still not observed in discrete sources. The slope of the galaxy luminosity functions may increase for galaxies fainter than the 30th magnitude, in which case the counts would agree with the CIB values. However, models predict other possible contributions to the background at these wavelengths (Bond et al., 1986) such as a burst of star formation either in primeval galaxies or in Population III stars ($`z10`$), very massive black holes (accreting from an uniform pre-galactic medium at $`z40`$), massive decaying big bang relic particles ($`z300`$). We have put new constraints in the near-infrared that encourage revisiting the importance of those contributions to the CIB in cosmological models.

The Space Infrared Telescope Facility (SIRTF) may extend this result for the $`L`$ band at 3.5 $`\mu `$m and might be able to constrain the contributions of the different components to the background.

The final all-sky 2MASS catalogs will permit a more accurate study of the near-infrared excess: 100% of the sky will be still available after integration on the DIRBE beam, $`75\%`$ after bright star removal (only 18% for the 2MASS 2nd Release presently used) and about 5000 deg2 will remain at high galactic and ecliptic latitudes. This coverage should allow us to investigate structures in the zodiacal light and may help to improve interplanetary dust models.

We are grateful to Roc Cutri for his help in accessing the 2MASS data and for his critical reading of this paper and to George Helou for helpful discussions.
This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.
The COBE datasets were developed by the NASA Goddard Space Flight Center under the guidance of the COBE Science Working Group and were provided by the NSSDC.
L. Cambrésy acknowledges partial support from the Lavoisier grant of the French Ministry of Foreign Affairs.


format: text

text:
# Monte Carlo simulation of subsurface ordering kinetics in an fcc–alloy model

## 1. Introduction

Near a surface, kinetic processes connected with first–order phase transitions are generally modified in a fundamental manner relative to the corresponding bulk process. Well–known examples are surface–induced spinodal decomposition and nucleation phenomena, whose understanding and control is of major concern in modern thin film and nanostructure technologies. Surface kinetic effects during ordering transitions in metallic binary alloys  acquire special interest from the scientific point of view as their (low–temperature) ordered phase normally is characterized by a multicomponent order parameter, whose components in general couple differently to a crystal surface with a particular symmetry and hence will display different relaxational behaviors. In such cases, the surface is expected to induce rich dynamical behavior of the order parameter, predominantly at short times.

A well–studied system is the (001) surface of $`\mathrm{Cu}_3\mathrm{Au}`$, an fcc–alloy which undergoes a first–order bulk transition from the disordered phase to the ordered $`L1_2`$ structure at a temperature $`T_0=663K`$. In the bulk, $`\mathrm{Au}`$–atoms preferentially occupy one of the four simple–cubic sublattices of the underlying fcc–lattice. The ground state therefore is $`4`$–fold degenerate, and there exist two types of antiphase boundaries separating the $`4`$ types of energetically equivalent equilibrium domains. Several peculiar modifications of order near the (001) surface of $`\mathrm{Cu}_3\mathrm{Au}`$ have been reported. Below the ordering temperature $`T_0`$ the (001)–surface displays disorder–wetting. Moreover, above $`T_0`$, the tendency of Au–atoms to enrich in the outermost (first) layer leads to an oscillatory segregation profile with successive Au–depletion and enrichment in even and odd layers, respectively. This profile decays towards the bulk on a length scale given by the bulk–correlation length $`\xi (T)`$, which on extrapolation to lower temperatures appears to diverge at the spinodal temperature $`T_{sp}T_030\mathrm{K}`$. Intriguing non–equilibrium behavior of segregation amplitudes has been detected recently by time–resolved $`X`$–ray experiments. After a quench from an initial temperature $`T_i>T_0`$ to a final temperature $`T_f<T_0`$ the initial surface segregation profile induces a segregation front which rapidly progresses towards the bulk, while lateral order was found to develop on longer time–scales. This two–stage ordering process was interpreted recently in terms of time–dependent Ginzburg–Landau (GL) theory. Although restricted to one dimension and to the limit of vanishing thermal noise, this theory successfully yields the difference in time–scales for perpendicular and lateral ordering and identifies an anisotropy of the evolving domain structure within a characteristic penetration depth of the segregation front.

Motivated by these studies of the subsurface ordering kinetics of $`\mathrm{Cu}_3\mathrm{Au}`$ (001), our aim here is to gain more insight into possible surface–induced ordering scenarios in $`A_3B`$ fcc–alloys with the help of dynamical Monte Carlo simulations. As a minimal model with respect to equilibrium properties we choose a lattice with nearest–neighbor interactions among $`A`$ and $`B`$–atoms, known to account reasonably for the bulk phase diagram of the $`\mathrm{Cu}\mathrm{Au}`$ system. Our dynamic algorithm relies on the vacancy mechanism, where elementary atomic moves consist in an exchange of $`(A)\mathrm{Cu}`$– or $`(B)\mathrm{Au}`$–atoms with vacancies $`V`$. This kind of $`ABV`$–model, which has been used before by Frontera et al. in an investigation of the bulk ordering kinetics in $`A_3B`$ alloys, is physically more realistic than the conventional direct $`AB`$–exchange kinetics and will simultaneously allow us to scale the simulated surface–induced growth rates with the self–diffusion coefficient of $`A`$– and $`B`$–atoms.

The main results of our simulations are: (i) the occurrence of two distinctly different modes of progression of the segregation front, depending on the final temperature $`T_f`$ to lie significantly above the spinodal temperature, $`T_f>T_{sp}`$ (quench into the metastable regime) or below, $`T_f<T_{sp}`$ (quench into the unstable regime), (ii) an early–time linear (constant velocity) growth behavior of the segregation front in the case $`T_{sp}<T_f<T_0`$ and a characteristic time $`t^{}`$, where this linear growth decelerates. Below $`T_{sp}`$ it changes over to a growth roughly proportional to $`t^{1/4}`$ up to our longest computing times. (iii) A different behavior of lateral growth in the near–surface region. In connection with this we give a detailed description of the persistent anisotropy in the near–surface domain–structure.

The paper is organized as follows. Section 2 defines our model; its bulk properties and equilibrium surface properties are tested in Section 3. The subsequent Section 4 discusses in detail the processes of perpendicular and lateral ordering. Our findings have strong relevance to the temperature quench experiments on $`\mathrm{Cu}_3\mathrm{Au}`$ described above, although the present model based on spatially uniform parameters for nearest–neighbor interactions cannot be expected to work quantitatively. In fact it implies an enhancement of static lateral order near the surface as compared to the bulk, an effect not borne out in the real material $`\mathrm{Cu}_3\mathrm{Au}`$. Enhanced lateral order at the (001)–surface, however, seems to be observable in the fcc–alloy $`\mathrm{Cu}_3\mathrm{Pd}`$. Concluding remarks are presented in Section 5.

## 2. Model and Simulation Method

Consider a lattice of $`N\times N\times M`$ fcc-cells with periodic boundary conditions in the $`x`$– and $`y`$–direction. Along the $`z`$–axis our system is bound between two free surfaces. Atomic layers normal to the $`z`$–direction with distance $`a`$ are labeled by $`n`$ with $`1n2M+1`$. Each lattice site $`i`$ can be occupied either by an $`A`$–atom, or a $`B`$–atom, or a vacancy so that the corresponding occupation numbers $`c_i^A`$, $`c_i^B`$ and $`c_i^V`$ satisfy $`c_i^A+c_i^B+c_i^V=1`$. In our simplified alloy model only pairwise nearest–neighbor interactions are taken into account. The configurational Hamiltonian then reads

$$H=\underset{i,j}{}\left[V_{AA}c_i^Ac_j^A+V_{BB}c_i^Bc_j^B+V_{AB}\left(c_i^Ac_j^B+c_i^Bc_j^A\right)\right],$$
(1)
where the summation is over nearest-neighbor pairs. In view of the very small concentration of vacancies in real alloys  it is natural to assume $`c_i^V1`$ such that macroscopic equilibrium properties of the alloy phases will remain essentially unaffected by the vacancies. For equilibrium considerations one can then ignore the vacancies and introduce spin variables $`s_i=2c_i^A1=\pm 1`$ by which equation (1) is mapped onto the spin–$`1/2`$ Ising model,

$$H_I=J\underset{i,j}{}s_is_jh\underset{i}{}s_ih_1\underset{i}{\overset{}{}}s_i$$
(2)
with

$$J=\frac{1}{4}\left(V_{AA}+V_{BB}2V_{AB}\right)$$
(3)
$$h=3(V_{BB}V_{AA})$$
(4)
The last summation in (2) is restricted to spins in the surface layers $`n=1`$ and $`n=2M+1`$, which are subjected to the surface field

$$h_1=V_{AA}V_{BB}$$
(5)
Bulk equilibrium properties at fixed composition only depend on the parameter $`J`$. Regarding the Cu($`A`$)–Au($`B`$) system, it is natural to take $`V_{BB}>0`$ because of the larger size of $`\mathrm{Au}`$–atoms relative to the $`\mathrm{Cu}`$–atoms, whereas $`V_{AA}<0`$, $`V_{AB}<0`$ such that bulk interactions will favor “anti-ferromagnetic” ordering $`J<0`$. Near the stoichiometric composition $`A_3B`$ the emerging bulk $`L1_2`$–structure is described by a $`4`$–component order parameter, one conserved concentration variable $`\psi _0`$ and three non-conserved components $`\psi _1,\psi _2`$, and $`\psi _3`$. These are defined via the differences $`m_\alpha =c_i^Ac_i^B;i\alpha ;`$ of averaged $`A`$ and $`B`$occupations of the four equivalent simple–cubic sublattices $`\alpha =1,\mathrm{}4`$ building the fcc lattice, namely

$$\begin{array}{cc}\hfill \psi _0& =m_1+m_2+m_3+m_4\hfill \\ \hfill \psi _1& =m_1m_2m_3+m_4\hfill \\ \hfill \psi _2& =m_1m_2+m_3m_4\hfill \\ \hfill \psi _3& =m_1+m_2m_3m_4\hfill \end{array}$$
(6)
Each non-conserved component describes a layerwise modulation of the $`B`$–atom concentration along one of the three cubic axes. The ground state shows $`4`$–fold degeneracy, where the four types of equivalent ordered domains have components $`(\psi _1,\psi _2,\psi _3)/\overline{\psi }=(1,1,1);(1,1,1);(1,1,1)`$ and $`(1,1,1)`$. In the ground state ($`T=0`$) the amplitude $`\overline{\psi }`$ takes the value $`\overline{\psi }=2`$.

The above conditions for the interaction parameters imply a negative surface field $`h_1<0`$, favoring an enrichment of $`B`$–atoms in the outermost layers $`n=1`$ and $`n=2M+1`$. This is to be expected in view of the $`BB`$–repulsion. As shown by Schweika and Landau  on the basis of the Hamiltonian (2) and single–spin flip dynamics, the value $`h_1/|J|=4`$ is consistent with the experimentally measured concentration of $`\mathrm{Au}`$–atoms of about $`50\%`$ in the surface layer of $`\mathrm{Cu}_3\mathrm{Au}(001)`$ near the ordering temperature $`T_0`$.

In our present work, all investigations will be based on the atom–vacancy exchange mechanism, where the energetics are given by the complete hamiltonian (1). For the final choice of interaction parameters see the next section. To affect an elementary atomic move in the Monte Carlo runs we randomly select a nearest–neighbor vacancy–atom pair and calculate the vacancy–atom exchange probability by using the Metropolis algorithm. Throughout most of our calculations the mean vacancy density is taken as $`c^V6.110^5`$, which amounts to $`128`$ vacancies in a system of typical size $`N=64;M=128`$ with $`4N^2M`$ sites. One Monte Carlo step (MCS) consists of $`4N^2M`$ attempted moves and is independent of $`c^V`$. Clearly, in an actual material the vacancy concentration generally will depend on temperature and the amount of $`A`$– and $`B`$–atoms. Our kinetic model with fixed $`c^V`$ therefore does not allow us to compare time–scales of vacancy–mediated processes taking place at different temperatures and composition. However, our subsequent studies of the ordering kinetics focuses on the exact stoichiometric composition $`A_3B`$ ($`c^B=0.25`$) and a narrow temperature range near $`T_0`$ where changes of $`c^V`$ are expected to be of minor importance. Local order parameters to be used in the analysis of our simulations are defined in terms of averages of (6) over one fcc–cell. Order parameter profiles and structure factors are obtained by averaging over 10 independent simulation runs.

## 3. Bulk properties and surface segregation

In this section we discuss some elementary properties of our model, part of which are already known in the literature but they are needed here as a test of our kinetic algorithm and as a basis for the subsequent investigations of surface–induced kinetic properties. Fig. 1a shows the phase diagram in the vicinity of $`c^B0.25`$ as deduced from simulations of a system with size $`M=N=64`$ and periodic boundary conditions in all directions. To determine the ordering temperature at stoichiometry ($`c^B=0.25`$), we analyzed the decay of a perfectly ordered initial state upon thermalization at different temperatures. This yields $`k_BT_01.83|J|`$, in good agreement with literature data ($`k_BT_01.85|J|`$). From a similar procedure, supplemented by computations of concentration ($`\psi _0`$)–histograms, we deduced estimates for the boundaries of the two–phase regions, which reasonably reflect the bulk phase diagram of the $`\mathrm{Cu}\mathrm{Au}`$ system near $`c^B=0.25`$.

Next we proceed to equilibrium properties connected with the ($`001`$)–surface. From our algorithm defined in Section 2 we computed equilibrated $`B`$–atom concentrations in the surface layers as a function of $`h_1`$, keeping $`J`$ fixed. As before, $`h_1`$ and $`J`$ are defined by (5) and (3). Good agreement was found with the work by Schweika and Landau. We consider this both as a positive test of our algorithm and as an indication that the small amount of vacancies used has practically no influence on these static results. In what follows we shall therefore choose the interaction parameters of our $`\mathrm{Cu}_3\mathrm{Au}`$–model such that $`h_1=4|J|`$. While $`J`$ and $`h_1`$ are now determined via the bulk ordering temperature $`T_0`$ and the $`\mathrm{Au}`$–segregation amplitude in the surface layer one more interaction parameter needs to be specified. Our final choice is $`V_{BB}=V_{AA}=V_{AB}>0`$ with $`|J|=V_{BB}/2`$. This can be shown to imply a weak effective attraction between vacancies.

Furthermore we remark that the present model displays a surface transition at a critical temperature $`T_{cs}>T_0`$, which we estimate as $`T_{cs}1.12T_0`$ for $`c^B=0.25`$. Below $`T_{cs}`$ the outermost layers with nearly equal amounts of $`A`$– and $`B`$–atoms develop a lateral (1,1)–superstructure. With respect to $`\mathrm{Cu}_3\mathrm{Au}`$, the existence of a surface critical temperature clearly is an artefact of the present model, connected with the fact that fcc lattice models based on spatially uniform short range interactions are unable to account for surface–induced disordering. The (001)–surface of the $`L1_2`$–structure of $`\mathrm{Cu}_3\mathrm{Pd}`$, on the other hand, appears to favor lateral ordering.

Moreover we estimate the spinodal temperature in a manner as done experimentally. In the disordered phase $`B`$–atoms are known to show an oscillatory segregation profile near the (001) surface, with enrichment and depletion in odd and even layers, respectively. When atoms in the outermost layers are identified with sublattices $`\alpha =3`$ and $`4`$, the segregation amplitude is given by the order parameter component $`\psi _3`$. At the surface, its magnitude is determined by the surface field $`h_1`$. Below the surface it decays towards the bulk on a length scale given by the bulk correlation length $`\xi (T)`$. Extrapolation of data for $`\xi (T)`$ to lower temperatures according to $`\xi (T)(TT_{sp})^{1/2}`$ allows one to estimate the spinodal temperature $`T_{sp}`$. Experimentally, $`T_0663K`$ and $`T_{sp}T_030K`$ near the stoichiometric composition of $`\mathrm{Cu}_3\mathrm{Au}`$. Application of this procedure to our simulation data is illustrated in Fig. 1b and yields the ordering spinodal displayed in Fig. 1a. At $`c^B=0.25`$ we find $`T_{sp}/T_00.967\pm 0.003`$. This is somewhat larger than the experimental value but reasonably close in view of the simplicity of our model. Simultaneously we get (see Fig. 1b) $`\xi (T_0)6a`$, which again agrees fairly well with experiment.

For metallic alloys where short range interactions prevail, it is well known that the spinodal and the nucleation regime connect smoothly. The qualitative significance of the spinodal temperature, however, becomes apparent from Fig. 2a and b. There we compare bulk patterns of the order parameter component $`\psi _1`$ at a time $`t=710^3`$ MCS after a sudden quench from a purely random initial state to final temperatures $`T_f0.972T_0>T_{sp}`$ (Fig. 2a) and $`T_f=0.963T_0\stackrel{<}{}T_{sp}`$ (Fig. 2b), respectively. Snapshots were taken in a section parallel to the ($`xy`$)–plane. In the first case only small ordered domains have nucleated in an essentially disordered background (Fig. 2a). By contrast in Fig. 2b we observe a continuous domain pattern typical of spinodal ordering. At that temperature ($`T_f=0.963T_0`$) the same type of patterns is seen also at shorter observation times. For a complete characterization of the evolving domains in Fig. 2b we need in addition the $`\psi _2`$– and $`\psi _3`$–patterns displayed in Figs. 2c and 2d. From the grey–scales in Figs. 2b–2d one can easily see that inside a domain we have $`\psi _1\psi _2\psi _3<0`$. Furthermore, the well–known two types of antiphase boundaries  are recovered. For example, boundaries where $`\psi _1`$ and $`\psi _3`$ simultaneously change sign in the $`y`$–direction are the “low–energy” (type–1) domain walls and in fact cost no energy in our nearest–neighbor model. These walls are almost completely flat and turn out to be rather stable. On the other hand, boundaries where $`\psi _1`$ and $`\psi _3`$ simultaneously change sign in the $`x`$–direction are “high–energy” (type–2) walls where curvature–driven coarsening is much more effective than for type–1 walls. Bulk coarsening effects at quench temperatures lower than those considered here have been simulated in detail by Frontera et al. These authors obtained coarsening exponents $`n`$ smaller than that in the conventional Lifschitz–Allen–Cahn coarsening law, where $`n=1/2`$, and interpreted their findings in terms of the presence of those stable low–energy walls.

In addition we extract tracer diffusion coefficients from our model, restricting ourselves again to $`c^B=0.25`$ and to temperatures above but close to $`T_0`$. Diffusion coefficients $`D_A`$ and $`D_B`$ for $`A`$– and $`B`$–atoms are obtained in the usual way from time–dependent mean–square displacements and are plotted in Fig. 3 in an Arrhenius representation. Since in our model $`V_{AA}=V_{AB}`$, the energy remains essentially unchanged by a jump of an $`A`$–atom so that $`D_A`$ is found to be nearly temperature–independent. On the other hand, because of the $`B`$$`B`$–repulsion the $`B`$–atoms will migrate more slowly than the $`A`$–atoms, and $`D_B`$ decreases upon cooling. As indicated already, in a real system $`c^V`$ will grow with temperature. Therefore the individual experimental curves $`\mathrm{log}D_A`$ or $`\mathrm{log}D_B`$ versus $`1/T`$ will be steeper than those in Fig. 3. However, the ratio $`D_A/D_B`$ roughly agrees with experimental data for $`\mathrm{Cu}_3\mathrm{Au}`$, where $`D_A/D_B1.45`$ at $`T/T_01.7`$, while we find $`D_A/D_B2`$ at the same $`T/T_0`$–ratio.

Let us close this section by a brief remark on the distribution of vacancies in our system with free surfaces. First of all, vacancies tend to enrich in the outermost layers. For example at $`T/T_0=0.972`$ we find a surface vacancy concentration of about $`c_{\mathrm{surface}}^V2.510^3`$, much larger than the overall concentration. The fact that vacancies are expelled from the bulk suggests that in fact we are dealing with a model for a stable solid which, when regarding the complete ternary ($`ABV`$)–system, would be phase–separated from the “vacuum” (vacancy–rich) phase. Furthermore, when analyzing the above–mentioned layered structure along the $`z`$–axis caused by surface segregation, we find that vacancies deplete in the $`AB`$–layers and enrich in the $`A`$–layers. The reason is that the $`BB`$ repulsive energy in $`AB`$–layers becomes minimized when vacancies avoid $`B`$–sites. The resulting oscillations in the vacancy–density along the $`z`$–axis have an amplitude of about $`30`$ percent relative to the average vacancy density.

## 4. Surface–induced ordering kinetics

Now we turn to the question of how order evolves in a system with free surfaces. Quench conditions are chosen as in Section 3 with final temperatures below $`T_0`$ but in the vicinity of $`T_{sp}`$ ($`0.9\stackrel{<}{}T_f/T_0\stackrel{<}{}1`$). First, within few MC–steps $`B`$–atoms enrich in the surface layer to a concentration of about $`50\%`$ and deplete in the adjacent layer, thereby establishing local equilibrium as enforced by the surface field $`h_1`$. This implies $`\psi _31.85`$ to about $`1.9`$ at $`z=0`$ in the temperature interval considered. This process in turn induces moves of $`B`$–atoms from even to odd layers within an increasing depth from the surface, thus forming a $`\psi _3`$–front which penetrates into the bulk. As an example, we show in Fig. 4 the $`\psi _3`$–distribution in a section parallel to the ($`y,z`$)–plane at $`T_f=0.972T_0>T_{sp}`$ for three different times after the quench. The patterns clearly demonstrate a systematic progression of an ordered state with fixed (positive) sign of $`\psi _3`$. Taking averages over sheets parallel to the surface, we obtain profiles $`\psi _3(z,t)`$ which are plotted in Fig. 5.

At $`t=10^4`$ a plateau region for smaller $`z`$ has evolved with $`\psi _3`$–values in the vicinity of the equilibrium order parameter $`\overline{\psi }1.79`$ at that temperature, before the profile decays to zero for larger $`z`$. Evidently, the thickness of the interface between the partially ordered state and the bulk (not penetrated by the $`\psi _3`$–front) increases with time, an effect expected to arise both from growing intrinsic capillary wave fluctuations  and from bulk fluctuations in all three order parameter components with the character as in Fig. 2a, to be “incorporated” into the front during propagation.

The profiles $`\psi _3(z,t)`$ allow us to define a time–dependent penetration depth $`z_3(t)`$ by setting $`\psi _3(z_3(t),t)=\frac{1}{2}\psi _3(0,t)`$. Plots of $`z_3(t)`$ for a series of final temperatures $`T_f`$ are shown in Fig. 6a for a system large enough in the $`z`$–direction ($`N=64,M=128`$) to avoid the $`\psi _3`$–fronts induced by both surfaces to overlap. (In fact, data in Fig. 5 and 6 contain an average over those two fronts.) As expected, after a quench not passing the phase boundary $`z_3(t)`$ rapidly saturates. For example, at $`T_f=1.006T_0`$ the final penetration depth is about 6 atomic layers, which is of the order of the bulk correlation length $`\xi `$ in the disordered phase at that temperature. On the other hand, when $`T_{sp}<T_f<T_0`$, we see that after a short transient $`z_3(t)`$ grows more or less linearly with time within an extended time interval. This allows us to introduce a penetration velocity $`v(T)`$ which increases upon lowering the temperature. For example, at $`T_f=0.972T_0`$ we find that $`z_3(t)`$ grows linearly up to at least $`10^4`$ MCS, where it reaches about $`70`$ atomic layers, see also Fig. 4. In determining $`v(T)`$, shown in Fig. 7, care has to be taken with respect to the lateral size $`N`$ of the system. For small $`N`$ the system is likely to develop single (closed) domains in the lateral direction within the front region. Increased order of this type will speed up the penetration process. This $`N`$–dependence of $`v`$ is exemplified in Fig. 6b for $`T_f=0.972T_0`$. Data up to $`t=10^4`$ MCS, however, seem to converge sufficiently at our largest $`N=128`$. Moreover, using again $`N=64`$, we have also tested the influence of the vacancy concentration on $`v`$. At $`T_f=0.972T_0`$, results with $`c^V=3.0510^5`$ ($`64`$ vacancies) are indistinguishable from those shown in the figure. Hence the vacancy system appears to be in the dilute limit also with respect to growth processes of this type. Note that for $`\mathrm{Cu}_3\mathrm{Au}`$ the experimental vacancy concentration at $`T_0`$ is even lower.

As seen from Fig. 7, when temperature is decreased from $`T_0`$, the velocity $`v(t)`$ first increases with upward curvature, but decelerates when $`TT_{sp}`$. Sufficiently below $`T_{sp}`$, linear growth can no longer be identified and $`v(T)`$ becomes undefined, see the case $`T_f=0.911T_0`$ in Fig. 6a. Some estimates based on mean–field arguments concerning the order of magnitude of $`v(T)`$ are presented in the Appendix.

When temperature falls near or below $`T_f0.963T_0`$, a new regime of slower growth beyond some crossover time $`t^{}`$ becomes detectable within the time window of Fig. 6. $`t^{}`$ decreases with temperature, so that the regime of linear growth shrinks. Notice in addition the non–monotonous behavior of $`z_3(t)`$ as a function of temperature for fixed observation time $`t`$. As indicated already, the strong shrinkage of the linear regime, where $`t^{}`$ tends to merge with the short time transient, occurs at temperatures slightly below but very near the spinodal temperature $`T_{sp}0.967`$ discussed in Section 3 on the basis of static considerations. Regarding these surface phenomena we conclude that the spinodal temperature as estimated via the static segregation profiles separates two distinct dynamical regimes, roughly characterized by the presence or the absence of surface–induced linear (constant velocity) growth in an intermediate time domain.

The slowing down of the front motion beyond the linear regime, i. e. for $`t>t^{}`$, clearly originates from its competition with nearly saturated domains that have already nucleated in the bulk. Crudely speaking, $`t^{}`$ will be determined by the interplay of two dynamic quantities, the penetration velocity $`v(T)`$ and the nucleation rate of ordered domains. Moreover, unstable growth of bulk domains at temperatures $`T_f<T_{sp}`$ totally suppresses the linear regime. To investigate this situation more closely, we carried out even longer runs up to $`t=210^4`$ for the case $`T_f=0.911T_0`$. Results at that temperature for $`z_3(t)`$ in a log–log representation are shown in Fig. 6c. The data in the range $`t\stackrel{>}{}210^3`$ MCS appear to follow a power–law $`z_3(t)t^n`$ with $`n0.27\pm 0.02`$.

The appearance of such small growth exponent is confirmed by additional simulations at a temperature $`T_f=0.759T_0`$ considerably below $`T_{sp}`$, also shown in Fig. 6c. These data indicate in addition that the onset of this type of power law is moved to earlier times when temperature is lowered. Conversely, data for $`T_f=0.963T_0`$ display a crossover to a slow temporal growth compatible with $`n1/4`$ at a larger crossover time, indicated as $`t^{}`$ in Fig. 6c. Since the penetration depth up to $`t=210^4`$ stays significantly smaller than the system size (in this case $`N=M=128`$), we expect that finite size effects do not play an essential role in estimating that exponent $`n`$. Preliminary results for the structure factor of a bulk system seem to confirm that value of $`n`$, although no definitive conclusions concerning the asymptotic long–time behavior can be drawn from these studies so far. Nevertheless it is interesting to compare our findings in Fig. 6c with recent work by Castán and Lindgård. These authors studied coarsening processes in a $`2`$–dimensional system involving both curved and flat (zero curvature) domain walls. The latter were essentially immobile and could be removed only by the progression of an intervening curved wall. Such a situation basically is encountered also in the present model, where type–1 antiphase boundaries turn out to be essentially immobile. This stability of flat walls induces an overall growth exponent which is smaller than the classical Lifschitz–Allen–Cahn exponent $`n=1/2`$. Castán in fact gave arguments in favor of an exponent $`n=1/4`$. A similar conclusion with respect to $`n`$ was reached by Deymier et al. from studies of a different model with two types of contacts between domains.

Next we study lateral ordering within layers parallel to the ($`xy`$)–plane at a depth $`z`$ from the surface. This process proceeds in different ways depending on whether $`z<z_3(t)`$ or $`z>z_3(t)`$. When analyzing the layer $`z=M`$ midway between the two surfaces, representative for the latter case, we recover the characteristics of bulk ordering as described briefly in Section 3, see also Fig. 2. Our main objective, however, is the near–surface region $`z<z_3(t)`$ and consequences of the surface–induced segregation wave. As emphasized before, for $`z<z_3(t)`$ the sign of $`\psi _3`$ is essentially fixed to $`\psi _3>0`$ so that only two types of domains can appear: $`(\psi _1,\psi _2,\psi _3)/\overline{\psi }=(1,1,1)`$ and $`(1,1,1)`$. These correspond to the two ways $`B`$–atoms can form a $`(1,1)`$–superstructure on the square lattice building an odd ($`AB`$–type) layer. Walls between those domains encountered when going in the lateral direction are exclusively of type 2, whereas type–1 walls are parallel to the surface. To analyze this situation further, we have enlarged the system in the lateral direction to $`N=256`$, but restricted ourselves to $`M=8`$ so that $`\psi _3`$ is nearly uniform along the $`z`$–axis. Patterns of all three order parameter components were extracted from atomic layers $`7`$ and $`8`$ midway between the two surfaces. An example with $`T_f=0.911T_0`$ and $`t=10^4`$ MCS is shown in Fig. 8. In fact, domain walls in Figs. 8a and 8b with simultaneous sign changes of $`\psi _1`$ and $`\psi _2`$ are always curved, indicative of type–$`2`$ walls, and the domain structure appears completely isotropic, in contrast to the anisotropic shape of bulk domains (cf. Fig. 2). Fig. 8c depicts $`\psi _3`$ which is nearly constant.

The time–dependent average size of domains in the ($`xy`$)–section of Fig. 8 is analyzed in Fig. 9, where we present results for the first moment

$$k_\alpha (z,t)=\frac{_k_{}k_{}S_\alpha (k_{},z,t)}{_k_{}S_\alpha (k_{},z,t)}$$
(7)
of the time–dependent lateral structure factor, defined by

$$S_\alpha (k_{},z,t)=|\psi _\alpha (k_{},z,t)|^2$$
(8)
where $`\alpha =1`$ to $`3`$ and $`\psi _\alpha (k_{},z,t)`$ is the lateral Fourier transform of the local order parameter $`\psi _\alpha (\stackrel{}{r},t)`$. Note that vacancy depletion in $`AB`$–type layers and enrichment in $`A`$–type layers (see Section 3) expands the overall time–scale for lateral ordering. While the behavior of $`k_3(t)`$ reflects a rapid transition to a near–uniform $`\psi _3`$–distribution (see Fig. 8c), the data points for $`k_\alpha (t)`$ with $`\alpha =1,2`$ in Fig. 9 seem to approach a growth behavior $`k_\alpha ^1(t)t^n`$ described by the Lifschitz–Allen–Cahn exponent $`n=1/2`$. In fact, our results indicate curvature driven coarsening in the lateral directions, distinctly different from the perpendicular growth displayed in Fig. 6c.

## 5. Summary and concluding remarks

Surface–induced ordering processes in an fcc alloy model were studied on the basis of the atom–vacancy exchange mechanism. One essential result of our simulation studies is the appearance of a segregation wave induced at the surface and propagating into the bulk immediately after a temperature quench across the ordering temperature $`T_0`$. This feature of our model qualitatively agrees with peculiar subsurface ordering phenomena at $`\mathrm{Cu}_3\mathrm{Au}`$ (001) observed in recent $`X`$–ray experiments. In the region $`z<z_3(t)`$ covered by the segregation wave, the present model with nearest–neighbor interactions only shows an increased tendency to build up lateral order. The experimental verification of this effect may depend on the alloy material under consideration.

A second important conclusion confirmed by our simulations is the existence of an anisotropic domain structure in the near–surface region $`z<z_3(t)`$. Only two types of domains and, for $`z`$ fixed, only high–energy domain walls appear, as opposed to the well–known domain structure in the bulk with $`4`$ types of domains and both high– and low–energy walls in either direction. Since in $`\mathrm{Cu}_3\mathrm{Au}`$ $`z_3(t)`$ can become larger than $`200`$Å or about $`50`$ atomic layers (in our simulations more than $`70`$ atomic layers) this should offer the possibility to grow films of mesoscopic thickness with an anisotropic domain structure not realized in the bulk. An interesting question for future work concerns films with odd and even numbers of atomic layers and overlapping $`\psi _3`$–fronts induced by both surfaces.

Moreover, our present studies exemplify different ordering scenarios that may emerge more generally in systems whose bulk structure requires a multicomponent order parameter, and where the surface favors some kind of partial ordering. Concerning domain growth in our model, we have to distinguish between perpendicular and lateral growth modes, between a near–surface partially ordered region ($`z<z_3(t)`$) and the bulk region ($`z>z_3(t)`$), and between quench temperatures $`T_f`$ above and below the spinodal temperature $`T_{sp}`$. At least three distinct coarsening schemes at intermediate or long times were observed. If $`T_{sp}<T_f<T_0`$, perpendicular order initially grows linearly with $`t`$ such that $`z_3(t)vt`$ for $`t<t^{}`$, while for $`t>t^{}`$ our data indicate a temporal regime where $`z_3(t)t^{1/4}`$. On the other hand, for $`T_f<T_{sp}`$ the regime of linear perpendicular growth is suppressed. Analysis of the typical domain structure in the lateral direction shows that in the region $`z<z_3(t)`$ it is isotropic and coarsens according to $`k_\alpha ^1(t)t^{1/2}`$ at long times ($`\alpha =1;2`$) while for $`z>z_3(t)`$ one recovers bulk behavior.

Besides some open questions addressed already above we like to point out that non-stoichiometric alloys should in principle display even more rich and interesting behaviors, in particular in cases where ordering and phase separation occur simultaneously.

## Acknowledgments

The authors are indebted to H. Reichert for several discussions concerning the experimental situation and to P. Maass for helpful comments. A useful conversation with M. Maret, M. Albrecht and G. Schatz about surface segregation effects in the fcc alloy $`\mathrm{CoPt}_3`$ is also gratefully acknowledged. This work was supported in part by the Deutsche Forschungsgemeinschaft, SFB 513.

## Appendix

We attempt here to give some qualitative estimates concerning the magnitude of the penetration velocity $`v(T)`$, see Section 4, from simplified mean–field considerations. Our starting point are time–dependent Ginzburg–Landau equations, based on the free energy density $`f(\psi _1,\psi _2,\psi _3)`$ as given by Lai  for the $`\mathrm{Cu}_3\mathrm{Au}`$ structure,

$$f=\underset{\alpha }{}\left(\frac{r}{2}\psi _\alpha ^2+\frac{v}{4}\psi _\alpha ^4\right)+\frac{u}{4}\left(\underset{\alpha }{}\psi _\alpha ^2\right)^2+w\psi _1\psi _2\psi _3$$
(9)
with $`r=r_0(TT_{sp});r_0>0;u>0;|v|<u;w>0`$. The ordering temperature $`T_0`$ is determined by $`r(T_0)=(2/9)w^2/(3u+v)`$ and the order parameter at $`T_0`$ has the magnitude $`|\psi _\alpha |=\overline{\psi }`$ with $`\overline{\psi }=(2/3)w(3u+v)`$. Regarding the penetration of the segregation wave as a one–dimensional process along the $`z`$–axis (perpendicular to the surface), we assume an equation of motion for $`\psi _3(z,t)`$ of the form

$$_t\psi _3=\mathrm{\Gamma }\left(\frac{f}{\psi _3}K\psi _3^{\prime \prime }\right)$$
(10)
with a kinetic coefficient $`\mathrm{\Gamma }`$ and $`K>0`$. Two limiting cases concerning the coupling of $`\psi _3`$ with $`\psi _1`$ and $`\psi _2`$ can be considered such that (10) becomes a closed equation for $`\psi _3`$ only. Thereby it turns out that the simulated velocity falls between these two limits. First, if the parameters are such that the evolution of lateral order, expressed by $`\psi _1`$ and $`\psi _2`$, considerably lags behind the progression of the $`\psi _3`$–front, one may set in (10) $`\psi _1=\psi _2=0`$. In this case, the form (9) of the free energy density allows a nontrivial solution for a propagating front only at temperatures $`T<T_{sp}`$, with a velocity of propagation $`v(T_{sp}T)^{1/2}`$. Correspondingly, we obtain $`v=0`$ for $`T>T_{sp}`$, contrary to Fig. 7. However, simulations of our nearest–neighbor interaction model indicate that lateral order at least on small length scales evolves almost simultaneously with the front motion in the perpendicular direction. This suggests to consider a second, in fact opposite limit, namely that $`\psi _1`$ and $`\psi _2`$ instantaneously relax towards the “local equilibrium condition” $`\psi _1=\psi _2=\pm \psi _3`$. This yields the approximation $`w\psi _1\psi _2w\psi _3^2`$, to be used in (10). Introducing $`\phi =\psi _3/\overline{\psi }`$ and the “potential”

$`V(\phi )`$ $`=`$ $`f/(3r(T_0)\overline{\psi }^2)=`$ (11)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(1r/r(T_0)\right)\phi ^2{\displaystyle \frac{1}{2}}\phi ^2(\phi 1)^2`$ (12)
we obtain

$$_t\phi =\mathrm{\Gamma }r(T_0)(dV/d\phi +(\xi (T_0))^2\phi ^{\prime \prime })$$
(13)
where $`\xi _0=(K/r(T_0))^{1/2}`$ denotes the correlation length at $`T_0`$. The ansatz $`\phi (z,t)=\eta (zvt)`$ leads to an ordinary differential equation for $`\eta (z)`$. The velocity $`v`$ is determined in the usual manner  by requiring that appropriate boundary conditions can be fulfilled. These are $`\eta (z)\eta _{\mathrm{max}}`$; $`\eta ^{}(z)0`$ as $`z\mathrm{}`$, where $`\eta _{\mathrm{max}}>0`$ corresponds to the absolute maximum of $`V(\eta )`$ and $`\eta (z)0`$; $`\eta ^{}(z)0`$ as $`z\mathrm{}`$. Numerical solution of this problem shows that, at $`T_0`$, $`v(T)`$ starts to grow linearly in $`T_0T`$ and that $`v(T_{sp})1.4\mathrm{\Gamma }r(T_0)\xi _0`$. To estimate $`\mathrm{\Gamma }`$, one can formulate mean–field arguments for the diffusion of the slower atomic species, in this case the $`B`$–atoms, along the $`z`$–axis, which yields $`D_B/a^22\mathrm{\Gamma }r(T_0)`$. Using $`\xi _06a`$, this gives the relation $`v(T_{sp})4(D_B/a)`$, independent of the Monte Carlo time unit. From the simulated value $`D_B0.04a^2/\mathrm{MCS}`$ at $`T_0`$, see Fig. 3, we obtain $`v(T_{sp})0.16a/\mathrm{MCS}`$. Comparison with Fig. 7 shows that these simple arguments overestimate the velocity at $`T_{sp}`$ by about one order of magnitude. This is not surprising in view of the simplifications made, especially in view of the neglect of lateral fluctuations. These may strongly reduce the penetration velocity, an effect already inferred from the $`N`$–dependence of the velocity displayed in Fig. 6b.