Microscopic processes accompanying 2D condensation in metallic submonolayers

surface science ELSEVIER

Surface Science 355 (1996) 248-254

Microscopic processes accompanying 2D condensation in metallic submonolayers J.M. R o g o w s k a * Institute of Experimental Physics, University of Wroclaw, Pl. Maxa Borna 9, 50-204 Wroctaw, Poland Received 26 September 1995; accepted for publication 6 December 1995

Abstract A continuous space model of a submonolayer with non-additive lateral interactions is applied to a Monte Carlo simulation of two-dimensional (2D) condensation phenomena. The system of Cu/W(110) is chosen to represent metallic adsorbates with attractive interactions. Equilibrium properties of adlayer as a function of temperature at a constant concentration of adatoms are analysed. Morphology and energetics of Cu submonolayers are studied by analysing island sizes, number of islands, and number of monomers as well as total and lateral configurational energy of the system. The simulations show that the 2D condensation strongly depends on the concentration of adatoms as a result of changes in the degree of association in a 2D imperfect gas. The obtained heat of 2D condensation is in agreement with experiment and has its maximal value at relatively low coverage. Several stages of condensation are distinguished and their evolution with coverage and temperature is discussed.

Keywords: Adatoms; Chemisorption; Computer simulations; Construction and use of effective interatomic interactions; Copper; Equilibrium thermodynamics and statistical mechanics; Semi-empirical models and model calculations; Surface thermodynamics; Tungsten

I. Introduction

The aim of this paper is to contribute to a better understanding of the phenomenon of submonolayer condensation of noble metals on W(ll0) reported in papers of Kotaczkiewicz and Bauer [ 1]. The phase boundaries between single-phase (2D gas) and two-phase (2D gas + 2D condensate) regions are determinated in Ref. [ 1] from discontinuities of the work function slope with temperature caused by a rapid change in local adatom coordination. These phase boundaries are calculated via the 2D van der Waals equation of state assuming that condensing particles are mobile, but very little * E-mail: [email protected]. 0039-6028/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved PII S0039-6028 (95) 01373-3

is known about the morphology of the submonolayer during the process of condensation. Recently, there is much interest in the microstructure of metallic submonolayers [2] but most of the theoretical work is focused on modelling non-equilibrium processes. Equilibrium processes like 2D condensation were studied in some earlier papers [3], however, the commonly used lattice gas model is too simplified for the relatively smooth W(110) surface. In this paper, the more realistic continuous space model, studied in Ref. [4], is extended to include non-additivity of lateral interactions and to determine the island size distribution function, which is widely used in investigations of submonolayer morphology [2]. Using a standard canonical Monte

J.M. Rogowska/Surface Science 355 (1996)248-254

Carlo (MC) procedure [5] the simulations are performed for the Cu/W(ll0) system for which the complete set of parameters of the intermolecular potential is available from Gollisch [6]. To the best author's knowledge, this system has not been modelled recently by a more accurate method and this is the only set of parameters of the effective potential available for the W(110) substrate.

2. The model In the present work the substrate is assumed to be unreconstructed and acting as an external potential W(x,y) on the adsorbate. W(x,y) is defined for the adatom-surface coordinate z where the adatom-surface potential is minimal at given x,y coordinates. The values of this potential in the Cu/W(ll0) system are calculated from analytic expressions given in Ref. [6]. The lateral interaction in the submonolayer of N adatoms is described in terms of the generalised Morse potential [7]

Vi(xi,yl)=b

e-2Z'iJ+ a j~i

e -z'

,

(1)

\j~,i

