text:
# ISOCAM-CVF Spectroscopy of the Circumstellar Environment of Young Stellar ObjectsBased on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA.

## 1 Introduction

The circumstellar environment of young stars is difficult to study, as in its early phases a young accreting and evolving star can be surrounded by sufficient gas and dust to cause a large extinction. However, at visible and near-IR wavelengths this continuum extinction is smooth and contains little information on the nature of the surrounding gas and dust. The re-radiated energy is emitted mainly at wavelengths longer than 2$`\mu `$m. Molecular vibrational-rotational gas-phase and solid-state transitions occur in the mid-IR, and a number of different molecules have been detected towards YSOs (see the review by van Dishoeck & Hogerheijde vDH99 (1999)). Spectroscopy of YSOs has also confirmed the presence of a number of molecular ices, such as ices of water, CO2 and CH3OH (eg. Gibb et al. gibb00 (2000); Boogert et al. boogert00 (2000)). In addition we also see both absorption and emission from refractory solids such as silicates and PAHs, and infrared continuum emission. This continuum emission results from the reprocessing of stellar radiation by the dense dust envelopes of the YSO, and this infrared excess is often used as a criterion for the identification of YSOs in photometric surveys.

Unfortunately, studies from the ground are constrained to work within atmospheric windows. In some cases, such as CO2, telluric absorption by the molecule at 15.2$`\mu `$m makes its study in astronomical objects impossible. Furthermore, the sensitivity of observations in the mid-IR is greatly reduced because telescope and atmospheric thermal emission greatly increase the photon shot noise. Consequently, ground-based studies longward of K band are generally restricted to observing relatively bright objects. The only solution to these problems is to move to an airborne or space platform, where even a modest sized telescope can outperform ground-based systems.

Mid-IR spectra of young stars at high resolution were successfully taken with the Infrared Space Observatory (ISO, see Kessler et al. iso (1996) for an overview of the ISO mission) by various investigators (eg. W33A by Gibb et al. gibb00 (2000); Elias 29 by Boogert et al. boogert00 (2000)). In an attempt to push mid-IR observations to even higher levels of sensitivity we made use of ISOCAM (Cesarsky et al. isocam (1996)) and the Circular Variable Filters (CVF) to image 42 YSO candidates in four well-known low-mass star formation regions: RCrA, $`\rho `$ Ophiuchi, Serpens and Chamaeleon I. The high sensitivity and good sample size was achieved by using the low CVF resolution and by observing dense clusters of YSOs, so that ISOCAM’s imaging capability allowed simultaneous spectra to be obtained. It should be noted that the low spectral resolving power of the CVF ($`\lambda `$/$`\mathrm{\Delta }\lambda 40`$) essentially prohibits the detection of gas-phase spectral features; this is a study of the molecular ices and dust that lie around these young stars. However it is possible to draw a number of conclusions about the physics and chemistry around these YSOs, with information about their composition, thermal behaviour and evolution readily obtainable from such spectroscopic observations.

## 2 Observations and Data Reduction

### 2.1 Observations

5–16.5$`\mu `$m imaging spectroscopy was performed on six target fields in four low-mass star formation regions using the ISOCAM-CVF, on various dates between April 1996 and February 1997. This instrument made use of 2 variable width filters to take images at a number of different wavelengths, using a 32$`\times `$32 pixel gallium-doped silicon array. Exposure times of 2.1s per CVF step were used. Each was exposed repeatedly and the images were sequenced (in time) in decreasing order of wavelength, resulting in a total exposure time of $``$2500s per target region for each CVF set. This procedure produced 2 datacubes for each target region. The first covered the wavelength range 5.0–9.4$`\mu `$m (CVF1) and the other covered 9.3–16.3$`\mu `$m (CVF2), with a mean spectral resolving power of $`\lambda /\mathrm{\Delta }\lambda 40`$. A brief summary of the characteristics of the observed regions is given below.

#### 2.1.1 RCrA

The RCrA cloud is the name given to a molecular cloud of mass $``$120M (Harju et al. harju93 (1993)) around the variable star R Coronae Australis, at a distance of 130pc (Marraco & Rydgren mr81 (1981)). It is a region of low-mass star formation located approximately 18$`\mathrm{°}`$ below the galactic plane (Olofsson et al. olofsson99 (1999)), where it suffers little foreground obscuration. The previous K-band studies of Taylor & Storey (ts84 (1984)) and Wilking et al. (wilking97 (1997)) and the sub-millimetre survey of Harju et al. (harju93 (1993)) revealed a number of embedded objects and other objects classified as YSOs. One field of view was observed in the RCrA cloud, using the 6$`\mathrm{}`$ pixel field of view (pfov). It was centred on 19h01m45$`\stackrel{s}{.}`$1, $``$36$`\mathrm{°}`$58$`\mathrm{}`$34$`\mathrm{}`$111Please note that all coordinates given are J2000.0., immediately to the south-west of RCrA itself.

#### 2.1.2 $`\rho `$ Ophiuchi

The $`\rho `$ Ophiuchi dark cloud complex is a complicated structure of several large molecular clouds near the bright star $ρ$ Ophiuchi. Its relatively high galactic latitude ($`b`$=$`17^{}`$) combined with its relative proximity ($`d`$=160pc, Whittet whittet74 (1974)) result in comparatively little foreground obscuration, making it an ideal region in which to study star formation. Initially identified by Grasdalen et al. (gss73 (1973)), who identified 41 sources in a K-band map, it has since been studied extensively both in the infrared (Elias elias78 (1978); Greene & Young gy92 (1992); Strom et al. sks95 (1995); Abergel et al. abergel96 (1996); Bontemps et al. bontemps01 (2001)) and at millimetre wavelengths (Loren et al. lww90 (1990); André & Montmerle am94 (1994)). In particular, Loren et al. (lww90 (1990)) showed that the cloud consists of a number of smaller sub-structures. They identified 6 dense cores (labelled A-F) and observed significantly higher clustering of YSOs in three of these cores (A, B, and E/F). Two fields were observed in $`\rho `$ Oph, both with the 6$`\mathrm{}`$ pfov: one centred on 16h26m21$`\stackrel{s}{.}`$4, $``$24$`\mathrm{°}`$23$`\mathrm{}`$59$`\mathrm{}`$ (field $`\rho `$ Oph A) and one at 16h27m23$`\stackrel{s}{.}`$7, $``$24$`\mathrm{°}`$40$`\mathrm{}`$39$`\mathrm{}`$ (field $`\rho `$ Oph E). The positions correspond (approximately) to cores A and E/F above.

#### 2.1.3 Serpens

The Serpens molecular cloud was initially identified by Strom et al. (svs76 (1976)). It has since been very well surveyed, both in the infrared (Churchwell & Koornneef ck86 (1986); Gomez de Castro et al. gel88 (1988); Eiroa & Casali ec92 (1992); Giovannetti et al. giov98 (1998); Kaas kaas99 (1999)) and the sub-millimetre (Casali et al. casali93 (1993)), with as many as 163 sources identified in the 6’$`\times `$5’ K-band map of Eiroa & Casali (ec92 (1992)). The cloud itself extends over a region of $``$20’$`\times `$15’ in CO line emission (Loren et al. loren79 (1979)) and is at a distance of $``$260pc (Straižys et al. scb96 (1996)). In addition to survey imaging, many individual objects in the cloud have been studied in detail: Eiroa & Casali (ec89 (1989)) studied the multiple outflow source SVS4; SVS2 has been shown to illuminate the Serpens Reflection Nebula (Gomez de Castro et al. gel88 (1988); Huard et al. hwk97 (1997)); SVS20 has been shown to be a close binary system (Huard et al. hwk97 (1997)). Two fields were observed in Serpens. Field Ser B, taken with the 6$`\mathrm{}`$ pfov, was of the region around SVS2 and SVS20 and was centred on 18h29m57$`\stackrel{s}{.}`$6, 1$`\mathrm{°}`$12$`\mathrm{}`$42$`\mathrm{}`$. Field Ser A was a higher resolution field (3$`\mathrm{}`$ pfov) of the multiple outflow source SVS4, centred on 18h29m57$`\stackrel{s}{.}`$2, 1$`\mathrm{°}`$14$`\mathrm{}`$25$`\mathrm{}`$.

#### 2.1.4 Chamaeleon I

The Chamaeleon dark cloud complex is a complicated structure consisting of 3 large molecular clouds (designated Cha I,II,III by Hoffmeister hoffmeister63 (1963)) and a number of smaller clumps and globules. Again, the relative proximity ($`d`$=160pc, Whittet et al. whittet97 (1997)) and high galactic latitude ($`b`$16$`\mathrm{°}`$) of the clouds result in little foreground obscuration and make them ideal for star formation study. The region has been well-surveyed at millimetre wavelengths (Mattila et al. mlt89 (1989); Henning at al. henning93 (1993)) and more recently in the near- (Kenyon & Gómez kg01 (2001)) and mid-infrared (Persi et al. persi00 (2000)). A single field was observed in the Cha I cloud, with the 6$`\mathrm{}`$ pfov, centred on 11h09m39$`\stackrel{s}{.}`$4, $``$76$`\mathrm{°}`$35$`\mathrm{}`$2$`\mathrm{}`$. The corresponds (approximately) to a previously observed dense core lying to the north of HD97300 (Jones et al. jones85 (1985); Persi et al. persi99 (1999)).

### 2.2 Data Reduction

Initial data reduction was performed using the auto-analysis package from the ISO data pipeline (OLP V10.0, Ott et al. pipeline (1997); Blommaert et al. handbook (2001)). This corrected the images for changes in the spectral response of the detector, as well as performing flatfielding and flux calibration. It should be noted that column 24 on the array was dead: this was corrected for by means of a simple linear interpolation from the adjacent pixels. Sources were identified “by eye” and everything representing a local maximum in intensity in more than 75% of the frames was treated as a source. For each source identified the flux in each frame was evaluated as the total flux contained within a 18$`\mathrm{}`$$`\times `$18$`\mathrm{}`$ (3$`\times `$3 pixel) aperture222A 15$`\mathrm{}`$$`\times `$15$`\mathrm{}`$ (5$`\times `$5 pixel) aperture was used for the 3$`\mathrm{}`$ pfov frames of field Ser A. less the sky flux (evaluated as the mean flux per pixel over a designated sky region). This produced a spectrum for each source from both filter sets. These spectra were smoothed to match the resolution of the CVF and the 2 “halves” taken from the 2 filter sets were then joined together. Unfortunately it was necessary to discard the data from positions 70–80 on CVF1 (8.94–9.43$`\mu `$m), as the array “memory” (a noted problem with ISO, see Blommaert et al. handbook (2001)) rendered these images useless. Consequently the spectra contain a small gap in their wavelength coverage at this point, from 8.94–9.33$`\mu `$m. The spectra were also corrected for a systematic under-estimation of flux at long wavelength, which arose from the use of a fixed-size aperture to measure the total flux from diffraction-limited images.