where V~ is the potential energy of the interaction of the ith adatom at r~(xi,y,) with ( N - 1 ) adatoms located at a distance rij=lri-rj[. In Eq. (1) the b, a, 2 and # are adjustable parameters. If # ~ 1 the pair contributions to V~ are influenced by the entire system under study. The total configurational energy U r of N adatoms is taken as

UT(N)=

N

W(x,,y,)

i=1 N

+ ~ Vi(xi,y~)=US(N)+UL(n),

(2)

i=1

where the interaction energy of the adatoms with substrate Us and the lateral interaction energy U L can be discussed separately. There should be given some comments concerning the non-additivity of the potential Vi(N ). Introducing the parameter of non-additivity P as the ratio of the strongest (one first neighbour) to the weakest (Z first neighbours) bond in the 2D system of adatoms, we have P = V(1)/V(Z), where Z is the number of nearest-neighbour sites and

249

V(N) is the energy per bond when an adatom has N first neighbours actually present. For the W(110) surface, which may be considered as a somewhat distorted hexagonal surface, Z = 6 is used. Applying Vi (Eq. (1)) with the parameters given by Gollisch [8] to calculate the value of V(N) for Cu on W(ll0) we obtain P = 1.28, while in the recent calculations of Breeman [9] for the Cu/Cu(111) system, where also Z = 6, the value of P is equal to 1.36. The above rough estimation enables us to believe that the non-additivity in the Gollisch potential is reasonable. Moreover, the non-additivity in this system contributes less than several per cents to the total energy, therefore the error introduced by the pairwise summation in Eqs. (1) and (2) does not change significantly the mechanisms underlying the phase transition we have been focusing on. 3. The Monte Carlo simulation procedure The canonical Monte Carlo simulation in continuous 2D space is performed using the conventional Metropolis algorithm [10] along with the Verlet tables [ 11]. Migration of adatoms is controlled by the total configurational energy difference AU T gained by the system on displacing adatoms. Adatom displacement is generated randomly with a maximal trial move allowing the adatom to jump from one substrate unit cell onto a neighbouring one. The changes in the configuration of adatoms are accepted if they lower the energy of the system, or rejected with a probability exp(-AUT/kT) if this energy is increased. T is the temperature of the system and k is the Boltzmann constant. The total configurational energy of the Cu/W(ll0) system is calculated from Eq. (2), where the effective pair potential Vi is cut off at a value of three equilibrium distances. The simulation is carried out typically with 80-100 adatoms confined in a rectangular box by periodic boundary conditions. The "computer experiment" is performed at a constant concentration of adatoms by a simulated cooling procedure to follow the experiment of Kotaczkiewicz and Bauer [ 1]. The temperature is reduced from 1900 K down to 100 K with steps

250

J. ~ Rogowska/Surface Science 355 (1996) 248-254

of AT= 100 K. In this simulation procedure the memory of the previous configuration is not lost, i.e. the equilibrium configuration at a temperature T is taken as the initial configuration at T-T - AT. All thermodynamic data are averaged over 12 x 103 MC steps after 4 x 103 equilibrating MC steps at each temperature. A record of the evolution of the averages is analyzed and it is assumed that the condensate is in equilibrium with its own vapour at every temperature under study. It means that the system at a given 0 is approximated by a sequence of equilibrium states as temperature decreases. The possibility of existence of metastable states in the twophase coexistence region is discussed in Section 4. On the other hand, the results of MC simulations obtained at low temperatures where the equilibration occurs slowly (in this system at T< 500 K) are always uncertain. However, the phase transition discussed in this work takes place at higher temperature and some uncertainty of low temperature results is irrelevant. The quantities sampled during the course of MC simulation are: (UT>, (U s> and (uL>, (compare Eq. (2)), the island size distribution function (n(s)>, the average size g and the areal density n of islands composed of more than one adatom, as well as the areal density of monomers nl. The angular brackets represent an average over the MC trajectory. The island size distribution (n(s)> denotes the average number of islands (per site) containing s adatoms at a given temperature and coverage. The coverage 0 is equal to N/Ns, where Ns is the number of adsorption sites available in the monolayer. The remaining quantities are defined as follows: g.~_ ENS = 2 s(n(s)>

<.(s)>

n 1 = (n(s

=

N

-3.4 1 ~o~ E -3.5

v -3.81

1)>,

'

n = ~ = 2 .

barrier between sites are equal to 3.50 and 0.30 eV, respectively. A more detailed description of adsorption sites on W(ll0) can be found in Refs. [6,4]. The total configurational potential energy (U T> sampled in the course of the MC simulation is shown in Fig. 1 as a function of temperature for different coverages. In this work, the study of the condensation phenomenon is restricted to intermediate low coverages. Two additional curves, for 0 = 0.0001 and 0 = 0.87, are shown in Fig. 1 in order to illustrate the asymptotic behaviour of the system. At extremely low 0 very few adatoms have a chance to meet one another so adatoms remain single. A more precise localisation of adatoms at the potential minima of the adsorption sites leads to a linear decrease of with lowering T. At coverages close to a monolayer a high degree of association in a very dense 2D gas does not change significantly during condensation, therefore, (UT(T)> returns back to the linear dependence on T. In Fig. 1 we can distinguish high temperature

(3)

~'~

!

-3.9I

4. Results and discussion

-4.0 5

The submonolayer of Cu on W(ll0) is pseudomorphical with a 13% lattice mismatch. Binding of a single Cu occurs in a hollow site with shallow dual minima. The binding energy and the energy

t t 10 15

T ( 102 K} Fig. 1. Total configurationa] energy (UT> versus temperature for the Cu/W(llO) system. Subrnonolayer coverage 0 is thc parameter of the curves.

J.M. Rogowska/Surface Science 355 (1996) 248-254

regions associated with a 2D gas phase, low temperature regions of a 2D condensed phase and transition temperature regions where the 2D condensate and the 2D gas phase coexist. In the transition region the (UT(T)) curves have points of inflection that can be regarded as phase boundary points. The results presented above are insufficient for the construction of a detailed phase diagram. However, the location of those points versus Tand 0 is in qualitative agreement with the coexistence curve given in Ref. [ 1]. The value of the abrupt change of the total energy IA(UT)I near the points of inflection may be taken as a measure of heat of 2D condensation. As is apparent from Fig. 1, the IA(UT)I depends on 0 and at 0,~0.05 has its maximal value equal to 0.32 eV, (see also Fig. 2). The experimentally obtained [12] heat of 2D condensation, at Tgl000 K and 0=0.1, was found to be equal to 0.30 eV. On the other hand, we have [A(UT)0=o.x[=0.31 eV. This agreement confirms the validity of the intermolecular potential and the model used in this work. At this point, some remarks concerning the smooth behaviour of (UT(T)) in the transition region have to be made. The phenomenon of condensation in an infinite 3D system is an example of first-order transition, which is characterised by a kink singularity in the internal energy at Tt,.

"E -

0.1

t _~ -03 ",t-0.4[ I

I

I

5

10

15

T ( lOZK}

Fig. 2. Lateral energy (U L) versus temperature for the Cu/W(110) system. Submonolayercoverage 0 is the parameter of the curves.

251

As it is seen in Fig. 1, in the considered finite 2D system, this singularity is rounded off over a wide temperature region ATtr(0). In MC simulations of finite systems the finite size effect [5] is responsible of rounding over 6 T~ 2k T2tr/[iA< U T)l/if], where L is the linear dimension of the simulation cell in units of lattice spacing and d is dimensionality of the system. This finite size rounding in the considered system gives f T ~ 1 0 K, taking /~=N. Moreover, in MC simulations of a first-order transition, it is difficult to avoid metastable states [5]. In order to estimate f'Tcaused by metastability, several additional runs in which the system was heated were performed. It was found, that (UT(T)) exhibits hysteresis with f'T,,'~100 K. However, the sum of fiT and f'T is still smaller than ATtr. Another source of broadening of ATtr may be a mobile-localised transition, that is hindered in this temperature region. At low temperature adatoms are localised with respect to the surface, whereas at higher temperature they are mobile. This problem needs further investigation. The total energy is the sum of ( U s) and (uL). We will discuss ( U L) separately (see Fig. 2) because the changes in the degree of association in the adsorbate are directly reflected by changes of (UL(T)). At the lowest coverage shown in Fig. 2, there is a perfect 2D gas of non-interacting monomers ( ( u L ) ~ 0 ) in the high temperature region, and a 2D condensate phase made of dimers ( ( U L ) ~ 0 . 1 1 e V = l / 2 of the binding energy of dimer) in the low temperature region. At higher coverages the 2D gas becomes more and more imperfect, therefore at the critical coverage 0o equal to 0.215 according to Ref. [ 1] (or 0c = 0.5 according to the Ising model) we cannot distinguish the 2D gas from the 2D condensate. Analysing our results not presented in this paper we can state that the critical coverage is rather closer to 0.3 than to 0.5. We consider the 2D gas to be perfect if adatoms are mobile and they do not interact, i.e. the value of ( U L) is negligible. For example, this condition is fulfilled for 0<0.001 and T>700 K. On the contrary, the imperfect gas consists of aggregates (unstable islands) of different sizes (compare Fig. 4) which are instantaneously created and annihilated by adatomic thermal motions. The decrease of high temperature values of ( u L ) , (see Fig. 2),

J.M. Rogowska/SurfaceScience 355 (1996) 248-254

252

results from the increase of the average size of the aggregates with increasing 0. It can be also noticed that the low temperature limit of the 2D gas phase shifts towards higher temperatures with increasing 0, i.e. from 700 K at 0=0.001 to 1400 K at 0=0.2. The changes of ( U L) with T are governed by the decrease of the average number of monomers n, Ns. The striking similarity in the shape of nlNs(T) and (UL(T)) curves, presented in Fig. 3a for 0 = 0.012, is typical of all considered coverages. From Fig. 3a it is also seen that the 2D gas consists mainly of monomers while in the 2D condensate region the number of monomers is close to zero. Analysing the average island size g and the average number of islands nN~ shown in Fig. 3b, together with ( U L) and the average number of monomers n~N~ shown in Fig. 3a, the following characteristic stages in the sequence of equilibrium states as the temperature decreases, at constant 0, can be distinguished: Reversible nucleation of mobile monomers into

o)

80

-0.1 E o

60

-o -0.2 O

40

QJ

-.. t-

-0.3

small aggregates: ( u L ) , nl, n and g are nearly constant. Irreversible nucleation of mobile monomers into aggregates: ( U L) and nl slowly decrease, n increases, g is nearly constant. Aggregation of mobile monomers into existing islands: ( U L) and n, rapidly decrease, n reaches its maximal value, g slowly increases. Coalescence of existing islands into larger ones: ( u L ) , nl and n decrease, g rapidly increases. Internal ordering within large islands: ( U L) slowly decreases, n, ~0, n and g are constant. All the stages appear in the succession presented above only for relatively low coverages. Consider the situation shown in Fig. 4, where the nNs(T) curves shift to higher temperatures with increasing 0. At 0 >0.15, there is no nucleation stage in the considered temperature range, because even at a temperature comparable to the temperature of desorption (~ 1400 K) a considerable number of aggregates already exists in the 2D imperfect gas phase. On the other hand, for dilute adsorbates (0<0.01) the condensation never approaches the coalescence stage because of the large distances between islands. During the coalescence stage islands grow not by strict coalescence (direct island-island interaction) but mainly by ripening (island-island equilibration via vapour). Inspection of the slope of nlNs in Fig. 3a reveals that in this stage the decreasing number of mobile monomers is still high.

20 v

15 0.012

c

O/

, 5

, 10

I 15

T (10 2 K )

Fig. 3. System of N = 7 8 adatoms on the surface with Ns =6400 lattice sites (0=0.012). (a) Lateral energy ( U L) and average number of monomers nlNs versus temperature, (b) Average number of islands nN~ and average island size g versus temperature.

0j~02S

/

/

~

i T (102K1

Fig. 4. Average number of islands nN, versus temperature.

Submonolayer coverage is the parameter of the curves.

J.M. Rogowska/Surface Science 355 (1996) 248-254

As long as 0 is not exceeding 0.06 the system does not reach an absolute equilibrium, which corresponds to a single optimally coordinated island, but it ends up in a stationary state where a population of small islands is in local equilibrium with monomers. When 0>0.06, the adsorbate collapses in a single 2D aggregate, which with lowering T evolves into a compact island with a highly symmetric shape (compare the values of 650 K. This is in agreement with experiment, where the 2D roughening temperature in Cu islands on W ( l l 0 ) was reported [12] to be equal to 620 K.

5. Conclusions The obtained constant coverage MC simulation of a metallic submonolayer with non-additive lateral interactions shows that the 2D condensation is a first-order transition I-5] with a wide transition region and the latent heat depending on 0. The obtained heat of 2D condensation is in agreement with experiment and has its maximal value at

253

X (aw)

Fig. 5. Fragment of a snapshot of the adatomic configuration for 0=0.05 and T = 700 K. Lattice grid in the Y direction is compressed by a factor of x/2 and aw denotes the lattice constant of the W(110) substrate.

0~0.05. In the transition region the considered sequence of equilibrium states displays the characteristic stages of island growth known from nonequilibrium processes [2], namely nucleation, aggregation and coalescence. This succession is perturbed by the lack of a coalescence stage for 0<0.01 and the vanishing nucleation stage for 0>0.15.

Acknowledgements The author acknowledges stimulating discussion with J. Kotaczkiewicz and thanks M. Maciejewski from the Institute of Informatics, University of Wrodaw for his help in the computations. This work was supported by the Polish State Committee of Scientific Research with Grant 2 P 30211305.

References [1] J. Kotaczkiewicz and E. Bauer, Surf. Sci. 151(1985) 333; Phys. Rev. Lett. 53 (1984) 485. 12] G.S. Bales and D.C. Chrzan, Phys. Rev. B 50 (1994) 6057; C. Ratsch, A. Zangwill, P. ~milauer and D.D. Vvedensky, Phys. Rev. Lett. 72 (1994) 3194; P. Blandin, C. Massobrio and P. Ballone, Phys. Rev. B 49 (1994) 16 637; C. Ratsch, P. ~milauer, A. Zangwill and D.D. Vvedensky, Surf. Sci. 329 (1995) L599. [3] A. Milchev and K. Binder, Surf. Sci. 164 (1985) 1; L.D. Roelofs and R.J. Bellon, Surf. Sci. 223 (1989) 585. [4] N. Georgiev, A. Milchev, M. Paunov and B. Dtinweg, Surf. Sci. 264 (1992) 455; A. Milchev and E. Bauer, Surf. Sci. 222 (1989) 163.

254

J.M. Rogowska/Surface Science 355 (1996) 248-254

[5] K. Binder and D.W. Heermann, in: Monte Carlo Simulation in Statistical Physics: An Introduction, Ed. P. Fulde (Springer, Berlin, 1988). [6] H. Gollisch, Surf. Sci. 175 (1986) 249. [7] H. Gollisch, Surf. Sci. 166 (1986) 87. [8] The parameters for Vi (Eq. ( 1)) are extracted from Gollisch results (Surf. Sci. 175 (1986) 249) by T. Regenstein in MSc Thesis, University of Wrodaw, 1994. [9] M. Breeman, G.T. Barkema and D.O. Boerma, Surf. Sci. 323 (1995) 71.

[10] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth and A.H. Teller, J. Chem. Phys. 21 (1953) 1087. [-11] L. Verlet, Phys. Rev. 159 (1967) 98. [,12] J. Koiaczkiewicz and E. Bauer, Surf. Sci. 155 (1985) 700. [13] A. Zangwill, in: Physics at Surfaces (Cambridge Univ. Press, Cambridge, 1988) p. 263. [14] H. Shao and S. Liu, H. Metiu, Phys. Rev. B 51 (1995) 7827. [15] P. Blandin, C. Massobrio and P. Ballone, Phys. Rev. Lett. 72 (1994) 3072.