### 2.3 Spectral Features

An example of a typical spectrum is shown in Fig.1. The following spectral features were identified in the spectra:

* An absorption feature at 6$`\mu `$m attributed to the bending mode of H2O ice (Gerakines et al. gerakines95 (1995)).
* An absorption feature at 6.8$`\mu `$m, the well-known unidentified feature (Schutte schutte97 (1997)).
* A weak absorption feature at $``$7.6$`\mu `$m which was not resolved at the low spectral resolution of the CVF. It is tentatively identified as the bending mode of CH4 at 7.67$`\mu `$m (Hudgins et al. hudgins93 (1993); Boogert et al. boogert96 (1996); Boogert et al. boogert98 (1998)) but due to the low resolution of the CVF and the weak nature of the observed feature this identification is somewhat uncertain.
* A broad feature centred at $``$9.7$`\mu `$m, seen in both emission and absorption, identified as the Si–O stretching mode of “astronomical silicate” (Draine & Lee dl84 (1984); Simpson simpson91 (1991); O’Donnell odonnell94 (1994)).
* A broad absorption feature at $``$13$`\mu `$m identified as the libration (or hindered rotation) mode of H2O ice (Gerakines at al. gerakines95 (1995)).
* An absorption feature at 15.2$`\mu `$m attributed to the bending mode of CO2 ice (Ehrenfreund et al. ehren96 (1996)).

Other spectral features were observed in individual sources, such as the apparent 10.5$`\mu `$m emission feature in RCrA IRS5 (see Fig.2) or the apparent 14$`\mu `$m absorption feature in HH100-IR (see Fig.4). These rarer features have been noted in the results tables where they occur, but the goal of this study was to look for general trends rather than individual curiosities. With this in mind, and also considering the large number of spectra obtained, it was decided to fit only the strongest and most commonly occuring features highlighted above.

### 2.4 Fitting the Spectra

As many of these spectral features overlap, independent determinations of their strengths were not possible, so all of the features in the spectra were fitted simultaneously. This was achieved by fitting a spectrum of the form:

$$F_\nu (\lambda )=C(\lambda )\times \mathrm{exp}\left(\underset{i}{}s_i\alpha _i(\lambda )\right)$$
(1)
where $`C(\lambda )`$ is the continuum, $`s_i`$ represents the absorbing column density of each species $`i`$ and $`\alpha _i(\lambda )`$ are the absorption coefficients for each species $`i`$. This essentially treats the absorbing material as a cold absorbing screen. In order to account for possible silicate emission, the depth of the silicate feature, $`s_{\mathrm{Si}}`$, was allowed to run negative. Note however that this assumes that the same profile applies in both absorption and emission, so could result in errors if both are present. The absorption coefficients for H2O ice and CO2 ice were taken from the Leiden Observatory Database of Interstellar Ices (Gerakines et al. gerakines95 (1995); Ehrenfreund et al. ehren96 (1996)) and those coefficients evaluated at 10K were used, as thermal variations are negligible at the spectral resolution of the CVF. It should be noted, however, that the depths of the 6$`\mu `$m and 13$`\mu `$m water ice bands did not correlate as expected, so they were fitted independently. The unidentified feature at 6.8$`\mu `$m was fitted using an asymmetric Gaussian profile: a Gaussian profile with two different widths shortward and longward of the peak. The best-fitting profile for each 6.8$`\mu `$m feature was measured, and the widths were then fixed as the means of the widths measured in all of the observed profiles (FWHM=500nm). The poorly resolved 7.6$`\mu `$m feature was fitted using a Gaussian profile with a width equal to that of the CVF resolution element at that wavelength. A number of profiles exist in the literature for the silicate feature, obtained both theoretically and observationally, so the profile which provided the best fit to the observed data was used: profile 2 from Simpson (simpson91 (1991)). However the profiles obtained by Simpson (simpson91 (1991)) from IRAS observations extend shortward only to $``$7.5$`\mu `$m, and so shortward of this wavelength the optical constants for “circumstellar silicates” from Ossenkopf et al. (ohm92 (1992)) were used. It should be noted that the profiles from Simpson (simpson91 (1991)), Ossenkopf et al. (ohm92 (1992)) and Draine (draine85 (1985)) are almost identical away from the 9.7$`\mu `$m peak: only the shape of the absorption band is dependent on the choice of optical constants.

The fitting procedure was iterative, with repeated fitting converging to a final profile. The fitting algorithm took the form:

1. Fit initial continuum and divide out. The continuum was fitted as a 2nd-order polynomial through points at 5.6, 7.9 and 16.1$`\mu `$m, so chosen as they best constrained the resulting fits.
2. Measure the strength of each spectral feature at a number of points around the peak.
3. Use these strengths and the optical constants for the different species to divide out the features from the spectra.
4. Fit a new continuum to this “featureless” spectrum.
5. Iterate this procedure until the fit converged: typically 10-15 iterations were sufficient.

It should be noted that this allowed the continuum to drift away from the observed spectrum in some cases, due to absorption across the entire observed wavelength range. This is because the normalised silicate opacity over the range 5-7.5$`\mu `$m is approximately constant at $``$0.15, and towards the long wavelength end of the spectra the 18$`\mu `$m bending feature begins to “cut in”. (The normalised opacity at 16$`\mu `$m is approximately 0.3.) The effect of this can be seen in Fig.2: although the continua are seen to deviate significantly from the observed spectra, appearing to diverge from the observed spectra in some cases, this is primarily due to silicate absorption across the entire wavelength range. With only the silicate features included, these “continua plus silicate” fits are seen to follow the observed spectra well. As a result, any possible systematic errors arising from the manner of the continuum fit will only affect the evaluation of the silicate optical depth; the measured strengths of the narrower features are not affected by the choice of continuum.

On inspection of the residuals from this fitting procedure a further feature was observed: a broad feature centred on $``$11.2$`\mu `$m, seen in both emission and absorption. This has a number of possible identifications, such as crystalline silicates (Bregman et al. bregman87 (1987); Campins & Ryan cr89 (1989)), the Unidentified (PAH) Band at 11.3$`\mu `$m (Boulanger et al. boulanger96 (1996)), or a shoulder on the silicate profile. Given the variation in possible identifications, it was decided to fit this feature independently using an asymmetric Gaussian profile, once again taking the widths to be the mean of those observed (FWHM=1700nm). This feature was then included in the algorithm and the fitting procedure repeated: examples of successful fits are shown in Fig.2. Another similar residual around 8.5$`\mu `$m was also observed, but the properties and strength of this were greatly affected by the choice of continuum point around 8$`\mu `$m. Whereas the 11$`\mu `$m feature was robust against small changes in the continuum this apparent 8.5$`\mu `$m feature was not, so it was not included in the fitting procedure.

The 15.2$`\mu `$m CO2 ice feature is an interesting one and has been used in the past as a diagnostic of the ice environment. Gerakines at al. (gerakines99 (1999)) found that two distinct phases of CO2 ice (polar and non-polar) show significantly different absorption profiles. The non-polar phase is characterised by a double-peaked profile, arising from the two-fold degeneracy of the bending mode, whereas the polar (H2O-rich) phase shows a significant long wavelength wing. The spectral resolution of the CVF is too low to resolve the double-peak structure, but the CO2 profiles were observed to vary noticeably from source to source. As a result it was decided to let the CO2 profile vary, so it was fitted using an asymmetric Gaussian profile with the widths as variable parameters.

There are two forms of error associated with this fitting procedure. Firstly there are the ordinary statistical errors associated with the fits, evaluated by combining the statistical errors on the data points used to fit the features and the intrinsic uncertainties in the individual feature fits. These are typically around 5–10% for the stronger features (silicate, H2O, CO2) and around 15–20% for the less well constrained features, although in the case of some poor S/N spectra they are much larger. In addition to this, there are systematic errors associated with poor fits. These are usually due to the presence of “extra” features in the spectra, or because of poorly fitting feature profiles, and are less easy to quantify. Consequently, the errors reported in the results tables are the statistical errors only; where fits are flagged as poor additional systematic errors also apply.

### 2.5 Dealing with poor fits

Generally the procedure described above proved successful, but in some cases it led to poor fitting of certain features. This was primarily due to the manner in which the continuum was fitted, as shown in the examples in Fig.3. The continuum was fitted as a 2nd-order polynomial through three sections on the spectrum expected to show the least absorption, as described above. However, constraining the continuum through these three points led in some cases to poor fitting elsewhere, particularly at the short end of the CO2 ice band and immediately shortward of the silicate feature.

In such cases, the narrower features at the short end of the spectrum (6.0, 6.8 and 7.6$`\mu `$m) were fitted in exactly the same way as before, only this time using a linear continuum through the points at 5.6 and 7.9$`\mu `$m. As absorption due to the silicates is approximately constant over this wavelength range this removed complications due to the quadratic continuum fit. When tested on spectra which were fitted well, such as those in Fig.2, this produced results which were the same as before, within the fitting errors. Also, the manner in which the continuum was fitted often led to the CO2 feature fitting poorly, as seen in Fig.3. As a result it was decided to measure the CO2 absorption in all the sources by fitting a linear continuum from 14.5–16.1$`\mu `$m, as the silicate absorption does not vary much over this range and the CO2 is not blended with any other spectral feature. This procedure was only applied to these narrower features, however, as the broad silicate, 11 and 13$`\mu `$m bands were deemed to be too heavily blended to deconvolve them meaningfully over a short section of the spectra. In addition, all features fitted in this manner are flagged as such in the results tables.

### 2.6 Spectral Index

In order to take a standard measurement of the SED, the spectral index $`\alpha `$ was evaluated as a traditional IR spectral index (eg. Kenyon & Hartmann kh95 (1995)):

$$\alpha =\frac{\mathrm{log}\lambda _2F_{\lambda _2}\mathrm{log}\lambda _1F_{\lambda _1}}{\mathrm{log}\lambda _2\mathrm{log}\lambda _1}$$
(2)
This sets $`\alpha =0`$ as the index of a flat $`\nu F_\nu `$ and $`\alpha =3`$ as the index of a Rayleigh-Jeans spectrum. Equivalently, class I YSOs have $`\alpha >0`$, class II have $`3<\alpha <0`$ and class III $`\alpha 3`$, using the conventional definition of spectral class (Lada lada87 (1987)), although more recent studies (eg. André & Montmerle am94 (1994); Greene et al. greene94 (1994)) have found that the class II/III border lies at a higher value ($`\alpha 1.5`$). Where previous K-band data exist in the literature the wavelengths were taken as $`\lambda _2=14.0`$$`\mu `$m (our data) and $`\lambda _1=2.2`$$`\mu `$m (K-band data from references in Table 1) and this was applied to both the observed spectrum and the calculated continuum for each source. (Where no K-band data exist, for sources not previously observed or with uncertain identifications, $`\lambda _1=8.0`$$`\mu `$m was used.) This second (“continuum”) measurement of the spectral index is important, as many of the sources suffer significant absorption (including foreground absorption) and a straight application to the observed spectrum could result in sources being mis-classified in cases where deep absorption is present. In cases where the flux suffers little extinction the two give the same results, but where deep absorption is present the continuum spectral index derived is greater than that obtained from the observed spectrum. In such cases the continuum spectral index provides a truer measurement of the shape of the underlying continuum.

## 3 Results

### 3.1 Detected Sources

A total of 43 sources were detected: a complete list is presented in Table 1. The 6$`\mathrm{}`$ pfov used here results in coordinates which are accurate only to $`\pm 3`$$`4\mathrm{}`$ at best, so where sources have been identified with previous studies the coordinates presented are those from previous observations, where higher spatial resolutions or better S/N have resulted in more accurate positions. The coordinates evaluated here are only presented for sources which represent new or uncertain detections. In the 3 cases where sources have been identified as unresolved doubles the coordinates presented are those of the brighter source. Of the 43 sources detected 6 were not identified with detections from previous surveys (see Sect. 3.1.1) and 2 (Ser A 7 and Ser B 11) were determined to be different images of the same source, as the two Serpens fields overlap slightly. Consequently we were able to obtain complete 5.0–16.3$`\mu `$m spectra for a total of 40 sources, along with a further 2 partial spectra333These partial spectra were the result of various observational problems, such as sources lying near to the edge of the frame or lying on the bad column.: a complete atlas of these 42 spectra is presented in Appendix A.

#### 3.1.1 Ghost Images and New Detections

As noted by Blommaert et al. (handbook (2001)) ghost images and stray-light are a significant problem when using the CVF, so care must be taken when identifying sources. Ghost images arise in two ways and are found either around the true image or directly opposite it (relative to the optical axis). Potential ghost images were identified by looking at the spectra of the objects, their flux densities and their positions relative to the instrument optics. In this manner 7 objects initially identified as sources were reclassified as ghost images (4 in $`\rho `$ Oph E, 2 in Cha I and 1 in RCrA). This leaves 6 “new” detections in our survey, which we now address:

* 3 sources in RCrA are identified as new detections: sources RCrA 2, 6 and 7. Source RCrA 2 lies close to RCrA 1, separated by $``$30$`\mathrm{}`$, but shows sufficiently different spectral characteristics that it cannot be a ghost of this source. It is therefore considered to be a real source, probably a YSO. Both RCrA 6 and 7 lie in a region of nebulosity between the bright sources HH100-IR (RCrA 5) and RCrA (which lies just outside the image frame). These are not considered to be ghost images, but may well represent peaks in the nebulous emission rather than true YSOs. It should also be noted that this region of the RCrA cloud was not included in the ISOCAM survey of Olofsson et al. (olofsson99 (1999)), due to detector saturation, so it is not unreasonable to expect new detections here.
* A large number of ghost images were found here and after they were removed 1 unidentified source remained: $`\rho `$ Oph E 4. However, this is the faintest source in the survey, with a flux rising to just 0.17Jy at 16$`\mu `$m, so is probably a peak in the nebulosity rather than an “new” YSO.
* 2 sources in Ser B were not identified with previous observations: Ser B 4 and 7. Both these sources lie near to SVS20 and so could be ring-like ghost images associated with it. However they also lie in a previously identified region of nebulous emission (Kaas kaas99 (1999)) and so may instead be peaks in this nebulosity. It is impossible to resolve this ambiguity completely, but in either case these 2 sources are probably not “new” YSOs.

#### 3.1.2 Background Stars

Of the 36 sources identified with previous studies 3 have previously been identified as stars which are not associated with the star-forming clouds: CK2 (Ser B 12) by Chiar et al. (chiar94 (1994)) and Casali & Eiroa (ce96 (1996)); Ser B 1 and Ser B 11 / Ser A 7 by Giovannetti et al. (giov98 (1998)).

CK2 is the only one of these 3 sources which has been identified as a background field star (probably a background supergiant, Casali & Eiroa ce96 (1996)). Unfortunately it fell directly on the bad column in the array here and so no useful spectroscopic data regarding CK2 was obtained. Both of the other two sources in Serpens show extremely deep absorption features due to ices and silicates and both have very red SEDs. Taking the K-band magnitudes of Giovannetti et al. (giov98 (1998)) and assuming the SEDs of these stars to be Rayleigh-Jeans spectra ($`F_\nu \lambda ^2`$) results in predicted mid-IR flux densities of $``$0.1mJy for both sources, approximately 1000 times less than what is observed. Consequently it seems unlikely that these are background sources, so we interpret them to be deeply embedded objects: with this interpretation these 2 sources are 2 of the 3 most deeply embedded objects in this survey.

### 3.2 Spectral Fitting

The spectra of the 43 sources identified in Table 1 were fitted using the method described above: the results of this procedure are presented in Table 2 and the fits and continua obtained are presented in Appendix A. The fitting procedure is both robust and unambiguous, converging to a good fit in most cases, but a few weaknesses exist. Firstly, as discussed in Sect. 2.3, the procedure only measures the 7 spectral features observed in almost all of the sources and consequently does not fit some rarer features, such as an apparent 14$`\mu `$m absorption band observed in source $`\rho `$ Oph A 2. Further, the profile fitted to each feature (except the CO2) is not allowed to vary from source to source. While this is broadly valid, a few exceptions led to some poor fits, such as the broadened 6.0 and 6.8$`\mu `$m bands in SVS2 (see Fig.3), or a handful of sources which show broadened silicate features. However, given the low spectral resolution of the CVF and the variation in the spectra over the large number of sources

detected, this was found to be the most reliable and consistent method of spectral fitting.

In order to check the instrumental calibration and the validity of this fitting method, comparisons were made with previous observations of the bright, well-studied object HH100-IR (RCrA 5, see Fig.4). Whittet et al. (whittet96 (1996)) find a silicate optical depth of 1.21$`\pm `$0.05, which is somewhat less than the value of 1.35$`\pm `$0.05 measured here: this is almost certainly due to the fact that the iterative fitting procedure used here allows the continuum to drift away from the observed spectrum slightly. As discussed in Sect.2.4, this is due to absorption across the entire wavelength range and is a benefit made possible by the broad wavelength coverage of the CVF. Whittet et al. (whittet96 (1996)) also find a H2O ice column density of 2.4$`\times `$1018cm-2, based on observations of the 3$`\mu `$m stretching mode. Adopting a band strength of $`A`$=1.2$`\times `$10-17cm molecule-1 (Gerakines et al. gerakines95 (1995)), we evaluate the column density $`N`$ as:

$$N=\frac{\tau (\lambda )d\lambda }{\lambda _{\mathrm{peak}}^2A}$$
(3)
where $`\tau (\lambda )d\lambda `$ is the equivalent width. This gives a value of 4.8$`\times `$1018cm-2, approximately double the value obtained from the stretching mode. This is consistent with the “6$`\mu `$m to 3$`\mu `$m paradox” (Gibb et al. gibb00 (2000); Dartois & d’Hendecourt dd01 (2001)), which commonly results in observations of these two bands producing column densities which differ by a factor of $``$2. The recent work of Gibb & Whittet (gw02 (2002)) suggests that the excess depth of the 6$`\mu `$m feature is the result of blending with a feature due to organic refractory matter. They also find a strong correlation between the excess absorption in the 6$`\mu `$m water ice feature and that of the 4.62$`\mu `$m “XCN” feature; unfortunately the wavelength range of the CVF prevented observation of the “XCN” feature here. Similarly, from observations of the 4.27$`\mu `$m stretching mode Nummelin et al. (nummelin01 (2001)) infer a CO2 ice column density of 6.2$`\pm `$0.6$`\times `$1017cm-2, which is broadly consistent with the value of 5.1$`\times `$1017cm-2 measured here (assuming a band strength of $`A`$=1.1$`\times `$10-17cm molecule-1, Gerakines et al. gerakines95 (1995)). Keane et al. (keane01 (2001)) measure optical depths of $`\tau _{6.0}`$=0.23 and $`\tau _{6.8}`$=0.09, which are consistent with our results of 0.23$`\pm `$0.01 and 0.13$`\pm `$0.01 respectively. It should be noted, however, that our value of the 8$`\mu `$m flux density, 1.05$`\pm `$0.02$`\times `$10-16W cm-2 $`\mu `$m-1, is somewhat less than the value of $``$1.7$`\times `$10-16W cm-2 $`\mu `$m-1 measured by Whittet et al. (whittet96 (1996)). This may indicate a possible error in the absolute flux calibration of the CVF data. However, the measurement of spectral features depends only on the relative flux calibration from pixel to pixel, which is considered to be accurate throughout.

#### 3.2.1 Water Ice: Bending and Libration Modes

The observed equivalent widths of the 6$`\mu `$m bending and 13$`\mu `$m libration modes of water ice showed no significant correlation at all. Further, the libration mode was not observed to correlate with any of the measured features. Whilst there is evidence for this in the literature (eg. Bowey at al. baw98 (1998)) the most likely explanation is that the libration mode is poorly fitted, due to the wide variation in possible profiles. The peak of this band has previously been found at $``$11$`\mu `$m in crystalline water, $``$12.5$`\mu `$m in amorphous water and at even longer wavelengths in mixtures with other molecules (Hagen et al. htg83 (1983); d’Hendecourt & Allamandola dha86 (1986); Cox cox89 (1989)). Here a single profile was used (that of amorphous ice) and this probably resulted in poor fitting. Unfortunately this section of the spectrum is strongly blended with the silicate feature, so this problem cannot be remedied without making further ad hoc assumptions about the silicate profile(s).

#### 3.2.2 The Unidentified Feature at 6.8$`\mu `$m

As seen in Fig.5, the strength of the unidentified 6.8$`\mu `$m feature is observed to correlate far more strongly with the neighbouring 6$`\mu `$m band of water ice than with the 15.2$`\mu `$m band of CO2 ice. This is consistent with previous observations, where the 6 and 6.8$`\mu `$m features usually correlate well (see Schutte et al. schutte96 (1996) and the review by Schutte schutte97 (1997)). This suggests that the carrier of the unidentified feature is probably a strongly polar ice, as its presence matches far more closely the strongly polar H2O ice than the markedly less polar CO2. The more recent work of Keane et al. (keane01 (2001)) found the peak position of this feature at different wavelengths towards different sources, and proposed that it consists of two (related) components. Further detailed study of the 6.8$`\mu `$m profiles could yield a great deal more information about their chemistry. However, at the low spectral resolution of the CVF the 6.8$`\mu `$m profiles are essentially constant across all the sources, and so such an investigation was not possible here. One individual source is also worthy of note: GY262 ($`\rho `$ Oph E 3) shows a deep 6.8$`\mu `$m feature with no corresponding 6.0$`\mu `$m feature (the “rogue” point to the upper left in Fig.5).

#### 3.2.3 11.2$`\mu `$m Feature

The depth of the measured 11.2$`\mu `$m feature was found to correlate strongly with the depth of the silicate feature, with a direct (negative) proportionality providing a good fit to the data. Consequently the 11.2$`\mu `$m feature was interpreted to be an emissive shoulder on the silicate feature, rather than an independent feature due to another species: the presence of this feature narrows the silicate absorption profile slightly. The silicate profile has previously been found to vary depending on the composition and structure of the silicate grains (eg. Demyk et al. demyk00 (2000)), so such a shoulder is not unexpected. However, 3 of the sources (RCrA 1, Ser A 6 and Cha I 3) show a significant 11.2$`\mu `$m emission feature without any significant silicate absorption feature: these features may be attributable to emission from crystalline silicates (Bregman et al. bregman87 (1987); Campins & Ryan cr89 (1989)).

#### 3.2.4 CO2 Ice Profiles

As mentioned in Sect. 2.4, the CO2 ice profile is a very sensitive diagnostic of the ice environment. Whilst the spectral resolution of the CVF was not sufficient to study the profiles in great detail, it is worth noting that a significant long wavelength wing was present in almost all the observed absorption profiles. This would seem to indicate that a large fraction of the CO2 ice observed exists in a polar (H2O-rich) phase, which is characterised by this long wavelength wing (Gerakines et al. gerakines99 (1999)). However quantifying the relative abundances of the different phases was not possible at this low spectral resolution.

#### 3.2.5 Spectral Classification and Silicate Profiles

If we look at the shapes of the spectra, it appears natural to divide the sources into 3 distinct groups based on the strengths of the spectral features observed. An example of each type is shown in Fig.6 and the classifications are included in Table 2.

* The first group show deep absorption features due to all of the ices and also due to the silicates (typically these have $`\tau _{\mathrm{Si}}0.6`$). These objects are interpreted to be heavily embedded objects, showing strong absorption due to cold foreground material.
* The second group shows weaker ice absorption features and weak silicate features ($`0.8\tau _{\mathrm{Si}}0.6`$) which are, in general, not especially well fitted. These are interpreted as less heavily embedded objects, showing less foreground absorption. The poor silicate fits are probably due to complications caused by combined silicate emission and absorption.
* The third group show strong silicate emission ($`\tau _{\mathrm{Si}}1.0`$) and little or no ice absorption. It should be noted that in the 2 good S/N spectra showing strong silicate emission, the observed emission profiles were significantly broader than both the observed absorption profiles and the profile used to fit the data (see Fig.6). (However it should also be noted that one of these 2 sources, SVS2, is an unresolved double.) These are probably young stars which have shed most of their circumstellar matter, retaining only a little hot, optically thin silicate dust around them (possibly in a disc). It therefore seems plausible that these three groups might represent an approximate evolutionary sequence, with the YSOs shedding circumstellar material as they evolve.

In order to test this the spectrum of SVS20 (Ser B 6) was re-fitted using a composite silicate profile, consisting of the standard absorption profile used above superimposed on the broader emission profile of SVS2 (Ser B 3 - see Fig.6). As seen in Fig.7 this provided a good fit to the observed silicate profile, indicating that both silicate emission and absorption probably do occur, with different profiles seen in absorption and emission. The physical interpretation of this is that these objects have a hot region of optically thin silicates around them but remain embedded in the cloud, which results in foreground absorption due to silicates and ices being superimposed on the “intrinsic” silicate emission feature in the spectra.

It should be noted that while it is possible to fit such profiles to these sources, they are not well constrained. As can be seen by comparing the 2 composite profiles in Fig.7 there is a degeneracy between the strength of this “flat-topped” silicate profile and the continuum strength, and a great deal of information is assumed about the silicate profiles themselves. Consequently this investigation was not pursued for all the sources. Where $`0.8\tau _{\mathrm{S}i}0.6`$ the values in Table 2 represent the correct optical depth at 9.7$`\mu `$m, but are usually in error towards the wings of the silicate feature. Such values are useful in studies of general trends, but they do not accurately describe the entire silicate feature: a single number cannot fully describe the complicated nature of these features.

This also provides more direct evidence, similar to that of Whittet et al. (whittet88 (1988)), that the extensive scatter in measurements of the $`A_V/\tau _{\mathrm{Si}}`$ ratio towards YSOs may be caused by complications due to silicate emission. Direct measurements of the silicate optical depth will significantly under-estimate the depth of foreground material in cases (such as SVS20) where optically thin silicate emission is also present. Rieke & Lebofsky (rl85 (1985)) derive $`A_V/\tau _{\mathrm{Si}}17`$, whilst noting the problem posed by “intrinsic” silicate emission, and such emission does indeed lead to errors when evaluating $`A_V`$ from the silicate depth. As an example, the ISOCAM survey of Bontemps et al. (bontemps01 (2001)) derives values of $`A_V28`$ for GY252 ($`\rho `$ Oph E 1) and $`24`$ for GY262 ($`\rho `$ Oph E 3); the X-ray survey of Imanishi et al. (ikt01 (2001)) derives slightly lower values of $`24`$ and $`18`$ respectively. However GY252, as seen in Fig.6, clearly shows a “composite” silicate profile whereas GY262 does not, and the respective silicate optical depths of 0.52 and 1.84 are clearly inconsistent with a single $`A_V/\tau _{\mathrm{Si}}`$ ratio. Unless the emitting and absorbing components of such silicate features can be separated unambiguously obtaining a single $`A_V/\tau _{\mathrm{Si}}`$ ratio will not be possible.

#### 3.2.6 Spectral Index

Of the 41 objects to which spectral indices were assigned 20 were found to be class I objects ($`\alpha _{\mathrm{cont}}>0`$) and 17 to be class II ($`1.5<\alpha _{\mathrm{cont}}<0`$). In addition to this 3 objects ($`\rho `$ Oph E 1 and 3, Ser B 2) were found to lie very close to the class I/II border (so-called transition objects: André & Montmerle am94 (1994); Greene et al. greene94 (1994)). The lack of any class III objects is a selection effect, but it seems to indicate that all of the observed sources are within the clouds, as field stars would show little or no infrared excess and would appear as class III objects in this survey. The distribution of class I and II objects is approximately constant throughout all four regions, and no significant trends involving spectral class were found. Evaluation of mid-IR spectral indices appears to support the finding of Bontemps et al. (bontemps01 (2001)) that mid-IR spectral indices do not discriminate strongly between class I and class II objects, and a classification scheme similar to that in Sect. 3.2.5 seems preferable when considering spectroscopic data.

Comparing the index-based spectral classes to the spectral feature grouping above, we see that the majority of the class I sources belong to group a) and the majority of the class II sources to group b). However it is notable that a significant number (6) of the 23 group a) objects are either class II or transition objects. This in keeping with the convention that class I objects are the youngest (Lada lada87 (1987)), but it may be significant that there are two distinct types of spectrum which both result in class II SEDs. Objects showing strong foreground absorption and objects showing little foreground absorption can both present class II continuum SEDs. This could be due to geometry of the YSOs relative to the absorbing clouds, or could be an effect intrinsic to the sources: spectroscopy along a pencil-beam cannot distinguish these. Whilst a larger sample size would constrain this problem further, it does seem that mid-IR observations alone do not provide strong constraints on the “true” spectral class of YSOs.

## 4 Variations Between Target Regions

In general, the most striking regional variation is between the sources in Cha I and those in the other three regions. All of the sources in Cha I, except for ISO-ChaI 192 (Cha I 2), show little or no absorption due to volatile ices and show silicates either in emission or an emission/absorption composite. As it lay outside the CVF1 images, due to rotation of the spacecraft between the two sets of observations, we only have a 9.3–16.3$`\mu `$m spectrum for ISO-ChaI 192. However it shows CO2 ice absorption and also appears to show deep silicate absorption ($`\tau _{\mathrm{Si}}2`$–4), so it is probably an embedded object. Persi et al. (persi99 (1999)) found that ISO-ChaI 192 lies in a dense core and that it is probably the engine for the bipolar outflow detected by Mattila et al. (mlt89 (1989)). Jones at al. (jones85 (1985)) found that several other YSOs have formed around the edge of this core, rather than at its centre: these were apparently formed during a burst of star formation triggered by a wind from the nearby star HD97300. Our results are consistent with this: ISO-ChaI 192 appears to be heavily embedded, and has a very red SED. The rest of the surrounding objects show similar spectral characteristics, and appear to be surrounded by far less circumstellar/foreground material than ISO-ChaI 192.

12 of the 20 objects in Serpens show deep absorption features, 7 show composite silicate profiles, and one (SVS2) shows silicate emission. It is clear that Serpens contains a wide range of YSOs, ranging from the most heavily embedded object observed (Ser B 11) to a source showing silicate emission. By contrast almost all of the sources in both RCrA and $`\rho `$ Oph show deep absorption features, with only 1 of the 7 sources in RCrA and only 2 of the 10 sources in $`\rho `$ Oph showing composite silicate profiles. Given that YSOs are generally expected to “sweep out” their circumstellar material as they evolve (Shu et al. shu87 (1987)), this seems to imply that the observed star formation is at a similar evolutionary stage in both RCrA and $`\rho `$ Oph, probably an earlier stage than that in Serpens. However, given the broad variation in the sources observed in Serpens, this conclusion remains somewhat tentative. It should also be noted that the observed fields (of 2$`\mathrm{}`$$`\times `$2$`\mathrm{}`$) are much smaller than the star-forming clouds, which typically extend over several square degrees of the sky, so the populations observed here may not be representative of the entire clouds.

### 4.1 Silicate Depth and H2O Bending Mode

Figure 8 shows the measured equivalent widths of the 6$`\mu `$m water ice feature plotted against the optical depth of the silicate feature: clear regional variations are apparent. As discussed above, all but one of the sources in Cha I are characterised by little or no ice absorption and by silicate emission; the other 3 regions show absorption from both silicates and ices with only a few examples of weak silicate emission. The sources in RCrA are characterised by similar absorption strengths throughout, whereas far deeper absorption and variation is observed in both $`\rho `$ Oph and Serpens. It is clear from Fig.8, however, that the sources in $`\rho `$ Oph show markedly less water ice absorption relative to silicates than those in Serpens. This could be due to sublimation effects, indicating a difference in the thermal structure and history of the two clouds, or could represent a real difference in the chemical composition of the two clouds; without measuring the corresponding gas-phase lines we cannot say which. Sub-millimetre observations of gas-phase H2O (Ashby et al. ashby00 (2000)) do show an under-abundance of gaseous H2O in $`\rho `$ Oph A relative to S140, but to the best of our knowledge no similar direct comparisons between the sources observed here exist in the literature. It is also worth noting that the silicate:water ice ratios appear to be approximately constant for all the different embedded sources in both $`\rho `$ Oph and Serpens. This indicates that in both clouds the composition of the absorbing intracloud material is roughly constant along different lines of sight, a result consistent with the previous Serpens work of Eiroa & Hodapp (eh89 (1989)).

### 4.2 H2O and CO2 Bending Modes

Figure 9 shows the ratio of the CO2:H2O ice equivalent widths (which is proportional to the corresponding ratio of column densities) plotted against the silicate optical depth. The results from Serpens and RCrA are fairly well correlated, with CO2:H2O ice column density ratios ranging from 0 to 0.16 (adopting the band strengths from Gerakines et al. gerakines95 (1995)). There is a possible trend showing the CO2:H2O increasing with decreasing silicate optical depth, but this is somewhat unclear in these data.

The data from $`\rho `$ Oph, however, show a far greater scatter with CO2:H2O ice column density ratios ranging from 0 to 0.4. This is consistent with previous observations: for example Chiar et al. (chiar94 (1994), chiar95 (1995)) found a greater scatter in column densities of CO and H2O ices against $`A_V`$ in $`\rho `$ Oph than in Serpens, RCrA or Taurus. The reasons for this are unclear but it seems likely that, as suggested by Chiar et al. (chiar95 (1995)), local conditions around the young stars in $`\rho `$ Oph play a role. Detailed study of the gas-phase abundances of CO2 and H2O towards these objects will further constrain this problem, but without such observations or further knowledge of the envelope structures we merely note the presence of this discrepancy here.

## 5 Conclusions

ISO-CVF spectroscopy of four low-mass star-formation regions has been performed, and a total of 42 different sources were observed. The resulting low-resolution 5–16.5$`\mu `$m spectra have been analysed, and a number of different conclusions drawn:

* The profile of the 9.7$`\mu `$m silicate absorption feature varies somewhat from source to source, but usually shows a weak emissive shoulder at $``$11.2$`\mu `$m. Significantly broader profiles are observed in emission. A number of sources appear to show “composite” profiles, with absorbing and emitting components superimposed.
* The presence of a long wavelength wing on the 15.2$`\mu `$m CO2 ice band in most of the sources suggests that a significant fraction of the CO2 ice observed exists in a polar (H2O-rich) phase.
* The strength of the unidentified 6.8$`\mu `$m band is observed to correlate far more strongly with the 6.0$`\mu `$m water ice band than with the 15.2$`\mu `$m CO2 ice band. This suggests, in a manner consistent with previous observations, that the carrier of this band is a strongly polar ice.
* Compared to the other 3 regions (which show similar spectral characteristics) the sources observed in Cha I are somewhat anomalous, showing comparatively little foreground absorption. However this is consistent with previous studies, which indicate that the observed Cha I sources represent an unusual cluster of YSOs which have formed at the edge of a dense cloud core.
* The $`\rho `$ Oph cloud appears to have a systematic under-abundance of water ice, relative to silicates, relative to both RCrA and Serpens.
* The CO2:H2O ice ratios towards the sources in $`\rho `$ Oph show considerably greater scatter than in the other regions, with values ranging from 0 to 0.4. This may be due to local conditions around the YSOs in $`\rho `$ Oph but more data is needed to constrain this problem further.

The most important limiting factor in this study is obviously the low spectral resolution of the CVF, which essentially prevents any detailed study of absorption profiles, particularly around the CO2 ice feature. However the use of the CVF enabled the study of a large number of objects, many of which would not have been observable at higher spectral resolution due to the low S/N. Consequently we have obtained a unique data set, with a large sample size and broad wavelength coverage. This enabled simultaneous study of the dominant mid-IR spectral features along the line-of-sight to a large number of YSOs, in a manner not possible using other spectroscopic methods. Further detailed study of the ice absorption profiles will undoubtedly yield more information regarding the chemistry around these objects, and study of the silicate profiles may yield more information about the dust composition and the structure of the dust envelopes.

###### Acknowledgements.

We thank Tom Greene for his work in planning the observations and for helpful comments. We thank Andy Longmore for useful comments about the manuscript and also thank Sylvain Bontemps for useful discussions. Part of this work formed part of the MPhys dissertation of RDA at the University of Edinburgh, and RDA also thanks the UKATC for funding a vacation studentship. Finally, we thank an anonymous referee for helpful advice which greatly improved the clarity of the paper.

## Appendix A Atlas of Spectra

All of the spectra are presented in this appendix. In each plot the data points are shown as crosses, with the fits as solid lines and the fitted continua as broken lines. The spectra are arranged in the same order as in Table 1. The error bars shown are the statistical errors on the data points and do not include the absolute flux calibration errors discussed in Sect. 3.2. Please also note that the fits shown are the “complete” fits, covering the complete wavelength range: in some cases (as indicated in Table 2 and discussed in Sect. 2.5) the narrower features were fitted independently and these fits are not shown here.


format: text

text:
# Kinetic Monte Carlo Study of Electrochemical Growth in a Hetheroepitaxial System

## 1 Introduction

The study of the kinetics of formation and growth of two dimensional phases in electrochemical systems is a classical problem in the field of physical chemistry. This sort of studies are experimentally performed by applying various perturbations to the system in order to obtain the response of a macroscopic observable which allows to infer the distinctive features of the nucleation and growth process. Over the last two decades, the improvement in in-situ nanoscopy techniques has also made possible to observe the time evolution of the surface morphology in electrodeposition experiments as it is done in the case of deposition experiments under ultra-high vacuum (UHV) conditions.

From a theoretical point of view, a large variety of models are generally used to describe nucleation and crystal growth phenomena . Many of these models apply continuum theories, while others use descriptions based on microscopic concepts, as does the DDA (Diffusion, Deposition, Aggregation) model which will be referred to in section 4. Most of these models have been formulated to describe experimental situations similar to the ones set up under UHV conditions, where the adsorption process occurs randomly and at constant flux at any place on the surface. However, during electrodeposition of a metal, the entrance of ad-atoms to the system is a thermally-activated process in which the activation energy for ion reduction is determined mainly by the reorganization energy of the solvent . This particular feature has raised a question about the type of site on the electrode surface where particle deposition occurs. It is possible that in some systems the reduction reaction occurs mainly on the edge of previously formed islands, whereas in other systems particles deposit preferentially onto the clean terraces and then diffuse in order to join the growing islands as proposed by Bockris . A conclusive answer to this question has not yet been found. This present work has two main aims. On one hand, we seek to study the most remarkable features of the electrochemical response and the evolution of surface morphology of the system for the two limiting cases mentioned above, in comparison with the trivial case where deposition occurs at the same rate on all free places of the surface. On the other hand, we shall give simple heuristic criteria that should allow the identification of the nature of the discharge process from experimental data.

## 2 Approach to the problem

As a model system, we consider here electrodeposition in the case of a negligible adsorbate/substrate lattice misfit through the dynamic Monte Carlo (DMC) method. In order to describe the interaction between particles in a realistic metallic system, we use semiempirical potentials corresponding to Ag deposition on Au(100). However, we believe that many features of the present model will also be valid for other hetheroepitaxy systems beyond the particular potential employed to describe the interaction between the particles in the system. Furthermore, the DMC method allows to study the system evolution in real time. Concerning the type of perturbation applied to the system, in the present work we shall simulate a potentiostatic experiment, were the electrode potential is changed from an initial value where the surface of the electrode is free from adatoms, to a value where adsorption occurs. Furthermore, we shall consider irreversible deposition, that is, we neglect adatom dettachment from the surface.

The relationship between the atom deposition rate on a given type of site $`x`$ and the potential applied is given by the following equation:

$$k_{dep,x}(E)=k_{dep,x}^0a_{Me^{z+}}exp\left(\frac{\mathrm{\Delta }G_{dep,x}^{(0)}}{RT}\right)exp\left(\frac{(1\alpha _c)zFE}{RT}\right)$$
(1)
where $`k_{dep,x}`$ is the atom deposition frequency on site $`x`$ , $`E`$ is the measured-electrode potential vs. a reference electrode, $`\mathrm{\Delta }G_{dep,x}^{(0)}`$ is the activation energy for ion transfer from the solution to the crystal at E=0, $`\alpha _c`$ is the charge transfer coefficient, $`a_{Me^{z+}}`$ is the activity of the metallic ion in the electrolyte, and $`k_{dep,x}^0`$ is the rate constant for the deposition reaction on a site $`x`$. According to equation (1), a simulation of adatom deposition appears a priori as a formidable task, since $`k_{dep,x}(E)`$ should be calculated for every possible type of site on the surface. This would require a detailed knowledge of the electron transfer itself. For example, the calculation of $`k_{dep,x}(E)`$ could involve on one hand the local electronic properties of the site where the discharge of the ion is taking place, as well as the nature of the rearrangement of the solvent concomitant with the charge transfer process. However, as we stated in the introduction, in this work we shall be concerned with some special cases of equation (1) and we shall think of the surface as essentially made of two types of surface sites: those which are next to at least one adatom and those which are not. Thus, in principle only two different values of $`k_{dep,x}(E)`$ may occur. Among the various possibilities that emerge even after this simplification, the following three possible mechanisms for atom attachment to the growing phase were studied, which are considered in Figure 1:

1. Particle adsorption occurs on all unoccupied sites on the surface, which are considered as equivalent, independently of whether or not the site is surrounded by atoms. In this case, a rate $`k`$ is assigned to the process of atom entrance into an unoccupied site, independently of its surroundings.
2. Particle adsorption occurs only at sites corresponding to the terraces, i.e. where no adsorbed atoms are around. In this case an adsorption rate $`k`$ is assigned to those sites. The adsorption rate is considered to be $`0`$ at those sites that have at least one of the nearest neighboring sites occupied, i.e. on the edge of steps or on kink sites. This model is based on Bockris’ idea that the ion reduces preferentially on the terraces, since it loses there the least part of its solvation sphere. After that, the ion diffuses towards the step edges and then towards the kink sites.
3. Particles are considered to be discharged preferably at step edges. With this purpose, a rate $`k_1`$ is assigned to the entrance into terrace sites, i.e. where no neighboring atoms are around, and a rate $`k_2`$ (where $`k_2`$ $`>`$ $`k_1`$) is given to the entrance of atoms into sites that have at least one atom as its nearest neighbor. This situation is comparable to the Avrami model  in which the monolayer grows from the island edges, assuming that atoms enter more easily there.

In the present model the parameter $`k`$ was changed between $`10^3s^1`$ and $`10^{+2}s^1`$ (at each adsorption site) and the nature of the site $`x`$ was determined by its atomic enviroment as stated above. In the present model no assumption is made concerning the evolution of the surface morphology. This is a result of the deposition rate $`k`$ and the interaction between the particles of the system that affects the diffusion properties. Thus emphasis is set on making a realistic model of atom diffusion on the surface, according to a potential that takes into account the many body effects and which is suitable for the study of metals.

This paper is organized as follows: the model and the simulation technique used are described in Section 3; the results are given in Section 4, which is subdivided into two parts: in the first, the system dynamic response is studied in terms of the evolution of the coverage degree and of the current (defined as $`i=d\theta /dt`$) as a function of time; in the second part, structural aspects are dealt with, that is, a description of the morphology and quantification of the formed islands is given. Finally, our conclusions are presented in Section 5.

## 3 Lattice Model and Simulation Technique

### 3.1 Dynamic Monte Carlo Method

Monte Carlo methods are used as computational tools in many areas of physical chemistry. Although traditionally applied to obtain equilibrium properties, they can also be used to study dynamic phenomena . In order to do this, the following conditions must be fulfilled: a) in addition to satisfying the detailed balance criterion, the probabilities of transition must reflect a ”dynamic hierarchy”; b) time increments between events must be correctly formulated in terms of the microscopic kinetics of the system; c) the events must be effectively independent. In the DMC method every step consists in a random selection of one of the possible processes. The probability of a process being selected is directly proportional to its rate. Once the randomly chosen method has been performed, all the possible processes are calculated and stored again in the corresponding vector and there is a time increment of $`\mathrm{\Delta }t=ln(u)/v_i`$ where u is a random number between 0 and 1, and $`v_i`$ is the sum of the rates of all the possible processes. This time increment is due to the assumption that we are dealing with a Poisson process.

Silver adsorption on a defect-free Au(100) surface was studied. The computational model worked with a square arrangement of $`n\times n`$ adsorption sites, under periodic boundary conditions. The simulations were performed with $`n=50`$ and $`n=100`$. Unless otherwise stated, the results presented correspond to the former case. The sites corresponding to the nearest neighbors on each side were kept in a matrix. The system was initialized: for the present studies, the initial state corresponded to a clean Au(100) surface, but the model could be easily extended to surfaces exhibiting previously adsorbed silver or gold atoms. The rates of every possible process were stored in a vector to perform the random selection described above. In this work since we consider the case of irreversible deposition (no desorption), such processes were the entrance of an atom into an empty site and the motion of an atom from a site to one of the four sites corresponding to first nearest neighbors (face 100). For example a fixed rate may be assigned to the first process, independently of the occupation of the neighboring sites. This would cause an increase in the coverage degree given by a law of the type $`\theta =1exp(kt)`$, where t denotes the time elapsed. In the most general case, the entrance rate of a particle depends on whether or not the site is surrounded by other atoms, as discussed below.

The algorithm developed by Hoshen and Kopelman was used to calculate the island number and size  and average was taken over 10 simulations for each rate. The interval was divided into 100 bins in order to calculate the average at each coverage degree.

### 3.2 Energy Calculation

To calculate the activation energies for adatom diffusion the Embedded Atom Method (EAM) was used . This method takes into account many body effects; therefore, it represents better the metallic bonding than a pair potential does. The total energy of the system is calculated as the sum of the energies of the individual particles. Each energy is in turn the sum of an embedding (attractive) energy and a repulsive contribution that arises from the interaction between nuclei. The EAM contains parameters which were fitted to reproduce experimental data such as elastic constants, enthalpies of binary alloys dissolution, lattice constants, and sublimation heats. We have employed the EAM within a lattice model as described in references .

### 3.3 Diffusion Rate Calculations

Diffusion rates used in the simulation were previously tabulated according to the different environments that an atom can find on the surface. Various configurations were generated, consisting of different arrangements of atoms near the starting and ending sites. For each configuration the path followed by an atom to jump from a site to the neighboring one was traced and the energy in each position was calculated, minimizing it with respect to z (coordinate perpendicular to the surface plane). The activation energy $`E_a`$ was calculated as the difference between the saddle point and the initial minimum in the energy curve along the reaction coordinate. The vibrational frequency $`\nu `$ of the atom in the starting site was calculated performing the harmonic approximation near the minimum of the curve. The diffusion rate was calculated as $`v=\nu exp(E_a/kT)`$. More details of this calculation procedure can be found in reference

## 4 Results

### 4.1 Evolution of the coverage degree and the current as a function of time

#### 4.1.1 Model A: Adsorption Independent of Surroundings. Dynamic response.

Model A is the simplest one and the evolution of the coverage degree responds to the equation d$`\theta `$/dt=$`k`$(1-$`\theta `$), where k denotes the adsorption rate per site. In this case $`\theta `$=$`1exp(kt)`$ and the current satisfies the law $`i=kexp(kt)`$. In a plot of the coverage degree as a function of $`kt`$, the curves corresponding to the different simulations overlap, i.e. the evolution of $`\theta `$ as a function of $`kt`$ is independent of the atom entrance rate. However, as we shall later see, the surface structure differs significantly regarding the island number and size. Note the difference of the present electrochemical (potentiostatic)conditions with respect to the UHV ones. In the latter case, the coverage degree evolves linearly with time, while in the potentiostatic case the rate is proportional to the available surface resulting in an exponential law.

#### 4.1.2 Model B: Preferential Adsorption on Terraces. Dynamic response.

For case B, that is, with atom adsorption only on the terraces, Figure 2 shows the evolution of $`\theta `$ as a function of $`kt`$ for several $`k`$ values. Here it can be observed that for higher entrance rates the surface tends to a lower coverage, i.e. it takes longer to reach a high coverage, while for lower rates the coverage tends to be closer to one. This is so because at higher rates more and smaller islands are formed, generating a larger amount of edge sites where deposition is not allowed. On the other hand, at low deposition rates the particles have enough time to diffuse reaching a more compact distribution, i.e. forming fewer and larger islands, therefore leaving more sites available for adsorption. Further analysis shows that at relatively short times $`\theta `$ initially evolves like $`\theta =1exp(\beta t^\alpha )`$, with $`\alpha `$$`<`$1. This can be verified by a plot where ln(-ln(1-$`\theta `$)) is represented as a function of ln(t) (Figure 3), since straight lines of slope $`\alpha `$ and y-intercept $`ln\beta `$ are obtained. Table 1 shows $`\alpha `$ and $`\beta `$ values for the various k. For higher rates we found $`\alpha `$ values considerably lower than those in case A (equivalent sites) where $`\alpha `$= 1. For lower rates $`\alpha `$ increases tending to 1. If the system is allowed a longer time to evolve, a change in the exponent is observed, i.e. the curve becomes much more horizontal. This effect is much more pronounced at high adsorption rates. The physical explanation for this is that at a given moment many small islands have already formed, and consequently few unoccupied sites are left on the terrace, which makes the entrance of new atoms more difficult. At this point diffusion starts playing a fundamental role, since the islands start rearranging in order to form larger ones. This process is slow compared with the diffusion of free atoms on the terrace. As diffusion occurs, adsorption sites are generated, allowing a slow entrance of new atoms, which produces a small current. Finally, at high coverage vacancy islands and simple hollows are left . There the diffusion rate of vacancies is higher and so is the generation of adsorption sites, which causes the curve slope to increase at the final stage.

Figure 4 shows the evolution of $`\theta `$ and of the current as a function of time for k = 1 s<sup>-1</sup>. Here a comparison can be made between model A (adsorption equal at all sites), and model B (adsorption only on the terraces). In the case of the latter, the coverage degree is found to be lower and tends to 1 more slowly for the reasons stated before. For the former case, current decreases exponentially, whereas for adsorption only on the terraces, it decreases faster at the beginning and more slowly afterwards, satisfying the law $`d\theta /dt=\beta \alpha t^{\alpha 1}exp(\beta t^\alpha )`$. The points obtained as the result of taking average over 10 simulations are shown as well as the curves predicted by this ansatz, generated with the values of $`\alpha `$ and $`\beta `$ obtained from the fit.

#### 4.1.3 Model C: Preferential adsorption on the Edges. Dynamic response.

The third case corresponds to atom adsorption that is preferential on step edges. Simulations with different adsorption rates on the terraces ($`k_1`$) and on the edges ($`k_2`$) were performed. Table 2 shows the rate constants of all the systems that have been studied here.

In these simulations the time at which the first particle deposits is taken as the initial time, since the time elapsed until this first process occurs may be, in some cases, significantly longer than the time taken by the simulation, especially when particle adsorption occurs on the edges of islands at rates much higher than those on the terraces. Figure 5 shows the evolution of the coverage degree as a function of time for the cases having $`k_2=1.0s^1`$ as entrance rate on the step edges, and $`k_1=0.01s^1`$, $`k_1=0.001s^1`$ and $`k_1=0.0001s^1`$ as entrance rates on the terrace. In these three cases, as well as in all cases where the adsorption rate is higher on step edges than on the terrace, it can be noticed that the initial evolution of the coverage degree has a curvature that first points upwards and then downwards, i.e. it presents a saddle point. This is due to the fact that the current depends mainly on the amount of edge sites, which initially increases and afterwards decreases as coalescence occurs.

In these cases, as before, the coverage degree was assumed to satisfy a law $`\theta =1exp(\beta t^\alpha )`$ and was represented in a plot as $`ln(ln(1\theta ))`$ as a function of $`ln(tt_0)`$ in order to obtain the coefficients $`\alpha `$ and $`\beta `$. These are summarized in Table 3. These values give a qualitative idea of the tendencies, since they correspond to a single simulation. It can be observed that, in all cases, $`\alpha `$ is larger than 1, and it tends to decrease for a given value of $`k_1`$ as $`k_2`$ decreases. Figure 6 shows the current as a function of time, along with the curves $`d\theta /dt=\beta \alpha t^{\alpha 1}exp(\beta t^\alpha )`$, with the parameters obtained from the fit. Three examples are shown. The plot at the bottom of the page corresponds to $`k_1=10^4s^1`$ and $`k_2=10^2s^1`$. Here the two rates are very different from each other and it can be clearly observed that initially the current is 0, then it increases up to a peak , and after that decreases. The plot in the middle corresponds to an intermediate case, where the rates are $`k_1=10^3s^1`$ and $`k_2=1.0s^1`$ and the curve is qualitatively very similar to that of the previous case. The plot at the top corresponds to $`k_1=10^2s^1`$ and $`k_2=2\times 10^2s^1`$. In this case the two rates are very similar and a current maximum is not observed in the simulated curves. This case is more similar qualitatively to the exponential decrease corresponding to equal adsorption at all sites.

### 4.2 Structural Aspects

In this section we shall mention first some of the quantities often used to describe the dynamic behavior of the surface structure in the so-called atomistic crystal growth models. We shall especially make use of the concepts used in the DDA model, which considers three basic processes as responsible for epitaxial growth: Deposition, Diffusion, and Aggregation . Then, we shall study the behavior of these quantities in each of the three cases considered in this work and their relationship with the corresponding current transients described in the previous section.

#### 4.2.1 Use of the Island-Size Distribution Function to Describe Surface Morphology

A fundamental quantity in describing kinetically the growth of epitaxial monolayers and submonolayers is the island-size distribution function $`N_s(t)`$, which is the density per site of islands of size $`s`$ at time $`t`$, $`s`$ being the number of atoms in the island. The function $`N_s(t)`$ is used in the DDA model to describe the surface structure and it has mainly been applied to the analysis of experimental results obtained under UHV conditions. In this case the coverage degree $`\theta `$ changes linearly with time and then the function of size distribution is generally expressed as a function of $`\theta `$. In this work we shall follow the same criterion, though we should keep in mind that in our simulations $`\theta `$ does not evolve trivially, because of which many of the relationships found do not strictly correspond to those of the dynamic scaling.

Defining the total number of islands $`N`$ and the coverage degree $`\theta `$ for a given instant by

$$N=\underset{s2}{}N_s\theta =\underset{s1}{}sN_s$$
(2)
then the average island size $`S(\theta )`$ can be written as a function of the zeroth and first moments ($`_{s2}N_s`$, $`_{s2}sN_s`$) as

$$S(\theta )=\frac{_{s2}sN_s}{_{s2}N_s}=\frac{\theta N_1}{N}$$
(3)
In analyzing the evolution of the surface structure and its relationship to the evolution of the current and of the coverage degree, we shall distinguish various dynamic regimes in the behavior of the total island density $`N`$ and of the monomer density $`N_1`$. Then, we shall present our work in a way similar to the one by Family et al .

It is important to highlight that the scaling theory is satisfied only under certain deposition conditions (e.g. relatively low deposition rates), and naturally, only during the aggregation regime. Thus only some of the performed simulations can be expected to satisfy completely the observations predicted by this theory. Such is the case for models A and B when the coverage degree is approximately at the interval (0.1,0.4) and the deposition rate is low (low overpotentials), as we shall see in sections 4.2.2 and 4.2.3.

#### 4.2.2 Model A: Adsorption Independent of Surroundings. Structural Aspects

As mentioned in Section 4.1.1, in model A the evolution of the coverage degree and of the current are given trivially by an exponential function. However, the observed variety of surface morphology depending on the deposition conditions indicates that no relationship exists between the surface structure and the electrochemical response.

Six different entrance rates of atoms were studied: $`10^{+2}s^1`$, $`10^{+1}s^1`$, $`10^0s^1`$, $`10^1s^1`$, $`10^2s^1`$ and $`10^3s^1`$. Figure 7 shows the evolution of the surface for two different adsorption rates: the left-hand column corresponds to $`k=10^{+2}s^1`$, and the right-hand column to $`k=10^2s^1`$. From top to bottom the coverage degrees are $`\theta =0.1`$, $`\theta =0.3`$ and $`\theta =0.6`$. As it can be observed, at higher adsorption rates more and smaller islands are formed. This is true when coverage is low or intermediate, since at high coverage and rates islands are branched so they coalesce in order to form one island. At high rates and low coverages a large number of monomers can be observed. This is due to the fact that, when particles have such a high entrance rate, they have little time to diffuse, whereas at a low entrance rate particles have time to join the already existing islands.

Figure 8 shows the total island number, and the monomer number divided into the total site number for the various entrance rates of atoms. The upper and the lower part of the figure corresponds to the island density and to the monomer desity respectively. As it can be observed, at higher deposition rates both the monomer density and the island density are larger. In all cases the monomer maximum is reached at lower coverages than the maximum in the island density. At low deposition rates three zones can clearly be distinguished during deposition. The first one consists in an initial increase in the island number, which corresponds to the nucleation process. Nuclei are formed at the beginning, so we can refer to an instantaneous nucleation. Then, the number of islands remains essentially constant and this corresponds to the aggregation regime in which every adatom joins an existing island whose size increases in this way, but whose number does not. Finally, the coalescence regime sets in, during which the number of islands decreases as they begin joining one another to form larger islands. Eventually, only one island is formed which grows up to form the complete monolayer. At higher deposition rates, the island density reaches a maximum at higher coverage degrees, and we can refer to a progressive nucleation, since the number of islands is increasing over a larger zone of $`\theta `$ and it is not possible to differentiate a specific aggregation stage.

Figure 9 shows the average size of islands (calculated as the average number of atoms of every island divided into the total number of sites) as a function of the coverage degree for the various adsorption rates. In all cases there is, initially, a linear increase in the island size, with higher slopes for lower rates. At higher rates, there is a very small increase in the average island size, since new small islands are continuously forming. On the other hand, at low rates few islands are formed, whose size increases considerably as new atoms join them. After that there is an abrupt change in size, which corresponds to coalescence, and the last stage observed corresponds to size increase of a single island.

Percolation properties of the system were also analyzed. The percolation coverage degree $`\theta _p`$ is defined as the coverage degree at which at least one island is formed that crosses over the system, at least in one direction. Figure 10 shows $`\theta _p`$ as a function of $`log(R)`$, with $`R=D/k`$, where D is the diffusion rate of the atoms on the surface ($`D=436s^1`$ in this case). This was calculated for a $`100\times 100`$ surface and averaged over several simulations. As it can be observed, $`\theta _p`$ increases as $`R`$ does, that is, as atom adsorption rate decreases. This is due to the fact that at lower rates islands are more compact, therefore they take more time to join one another. On the other hand, at high adsorption rates islands branch and they are more likely to join together to form a single island at lower $`\theta `$’s.

#### 4.2.3 Model B: Preferential Adsorption on Terraces. Structural Aspects.

As it was previously detailed, model B considers that adatom deposition occurs only at sites with coordination 0. In electrochemical terms, this means that the reduction reaction occurs only on defect-free terraces. In this case, we have shown above that the evolution of the coverage degree (and of the current) is non-trivial and is given by a ”stretched exponential” function which emerges from a scenario similar to that proposed by Palmer et. al. . The exponent $`\alpha `$ depends on the applied overpotential, and in all cases is less than one. Simulations performed under conditions analogous to those of Figure 7 yielded images that are qualitatively very similar to those of model A presented there. However, it was found that for model B the time at which each degree of coverage occurs is longer, that is to say, the filling in of the surface is slower. Comparative figures are given in Table 4. It can be noticed that this effect is stronger at higher coverages and at high deposition rates. As explained before this is due to the lesser availability of adsorption sites.

In Figure 11 we can observe the evolution of the total density of islands $`N`$ and the density of monomers $`N_1`$ as a function of $`\theta `$. A comparison between this plot and the one of the previous model in Figure 8 shows that at high rates the island and monomer density are larger, while at low rates the curves practically agree. This is due to the fact that at high adsorption rates atoms have little time to diffuse and reorder, and in turn those that adsorb do it in sites not surrounded by atoms. All this causes the formation of a larger amount of islands, whereas at low rates, in spite of adsorption only at neighbor- free sites atoms have sufficient time to diffuse and join the existing islands and the net growth occurs in a way similar to that of the previous model. In addition to the larger island density we can observe a shift at high coverage (between 0.4 and 0.6), which indicates that there are more islands in that zone, that is to say, the islands take longer to coalesce and do it at higher $`\theta `$s.

Figure 12 shows the average size of islands as a function of the coverage degree. If we compare this plot with the one of Model A we can observe that the abrupt change in size shifts to higher $`\theta `$s, especially for the case of an atom entrance rate of $`10^3s^1`$, which is the slowest rate. This is due to island compactness, which causes islands to take longer in joining and forming a single one.

#### 4.2.4 Model C: Preferential Adsorption on the Edges. Structural Aspects.

This model considers that atoms on island edges adsorb at a rate higher than that of atoms which adsorb on the terrace. Adsorption rates of $`10^2s^1`$, $`10^3s^1`$ and $`10^4s^1`$ on the terrace ($`k_1`$), and up to $`10^{+2}s^1`$ on the edges($`k_2`$) were studied. For very large rate differences a single island is formed in our simulation system, which grows from the edges. As adsorption rates become similar the formation of more islands is observed. Figure 13 shows the state of the surface at $`\theta =0.1`$ for different values of $`k_1`$ and $`k_2`$, entrance rate on the terrace and on the edges respectively. The left-hand column corresponds to $`k_2=10^{+2}s^1`$. Under these conditions, only for $`k_1=10^2s^1`$ two islands were formed, while at lower values of $`k_1`$ a single island appears. On the other hand, the results on the right column show that for similar values of $`k_1`$ and $`k_2`$ several islands are formed. In all cases, however, the island shape is quite compact and dendritic shapes are not observed, as they were found at high rates in the other models.

## 5 Conclusions

In this paper we considered several aspects of 2-D metal deposition by means of dynamic Monte Carlo simulations As a model system we considered $`Ag`$ deposition on Au(100) under various deposition conditions using interatomic potentials adequate for metals and a realistic diffusion simulation model.

The reduction of $`Ag^+`$ ions was modeled in three different ways: in the first and simplest case, deposition of atoms was assumed to occur on any unoccupied site of the lattice independently of the site surroundings. Here the evolution of the coverage as a function of time was given in a trivial way by an exponential function, while the structure of the surface depended on the deposition rate. In the second case, particle adsorption was allowed only on defect-free terraces, that is, only on those sites characterized by adsorbate-free surroundings. It was observed that the evolution of the coverage rate is given by a function $`\theta =1exp(\beta t^\alpha )`$, where the exponent $`\alpha `$ depends on the deposition rate, that is, on the overpotential applied, though in all the cases considered it turned out to be smaller than one. At high deposition rates a change in this exponent was also observed, which indicates different types of evolution of the surface structure. Thus a relationship exists between the coverage evolution - and therefore the current evolution - and the change in the surface structure. In the third deposition model, the ion reduction was considered to be possible on edge sites as well as on sites located on clean terraces. A higher deposition rate was always used on the edges. In this case, the coverage degree seems to follow again a $`\theta =1exp(\beta t^\alpha )`$ law, but with exponent $`\alpha `$ always larger than one. This kind of behavior can be compared with the predictions of the Avrami model , where the exponent $`\alpha `$ is two or three, depending on whether the nucleation is progressive or instantaneous. As regards the evolution of surface morphology, there was a remarkable reduction in the number of islands, and a single island was formed in our simulation system when there was a very large difference between the rates over the two types of sites. The exponent values reported must be considered to give a qualitative tendency, in the sense that we did not study systematically the finite-size effects and they correspond to a single simulation run. These features may be important in determining the exact values of the exponents.

## 6 Acknowledgements

We thank CONICET, SeCyt-UNC, Agencia Córdoba Ciencia, Program BID 1201/OC-AR PICT $`N^o`$ 06-04505 and Fundación Antorchas (M.G.D.) for financial support. We also thank to E. Albano and R. Moneti for Hoshen-Kopelman subroutine. Part of the present calculations were performed on a Digital workstation donated by the Alexander von Humboldt Stiftung, Germany. Language assistance by Karina Plasencia is gratefully acknowledged.

## 7 Tables

Table 1

Fitted $`\alpha `$ and $`\beta `$ values from simulations with various adsorption rates $`k`$ onto an unoccupied site . These parameters were obtained by fitting of the law $`\theta =1exp(\beta t^\alpha )`$ from simulations with model B described in the text.

| $`k(rate)`$ | $`\alpha `$ | $`\beta `$ |
| --- | --- | --- |
| $`10^2`$ | 0.722 | 12.72 |
| $`10^1`$ | 0.748 | 2.914 |
| $`10^0`$ | 0.799 | 0.588 |
| $`10^1`$ | 0.868 | 0.095 |
| $`10^2`$ | 0.909 | 0.012 |
| $`10^3`$ | 0.951 | 0.001 |

Table 2

Parameters employed in the simulation studies with model C. $`k_1`$ and $`k_2`$ denote the adsorption rates on sites without adsorbed neighbors( ’terrace site ’) and with adsorbed neighbors( ’border site ’) respectively. Crosses indicate the sets of values of rate constants considered .

| $`k1(terrace)\backslash k2(border)`$ | $`10^2`$ | $`10^1`$ | $`10^0`$ | $`10^1`$ | $`2\times 10^2`$ | $`10^2`$ | $`2\times 10^3`$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $`10^2`$ | x | x | x | x | x |  |  |
| $`10^3`$ | x | x | x | x |  | x | x |
| $`10^4`$ | x | x | x | x |  | x |  |

Table 3

Fitted $`\alpha `$ and $`\beta `$ values for various simulation conditions with model C. They are obtained by fitting of the law $`\theta =1exp(\beta t^\alpha )`$ .

| $`k1(terrace)`$ | $`k2(border)`$ | $`\alpha `$ | $`\beta `$ |
| --- | --- | --- | --- |
| $`10^2`$ | $`10^2`$ | 2.3 | 63 |
| $`10^2`$ | $`10^1`$ | 2.4 | 1.3 |
| $`10^2`$ | $`10^0`$ | 1.9 | 0.03 |
| $`10^2`$ | $`10^1`$ | 1.5 | 0.005 |
| $`10^2`$ | $`2\times 10^2`$ | 1.2 | 0.006 |
| $`10^3`$ | $`10^2`$ | 2.1 | 36. |
| $`10^3`$ | $`10^1`$ | 2.0 | 0.51 |
| $`10^3`$ | $`10^0`$ | 2.0 | 0.009 |
| $`10^3`$ | $`10^1`$ | 1.7 | $`7.1\times 10^4`$ |
| $`10^3`$ | $`10^2`$ | 1.3 | $`4.3\times 10^4`$ |
| $`10^3`$ | $`2\times 10^3`$ | 1.1 | $`7.4\times 10^4`$ |
| $`10^4`$ | $`10^2`$ | 2.0 | 27. |
| $`10^4`$ | $`10^1`$ | 1.9 | 0.24 |
| $`10^4`$ | $`10^0`$ | 2.0 | 0.002 |
| $`10^4`$ | $`10^1`$ | 2.1 | $`2.9\times 10^5`$ |
| $`10^4`$ | $`10^2`$ | 1.7 | $`4.9\times 10^6`$ |

Table 4 Time $`t_\theta `$ at which a given degree of coverage $`\theta `$ occurs in models A and B described in the text. $`k`$ denotes the rate assigned to the process of atom entrance onto an unoccupied site.

| $`\theta `$ | $`\begin{array}{c}t_\theta (k=10^2s^1)\\ \mathrm{m}od\text{el A}\end{array}`$ | $`\begin{array}{c}t_\theta (k=10^2s^1)\\ \mathrm{m}od\text{el B}\end{array}`$ | $`\begin{array}{c}t_\theta (k=10^2s^1)\\ \mathrm{m}od\text{el A}\end{array}`$ | $`\begin{array}{c}t_\theta (k=10^2s^1)\\ \mathrm{m}od\text{el B}\end{array}`$ |
| --- | --- | --- | --- | --- |
| 0.1 | 0.00094 | 0.0013 | 11.3 | 11.6 |
| 0.5 | 0.007 | 0.020 | 65.5 | 88.2 |
| 0.7 | 0.012 | 0.083 | 114. | 177. |
| 0.9 | 0.024 | 194. | 221. | 430. |

## 8 Figure Captions

Figure 1: Scheme of the mechanisms for atom attachment to the growing phase 2-D studied in the present work:

* Particle adsorption occurs on all unoccupied sites on the surface at the same rate $`k.`$
* Particle adsorption occurs only at sites where no adsorbed atoms are around ( ’terrace sites ’) at the rate $`k`$. The adsorption rate on sites that have one or more neighboring adatoms is $`0`$.
* Particles are considered to be discharged preferably at sites with neighboring adatoms ( ’edges ’). A rate $`k_1`$ is assigned to the entrance onto terrace sites and a rate $`k_2`$ ( $`k_2`$ $`>`$ $`k_1`$) is given to atoms entrance into sites that have at least one neighboring adatom.

Figure 2: Evolution of the coverage degree $`\theta `$ as a function of $`kt`$ corresponding to the adsorption model B described in the text (’terrace ’ adsorption). The adsorption rate $`k`$ is given in the figure in units of $`s^1`$.

Figure 3: Plot of $`ln(ln(1\theta ))`$ as a function of $`ln(t)`$ at various deposition rates $`k`$ for Model B described in the text (’terrace ’ adsorption). The adsorption rate $`k`$ is given in the figure in units of $`s^1`$. Slopes (exponents) associated to the linear relation found during the first stage of deposition are detailed in table 1.

Figure 4: Evolution of the coverage degree $`\theta `$ and of the current for models A and B for a deposition rate of $`1s^1`$ .

Black solid line: evolution of $`\theta `$ in model A. Gray dashed line: evolution of $`\theta `$ in model B.

Circles: current evolution in model A. The gray solid line represents the theoretical fitting according to the equation $`d\theta /dt=\beta \alpha t^{\alpha 1}exp(\beta t^\alpha )`$

Diamonds: current evolution in model B. The black dashed line represents the theoretical fitting according to the equation $`d\theta /dt=\beta \alpha t^{\alpha 1}exp(\beta t^\alpha )`$

The exponents $`\alpha `$ and $`\beta `$ are given in Table 1.

Figure 5: Evolution of $`\theta `$ as a function of time for the deposition model C. In these simulations, a deposition rate $`k_2`$ = $`1s^1`$ is considered on the edges and different deposition rates $`k_1`$ are assumed at the terraces. The $`k_1`$ values are given in the figure in units of $`s^1.`$

Figure 6: Current ($`d\theta /dt`$) as a function of time for deposition model C. The full line are simulation results, and the dashed lines show the curves predicted by the function $`d\theta /dt=\beta \alpha t^{\alpha 1}exp(\beta t^\alpha ).`$ Each plot employs the parameters obtained through minimum square fits of the simulation data. The adsorption rates employed for the simulations were the following: Top box: $`k_1=10^2s^1`$ ; $`k_2=2\times 10^2s^1`$; middle box: $`k_1=10^3s^1`$ ; $`k_2=1s^1`$; box at the bottom: $`k_1=10^4s^1`$ ; $`k_2=10^2s^1`$.

Figure 7: Frames showing surface morphology in model A at three different values of $`\theta `$. The deposition rates were $`k=`$ $`10^2s^1`$ for the left-hand column and $`k=10^2s^1`$ for the right-hand column.The coverage degrees are, from top to bottom, $`\theta =0.1`$, $`\theta =0.3`$ and $`\theta =0.6`$.

The associated times are $`t=9.4\times 10^4s`$, $`t=3.5\times 10^3s`$ and $`t=9.2\times 10^3s`$ from top to bottom for the left column and $`t=11.31s`$, $`t=34.22s`$ and $`t=87.26s`$ for the right one.

Figure 8: Plot of island density (top) and monomer density (bottom) as a function of the coverage degree for different adsorption rates according to model A. Rates used in each case are indicated in the body of the figure in units of $`s^1.`$

Figure 9: Average island size $`<s>`$ (defined as the average number of atoms per island divided by the total number of sites $`N_T`$ ) as a function of the coverage degree $`\theta `$, for different adsorption rates according to model A. The rates used in the simulations are included in the figure in units of $`s^1`$. $`N_T=2500.`$

Figure 10: Coverage degree at which percolation occurs, $`\theta _p`$ as a function of $`log(R)`$, where $`R=\frac{D}{k}`$, $`D`$ being the diffusion rate of an adatom on the clean surface and $`k`$ the particle entrance rate. Simulations were performed on a $`100\times 100`$ square lattice.

Figure 11: Plots of island density (top) and monomer density (bottom) as a function of coverage degree for different adsorption rates according to model B. Rates used in each case are indicated in the body of the figure in units of $`s^1.`$

Figure 12: Average island size $`<s>`$ as a function of the coverage degree $`\theta `$, for different adsorption rates according to model B. The rates used in the simulations are included in the figure in units of $`s^1.`$

Figure 13: Frames showing the surface structure in simulations with model C at a coverage degree $`\theta =0.1`$ for different adsorption rates on terraces and edges. Frame 1: $`k_1=10^2s^1`$; $`k_2=10^2s^1`$; Frame 2: $`k_1=10^2s^1`$; $`k_2=2\times 10^2s^1`$; Frame 3: $`k_1=10^3s^1`$ ; $`k_2=10^2s^1`$; Frame 4: $`k_1=10^3s^1`$; $`k_2=10^2s^1`$; Frame 5: $`k_1=10^4s^1`$; $`k_2=10^2s^1`$; Frame 6: $`k_1=10^4s^1`$ ; $`k_2=10^2s^1`$.