Rigid Molecular Model for the Assembly Characteristics and Optimal

1474

Langmuir 2003, 19, 1474-1485

Rigid Molecular Model for the Assembly Characteristics and Optimal Structure in Molecular Monolayers of Alkanethiols on Au(111)† Niels Grønbech-Jensen,*,‡,§ Atul N. Parikh,‡ Keith M. Beardmore,| and Rashmi C. Desai⊥ Department of Applied Science, University of California, Davis, California 95616, Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, Computational Nanoscience Group, Motorola Incorporated, Los Alamos Research Park, Los Alamos, New Mexico 87545, and Department of Physics, University of Toronto, Toronto, Ontario M5S 1A7, Canada Received August 2, 2002. In Final Form: December 9, 2002 We present a simple molecular model for the self-assembly of alkanethiols on gold. The model, a rigid rod representation of a molecule which is constrained to a given distance from the gold surface, allows direct long-time simulations of large-scale molecular ensembles (104-105 molecules) on desktop workstations. As a result, the model allows for efficient studies of evolution and ordering of, for example, orientational order and domain patterns in a full range of monolayer coverages. The model is parametrized entirely by existing literature on atomic and molecular scale interactions. Extensive simulations of molecular selfassembled monolayer domain formation at optimal (packing commensurate with a gold surface (111) structure) and suboptimal packing are conducted and presented. The results show close correspondence between the model features and the existing, and emerging, picture observed through experimental characterization of self-assembled monolayers on Au(111). This strong experimental validation of the conformationally insensitive molecular model suggests that the conformational degrees of freedom are not essential for the self-assembly of alkanethiols on gold. It appears that the interplay between the substrateheadgroup and chain-chain interactions determines the self-assembly process and the emergent molecular structures. The presentation of simulation results for different molecular surface coverages is used to derive a primitive two-dimensional isothermal phase diagram. The latter was found to be in good general agreement with available experimental data and provides insight into the formation and growth mechanism of monolayers. The work suggests directions for a minimal approach to studying at least some of the complexity contained in molecular self-assembly processes.

I. Introduction Interest in surface self-assembly1,2 persists for many reasons. First, it offers a simple approach to constructing molecularly defined interfaces with desired surface characteristics (e.g., surface free energy and chemical composition) relevant to a host of applications ranging from biological sensing3 to micro/nanofabrication.4 Second, owing to their high structural definition and molecularlevel control in surface chemical composition, surface selfassembled monolayers (SAMs) provide ideal candidates for exploring the microscopic basis of a variety of interfacial processes, for example, wetting, adhesion, and interfacial * Corresponding author. † Part of the Langmuir special issue entitled The Biomolecular Interface. ‡ University of California. § Lawrence Berkeley National Laboratory. | Motorola, Inc. ⊥ University of Toronto. (1) Bigelow, W. C.; Pickett, D. L.; Zisman, W. A. J. Colloid Sci. 1946, 1, 513. (2) Ulman, A. Introduction to Ultrathin Organic Films: From Langmuir-Blodgett to Self-Assembly; Academic: San Diego, 1991. (3) Rubenstein; Steinberg, S.; Tor, Y.; Shanzer, A.; Sagiv, J. Nature 1988, 332, 426. Haussling, L.; Ringsdorff, H.; Schmitt, F.-J.; Knoll, W. Langmuir 1991, 7, 1837. Prime, K. L.; Whitesides, G. M. Science 1991, 252, 1164. (4) See, for example: Kumar, A.; Abbott, N. L.; Kim, E.; Bieuyck, H.; Whitesides, G. M. Acc. Chem. Res. 1995, 28, 219. Lercel, M. J.; Craighead, H. G.; Parikh, A. N.; Seshadri, K.; Allara, D. L. Appl. Phys. Lett. 1996, 68, 1504.

charge transfer.5-7 Third, SAMs offer unique new opportunities in addressing fundamental issues in physics of low-dimensional systems owing to their pseudo-2D character and complex phase behavior.8 In this regard, a specific amphiphile/substrate system that has gained considerable attention over the past decade is based on organosulfur amphiphiles and gold surfaces.9,10 In a typical example, a simple immersion of a Au substrate into a dilute ethanolic solution of long-chain alkanethiols (CH3(CH2)nSH, with n > 10) results in a highly ordered SAM. The ease in preparation through versatility, high reproducibility, and a consistent degree of structural definition has led to recognition of this system as a model class of SAMs.11-13 (5) Nuzzo, R. G.; Fusco, F. A.; Allara, D. L. J. Am. Chem. Soc. 1987, 109, 2358. (6) Holmes-Farley, S. R.; Reamey, R. H.; McCarthy, T. J.; Deutch, J.; Whitesides, G. M. Langmuir 1985, 1, 725. (7) Dubois, L. H.; Zegarski, B. R.; Nuzzo, R. G. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 4739. (8) Zharnikov, M.; Grunze, M. J. Phys.: Condens. Matter 2001, 13, 11333. (9) Nuzzo, R. G.; Allara, D. L. J. Am. Chem. Soc. 1983, 105, 4481. Bain, C. D.; Whitesides, G. M. Angew. Chem., Int. Ed. Engl. 1989, 101, 506. (10) Dubois, L. H.; Nuzzo, R. G. Annu. Rev. Phys. Chem. 1992, 43, 437. (11) Ulman, A. Chem. Rev. 1996, 96, 1533. (12) Whitesides, G. M.; Laibinis, P. E. Langmuir 1990, 6, 87. (13) Delamarche, E.; Michel, B.; Biebuyck, H. A.; Gerber, C. Adv. Mater. 1996, 8, 719.

10.1021/la026344f CCC: $25.00 © 2003 American Chemical Society Published on Web 02/06/2003

Molecular Model for Monolayers of Alkanethiols

The static structure of alkanethiols on Au(111) has been exhaustively studied using a highly diverse range of experimental and computational characterization capabilities.8,10 Taken together, these studies serve to establish many structural attributes at the molecular and mesoscopic levels. A general picture supported by these studies suggests that the sulfur coordinates strongly to the gold surface, polymethylene chains adopt a predominantly alltrans conformation and tilt uniformly at ∼30-40° from the surface normal, and the methyl tailgroup is exposed at the film-ambient interface. Specifically, it has been estimated using thermal desorption measurements that S coordinates to the Au(111) surface with a chemisorption enthalpy of 126 kJ/mol.7,14 Precise placement of S headgroups relative to the gold lattice has been variously determined using a combination of experimental (grazingincidence X-ray diffraction,15 He diffraction,16-18 highresolution scanning tunneling microscopy (STM) measurements19) and computational (ab initio and density functional formalisms) approaches.20,21 These studies have consistently established that the S headgroup shows a strong preference for 3-fold hollow sites on the triangular Au lattice (spacing, 2.884 Å) at 1/3 occupancy in wellformed SAMs forming a c(4×2) superlattice of a commensurate based hexagonal arrangement. The state of the polymethylene chains of the adsorbates has been derived primarily using infrared and Raman spectroscopic measurements,22,23 but also using corroborating evidence from near-edge X-ray absorption fine structure (NEXAFS) spectroscopy,24 sum frequency generation (SFG) spectroscopy,25 grazing incidence X-ray diffraction (GIXRD),16-18,26 ellipsometry,27 and molecular simulations.28,29 These studies confirm a predominantly all-trans conformation for the polymethylene chain with low gauche defect content at room temperature; the chains are tilted 3040° uniformly in the direction of a next-nearest neighbor (NNN). The chains are further assessed to be organized into a crystalline-like packing composed of two chains per unit cell with their C-C-C backbones orthogonal to each other in a highly efficient interchain packing configuration. In-plane translational order determined using GIXRD and STM measurements has allowed the elucidation of domain structures, revealing the existence of tilt boundaries, S positional order, and extended correlations between the adsorbates in these SAMs. Recently, there has been some effort in reconciling the structural details that have (14) Lavrich, D. J.; Wetterer, S. M.; Bernasek, S. L.; Scoles, G. J. Phys. Chem. B 1998, 102, 3456. (15) Fenter, P.; Eisenberger, P.; Liang, K. S. Phys. Rev. Lett. 1993, 70, 2447. (16) Camillone, N.; Chidsey, C. E. D.; Liu, G. Y.; Putvinski, T. M.; Scoles, G. J. Chem. Phys. 1991, 94, 8493. (17) Camillone, N.; Chidsey, C. E. D.; Liu, G. Y.; Scoles, G. J. Chem. Phys. 1993, 98, 3503. (18) Camillone, N.; Leung, T. Y. B.; Scoles, G. Surf. Sci. 1997, 373, 333. (19) Poirier, G. E. Chem. Rev. 1997, 97, 1117. (20) Sellers, H.; Ulman, A.; Shnidman, Y.; Eilers, J. E. J. Am. Chem. Soc. 1993, 115, 9389. (21) Beardmore, K. M.; Kress, J. D.; Grønbech-Jensen, N.; Bishop, A. R. Chem. Phys. Lett. 1998, 286, 40. (22) Laibinis, P. E.; Whitesides, G. M.; Allara, D. L.; Tao, Y. T.; Parikh, A. N.; Nuzzo, R. G. J. Am. Chem. Soc. 1991, 113, 7152. (23) Bryant, M. A.; Pemberton, J. E. J. Am. Chem. Soc. 1991, 113, 8284. (24) Hahner, G.; Kinzler, M.; Thummler, C.; Woll, C.; Grunze, M. J. Vac. Sci. Technol., A 1992, 10, 2758. (25) Bain, C. D.; Davies, P. B.; Ong, T. H.; Ward, R. N.; Brown, M. A. Langmuir 1991, 7, 1563. (26) Camillone, N.; Chidsey, C. E. D.; Eisenberger, P.; Fenter, P.; Li, J.; Liang, K. S.; Liu, G. Y.; Scoles, G. J. Chem. Phys. 1993, 99, 744. (27) Porter, M. D.; Bright, T. B.; Allara, D. L.; Chidsey, C. E. D. J. Am. Chem. Soc. 1987, 109, 3559. (28) Hautman, J.; Klein, M. L. J. Chem. Phys. 1990, 93, 7483. (29) Mar, W.; Klein, M. L. Langmuir 1994, 10, 188.

Langmuir, Vol. 19, No. 5, 2003 1475

emerged in the studies summarized above in terms of a balance of inter- and intramolecular interactions.8 It is now generally accepted that the dominance of the substrate-headgroup (Au-S) interactions with a compliant, efficient packing of polymethylene chains, determined by the chain-chain interactions, dictates the final, and quite possibly equilibrium, molecular structure in alkanethiolate SAMs. While the structure of alkanethiol monolayers on gold is well characterized, determination of the pathways leading to this structure is only beginning to emerge.30,31 The growth kinetics of alkanethiolate monolayers has been variously studied using ellipsometry, contact angle, a surface-acoustic wave device, second-harmonic generation, near-edge X-ray absorption fine structure, and quartz crystal microgravimetry.31 These studies conclude that the rate of monolayer growth follows Langmuir adsorption kinetics wherein the growth rate is directly proportional to the number of unoccupied sites. Using STM measurements of vapor-phase-generated alkanethiol monolayers at below-saturation coverages, Poirier and Pylant30 have proposed a molecular-level mechanism that parallels that of Langmuir films at the air-water interface. Here, it is suggested that the monolayer formation and growth follows a two-step process, first involving coverage of the Au surface by flat-lying (90° tilt) molecules forming a socalled striped phase followed by a solid-solid phase transition involving nucleation and growth of the tilted phase. The nature of the striped phase was further resolved in a recent study52,54 revealing three structurally distinct intermediates and the emergence of a precursor “liquid” phase just prior to the incipience of the tilted phase. The extensive body of the experimental work, only a small but relevant part of which has been outlined above, provides an excellent background to begin formulating a generic molecular model. It can be expected that such a model, which can accurately depict the formation process in terms of relevant intermolecular interactions, will not only provide a detailed understanding of the surface selfassembly mechanism but also allow design of efficient approaches to equilibrium structures and manipulation of the kinetic pathways to form useful variations of relatively stable, useful intermediate structures. Here, we present a simple coarse-grained molecular simulation model which reproduces the experimentally verified assembly characteristics and optimal structures of an archetypal class of molecular monolayers, namely, alkanethiols on the (111) surface of gold. The model quantitatively accounts for the dominant intermolecular interactions between the constituents of the alkanethiol adsorbate and the gold substrate but excludes the incorporation of the internal conformational flexibility of the alkanethiols. With only four degrees of freedom detailing the Au-S position, molecular tilt, and tailgroup orientation, the model has allowed long-time simulations of up to 105 molecules. The results of the simulations performed at saturation molecular density reveal the emergence of coexisting tilt domains and associated grain-boundary defects in six isoenergetic azimuthal rotations. Each domain was found to consist of hexagonally packed molecular rods at a (30) Poirier, G. E.; Pylant, E. D. Science 1996, 272, 1145. (31) Chailapakul, O.; Sun, L.; Xu, C. J.; Crooks, R. M. J. Am. Chem. Soc. 1993, 115, 12459. Chambliss, D. D.; Wilson, R. J.; Chiang, S. J. Vac. Sci. Technol., B 1991, 9, 933. Dimilla, P. A.; Folkers, J. P.; Biebuyck, H. A.; Harter, R.; Lopez, G. P.; Whitesides, G. M. J. Am. Chem. Soc. 1994, 116, 2225. Dubois, L. H.; Nuzzo, R. G. Annu. Rev. Phys. Chem. 1992, 43, 437. Edinger, K.; Golzhauser, A.; Demota, K.; Woll, C.; Grunze, M. Langmuir 1993, 9, 4. Terrill, R. H.; Tanzer, T. A.; Bohn, P. W. Langmuir 1998, 14, 845.

1476

Langmuir, Vol. 19, No. 5, 2003

uniform molecular tilt of ∼30° from the surface normal. A comparison of the S headgroup position domain structure at the substrate-headgroup interface with the tilt order indicates the discommensurate nature of the domain and defect characteristics between the tilt or tailgroup order and the headgroup positional order. Next, the simulations performed at variable surface densities provide a qualitative phase diagram illustrating the adsorption mechanism. The low-coverage regime, characterized by flat-lying molecules with their molecular axes aligned with the substrate surface, is replaced in the intermediate regime, revealing the coexistence of flatlying and tilted molecules, and the intermediate regime is subsequently replaced by a uniformly tilted monolayer at the saturation coverage. The simulation results were further organized as a primitive two-dimensional isothermal phase diagram for the formation of alkanethiolate monolayers, which provides important insight into the relative roles of intermolecular interactions in the dynamics of the self-assembly process. II. Model Atomic scale molecular dynamics (MD) simulations of alkanethiols on Au(111) have been carried out extensively over the past two decades,28,29,32 verifying many aspects of the packing and phase behavior of self-assembled monolayers. However, such simulations are computationally intensive due to the number of interacting atoms, the complexity of the interatomic potentials, and the time scales of the system dynamics compared to the intrinsic time step of computational MD. Consequently, several groups have explored avenues for simplifying the model systems to more efficiently explore large-scale, long-time behavior of the molecular ensembles. Such simplifications have been done at several different levels from, for example, uniting the CH2 groups of alkane chains into single united atoms32-34 to treating the entire alkanethiol as one35,36 or two36,37 particles with headgroup motion constrained to the gold surface, justified by the 28 kcal/ mol heat of adsorption of methyl sulfide.7 We will here describe a rigid, linear stick model of the alkanethiol. In this approach, the entire molecule is described by four essential degrees of freedom, the headgroup position on the gold surface and two angles describing the relative position of the tailgroup. However, the molecular interactions are described by point-to-point force fields along each molecule. Two major components of molecular properties have been sacrificed in this description. One is the internal, rotational, degree of freedom.28,29 The other is the molecular flexibility, including the swelling of the alkane chains through gauche defects. Justifying the latter simplification, we note that molecular dynamics studies have found only a very small gauche defect content present during simulations at 300 K and below28,29 in agreement with experimental observations.17,38 A linear rigid stick consists of N point particles with positions rji ) (xi, yi, zi)T (1 e i e N, N g 2) as illustrated in Figure 1. The normalized (32) Hautman, J.; Klein, M. L. J. Chem. Phys. 1989, 91, 4994. (33) Siepmann, J. I.; McDonald, I. R. In Computations for the NanoScale; Blo¨chl, P. E., Joachim, C., Fisher, A. J., Eds.; Kluwer Academic: Dordrecht, 1993; pp 49-62. (34) Siepmann, J. I.; McDonald, I. R. Mol. Phys. 1993, 79, 457. (35) Taut, C.; Pertsin, A. J.; Grunze, M. Langmuir 1996, 12, 3481. (36) Grunze, M.; Pertsin, A. J. J. Mol. Catal. A: Chem. 1997, 119, 113. (37) Pertsin, A. J.; Grunze, M. J. Chem. Phys. 1997, 106, 7343. (38) Delamarche, E.; Michel, B.; Kang, H.; Gerber, C. Langmuir 1994, 10, 4103.

Grønbech-Jensen et al.

Figure 1. Rigid model molecule. Particles are constrained to the fixed length, l12, between headgroup, r1, and tailgroup, r2. The headgroup is constrained in the z direction at height h1. The diameter of the jth particle at rj is given by the 6-12 Lennard-Jones parameter, σj. Due to the constraints, a molecule has four degrees of freedom, the x and y coordinates of the headgroup and the two angles, tilt θ and azimuthal rotation φ, of the tailgroup. A model molecule with N ) 5 particles is shown.

and overdamped Brownian dynamics equations of motion are N

rj41 ) hf 1 + s12(rj2 - rj1) + l12

∑ sjk1 + zˆ c1

(1)

k)3

N

rj42 ) hf 2 - s12(rj2 - rj1) + l12

∑ sjk2

(2)

k)3

rj4j ) hf j - l12(sjj1 + sjj2)

(3)

where 2 < j e N, l12 ) | rj2 - rj1|, and zˆ ) (0, 0, 1)T. Characteristic units are defined and discussed below. The head- and tailgroup positions are represented by i ) 1 and i ) 2, respectively, and all other participating objects are represented by i ) j. Forces external to the molecule, including intermolecular interactions and thermal noise, are represented by hfi. The scalars, s12 and c1, and the vectors, sjkl (k ) 3,‚‚‚, N and l ) 1, 2), are Lagrange multipliers ensuring that the following constraints are enforced:

|rj2 - rj1| ) l122

(4)

z1 ) h1

(5)

rjj ) rj1 +

lj1 (rj - rj1) l12 2

(6)

Here, h1 is the constrained headgroup height above the gold surface and lji ) |rji - rjj|, such that l12 ) lj1 + lj2. The Lagrange multipliers, sjki, are determined by the geo⊥ ⊥ metrical relationship lj1 sjj1 ) lj2sjj2 , where sj⊥ji is the component of sjji orthogonal to rj12 ) rj2 - rj1. For any given set of hfi and c1, the above equations can be efficiently solved39 to yield the Lagrange multipliers necessary for advancing time of the dynamical equations of motion of the rigid object. The remaining multiplier, c1, is then determined through an iterative process optimizing the vertical component of hf1 such that the resulting headgroup position is at the desired height, h1. (39) Grønbech-Jensen, N. To be published.


Langmuir, Vol. 19, No. 5, 2003 1477

The forces hfi are given by

hf i ) -∇iU + η j i(t)

(7)

where U is the total potential energy of the system and (y) (z) T η j i(t) ) (η(x) i , ηi , ηi ) is the thermal noise determined by the dissipation-fluctuation relationship,45

h 〈η j (n) i 〉 ) 0 kBT δ δ δ(t - t ) E0 n1n2 i1i2 1 2

〈ηi(n1 1)(t1)ηi(n2 2)(t2)〉 ) 2

(8) (9)

where n ) x, y, z; kB is Boltzmann’s constant; T is the imposed temperature; E0 is the characteristic energy of the normalized potential; and δij and δ(t - t′) are the Kronecker and Dirac delta functions, respectively. Due to the simplified model, constraining all internal structure and degrees of freedom into a cylindrically symmetric rod, we have not included the intertial terms into the equations of motion. Thus, the simulated time is normalized to a characteristic time, τ0, which will be addressed briefly in section III.A, by making direct temporal comparisons with published full MD results of relaxation. The potential energy function consists of three different parts, U ) U1s + Uks + Um, where U1s represents the potential energy of the lateral position of the headgroup relative to the Au(111) surface, Uks (1 < k e N) represents the energy due to the vertical position of the atoms relative to the gold surface, and Um is the energy due to intermolecular interactions. The energies are written

U1s )

∑p u1s(x1, y1)

Uks )

∑p k)2 ∑ uks(zk)

(11)

∑p q>p ∑ um(|rjr - rjs|)

(12)

(10)

N

Um )

where p and q represent the molecules, and rjr and rjs are the positions of all pairs of particles of the different molecules, p and q. The headgroup surface potential, u1s, is given by the energy surface found in ref 21. This energy surface, developed and fitted into an empirical surface potential for use in molecular dynamics of single thiol interactions with gold (111), is a strong function of the Au-S-C angle, which is not represented in our constrained model. We have therefore used the optimized parameters for Dbent as provided in Table 3 of ref 21, e providing a bridge site energy of about 1 kcal/mol for transitions between neighboring face-centered cubic (fcc) and hexagonal close-packed (hcp) hollow sites. We notice that the fcc hollow sites are energetically favored40 over the hcp sites by about 0.8 kcal/mol.42 We have adopted the 12-3 potential for uks given in ref 32 with the amplitude (40) The favored gold (111) fcc hollow site positioning of single thiols has been reconfirmed by more recent extensive DFT calculations in, for example, ref 41. Unfortunately, those publications have erroneously quoted ref 21 as suggesting the hcp hollow site is the energetically preferred location of thiols. (41) Yourdshahyan, Y.; Zhang, H. K.; Rappe, A. M. Phys. Rev. B 2001, 63, 081405. Yourdshahyan, Y.; Rappe, A. M. J. Chem. Phys. 2002, 117, 825. (42) Recent extensive DFT calculations, for example, refs 41, 43, and 44, suggest that thiol dimers may form on gold (111). This possibility is not directly included in the empirical surface potential used in this paper.

Figure 2. Minimized intermolecular energy surface between two adjacent model molecules of given headgroup distances: (a) minimized energy and (b) optimized tilt angle, θ, of the molecular pair. Solid lines represent the energy surface of the parametrized N ) 5, l12 ) 15 model used in this paper. Dashed and dotted lines represent corresponding results using full atomic scale models of dodecane alkanethiols. Three characteristic combinations of full atomic scale model configurations are shown.

of the potential scaled for each particle in order to give the entire molecule the desired surface interaction. For the interparticle potential, um(|rjr - rjs|), we use a smoothly truncated standard 6-12 Lennard-Jones (, σ) potential with the sulfur parameters given in ref 32 for i ) 1. For 1 < i e N, we use parameters determined such that the intermolecular energy reproduces the calculated energy between two atomic alkanethiols in terms of the “hard core” distance and minimum energy. We will attempt to represent a system of dodecane alkanethiol with model parameters l12 ) 15 (Å) and N ) 5 equally spaced particles per rigid rod. For this system, we obtain the following energy parameters: σj ) 4.25 (Å) and j ) 0.8 (kcal/mol) (j > 1) for the interparticle potential and a scaling parameter of 3 for the above-mentioned attractive 3-12 surface-particle interaction for 1 < i e N. These parameters result in a minimized intermolecular (optimally tilted molecules) energy of ∼6-7 (kcal/mol) at optimized distance of ∼4.75 (Å) in good agreement with all atomic models. All interatomic potentials are smoothly truncated at 2σ. Figure 2 shows the optimized intermolecular potential energy surface for a system of two dodecane alkanethiols at given separations. For each distance, the interaction energy is minimized by optimizing the tilt angles, θ1 ) θ2 ) θ for (φ1 ) φ2 ) 0). The result for our l12 ) 15, N ) 5 model is shown by a solid curve in Figure 2a, and the optimized tilt angle, θ, is similarly displayed in Figure 2b. Along with the solid curves representing parametrization of the rigid cylindrically symmetric model, we show that the corresponding energy surfaces form an all-atom model using the MM346 interatomic potential. Due to the atomic scale structure, and the noncylindrical symmetry, of the alkane chains in a molecular model,28,29 we observe that the minimized energy depends on the rotations of the two molecules. Thus, in parametrizing the rigid model, we have attempted to obtain a reasonable representation of the intermolecular energy surface; both the minimized energy and the optimized tilt angle suggest that the model may account reasonably well for the intermolecular behavior. (43) Gro¨nbeck, H.; Curioni, A.; Andreoni, W. J. Am. Chem. Soc. 2000, 122, 3839. (44) Hayashi, T.; Morikawa, Y.; Nozoye, H. J. Chem. Phys. 2001, 114, 7615. (45) See, for example: Parisi, G. Statistical Field Theory; AddisonWesley: Redwood City, CA, 1988. (46) Allinger, N. L.; Yuh, Y. H.; Lii, J. H. J. Am. Chem. Soc 1989, 111, 8551.

1478


Figure 3. Simulated time evolution of average tilt, θ, optimally packed (one thiol per A0 ) 21.6 Å2) samples of M ) 90, M ) 2016, M ) 10 044, and M ) 100 300 molecules at temperature kBT ) 0.6 kcal/mol in periodic boundary condition systems. Also shown is the corresponding evolution of the spatial standard deviation of the tilt. The initial conditions are uniform rectangular lattices of θ ) 0 molecules, randomly positioned in the two surface directions. The results for M ) 2016, M ) 10 044, and M ) 100 300 are nearly identical.

III. Simulation Results


Figure 4. Color map used in Figures 5 and 7. (a) Molecular color depending on azimuthal orientation, φ. Small tilt, θ, is represented by a light color, while larger tilt is shown darker. (b) Two-dimensional headgroup orientation relative to the hexagonal Au(111) surface. Each triangle represents a hollow site, each color edge represents a bridge, and each node is an Au atop site. The thiol-gold potential used differentiates the hollow sites as follows: hollow sites are energetically preferred positions of the thiols, and fcc 3 sites (dark-green, dark-red, and yellow) are energetically preferred over hcp 4 sites (lightgreen, orange, and purple). Optimal packing (one thiol per A0 ) 21.6 Å2) corresponds to one hollow site color being occupied by thiols.

We simulate rectangular systems, periodic in the x and y directions, with dimensions Lx and Ly. The gold surface is positioned at z ) 0, and initial conditions for all simulations are evenly rectangularly distributed headgroups of vertically oriented molecules. We conduct our simulations at T ) 300 K. At this temperature, a normalized time step of 0.005 e dt e 0.025 is adequate to represent the dynamics of the system. A. Optimal Packing Density. This section describes mainly results of molecular simulations for the dynamics of 10 044 molecules in a system dimension of (Lx, Ly) ) (162, 186x3/2)aAu, where aAu ) 2.884 Å. This provides each sulfur with a surface area of A0 ) 21.6 Å2 corresponding to three Au(111) hollow sites, which is the maximum commensurate packing density for sulfur on the (111) surface of gold. Figure 3 shows the average tilt angle, 〈θ〉, and the corresponding standard deviation as a function of simulated time in our model. We have shown results of samples of several different numbers, M, of molecules: M ) 90, 2016, 10 044, and 100 300. We find that the three larger simulated samples behave very similarly in both average tilt and spatial standard deviations thereof. The smallest sample, M ) 90, has a noticeably shorter relaxation time, which is due to the lack of domain walls (see below) between areas of directional tilt in small samples. The plot clearly shows that in only two decades of normalized time, the average chain in the assembly adopts the stable molecular tilt of 28° from its initial vertical orientation. The observed tilt behavior is further characterized by an extremely low standard deviation of only a few degrees for the molecular chains comprising the ensemble. Taken together, these results suggest the equilibrium (or frozen) monolayer structure to be characterized by a stable and single uniform molecular tilt of 28°, in excellent agreement with a previous pseudo-atom molecular dynamics simulation model where the average chain tilt of 28° was reported.28,47 Experimental support for this observation is also overwhelming. A variety of techniques, including ellipsometry,27 X-ray,15 SFG,25 infrared22 and Raman vibrational spectroscopy,23 and more recently STM,19 have all yielded average molecular chain tilts in a narrow range of 24-34° for the highest quality films. Such films,

obtained by extended substrate immersion (periods from 1 h to several days) in the self-assembly solution, yield limiting film coverages corresponding to 21.6 Å2 per molecule, the coverage used in the present simulations. Further, note the rapid and well-characterized approach to a stable equilibrium tilt angle in the simulation. We interpret this feature to be consistent with the notion that the initial discommensurate configuration of the molecular assembly, characterized by vertically oriented chains in a square lattice, very quickly adopts the stable tilt configuration commensurate with the Au lattice (see below). Comparing these data to MD simulations of the relaxation of optimally packed alkanethiols provides a rough estimate of the inherent time scale, τ0, of the simplified simulations. For example, typical relaxation times reported in MD studies48 of M ) 90 alkanethiols on gold at optimal packing, where the initial condition was also θ ) 0 as in this study, are about 50 ps for obtaining an optimally packed monolayer at 300 K. Making direct comparison between the studies of ref 48 and our simulations for M ) 90 displayed in Figure 3 provides a rough estimate of the normalized time scale in our model during the approach of forming a uniform tilt of θ ≈ 30°. We estimate the time scale to be τ0 ≈ 1 ps. Next, we followed the in-plane ordering in the monolayer structure. This can be most conveniently depicted by monitoring the dynamics of the azimuthal directional angle of the molecular axes in our model. To efficiently visualize this azimuthal lateral orientation, φ, of the molecules (Figure 1), we distinguish the molecules by their preference for the lateral orientation as depicted by the color map shown in the left panel of Figure 4. The continuous variation in the color scheme, based on the azimuthal orientation and tilt angles, allows capturing the full range of chain rotational and tilt angles in a single image. To visualize the headgroup position relative to the gold hollow sites, we use a hexagonal color map shown in the right panel of Figure 4. Each triangle in the lattice represents a single hollow site in the (111) surface of gold. The 6-fold symmetry in this depiction characterizes the

(47) Tarek, M.; Kechuan, T.; Klein, M. L.; Tobias, D. J. Biophys. J. 1999, 77, 964.

(48) Sprik, M.; Delamarche, E.; Michel, B.; Ro¨thlisberger, U.; Klein, M. L.; Wolf, H.; Ringsdorf, H. Langmuir 1994, 10, 4116.


two types of near-isoenergetic patterns of hollow sites offered by the Au lattice for S occupancy in the x3 × x3 R30° superlattice. It is further important to note that in the “rod” model, the S positions are intimately linked to tailgroup or (methyl group) position, offset only by the tilt of the molecular chains. Since the molecular tilt is found to be nearly uniform (see above) in our samples, it is expected that the methyl group ordering will follow. That this correspondence is often broken at room and elevated temperatures due to conformational unlocking of the terminal methyl groups has been suggested using He diffraction17 and infrared spectroscopy measurements.22 With this limitation of the model in mind, we do not attempt to characterize the in-plane order of the methyl tailgroups. A second limitation of the model arises from the use of the cylindrically symmetric chain approximation. As a result of this approximation, chain structural variations and domain structures, that result from precession of the C-C-C backbone plane about the chain molecular axis, cannot be captured in the model as it is described here. Figure 5 shows time-lapse snapshots at four different times representing the dynamic or “evolutionary” period for the in-plane ordering in our molecular rod ensemble. The left panels represent the progression of tilt-rotation boundaries, and the right panels the S positional order relative to the 3-fold hollow sites. A casual glance at the images shows that even at very early times in the molecular simulation, domains of preferred tilt and azimuthal angles begin to emerge with well-defined domain boundaries at the molecular scale. Corresponding images on the right panels reflect similar evolution of the S positional order. It is obvious from the headgroup position domains that all six colors of hollow sites are represented initially, while later times clearly select domains of the three lower energy, fcc 3, sites (see Figure 4). Initially, both tilt and S positional domains are small, and they rapidly coalesce into larger domains as seen in Figure 5a-d. Along with the molecular conformation pictures, Figure 5 displays the distribution of azimuthal angles in the sample. We can here follow the tilt domain growth as a strong selection of the density of tilt orientation with selective preference for the NNN orientations, φ ) n60°. At longer simulation times, the domain patterns become quasi-static, presumably indicating an approach to the final (glassy) or equilibrium packing. That an approach to the optimal structure occurs via domain growth and coalescence at the highest coverage has also been suggested in many experimental studies. For example, Delamarche and co-workers have observed, using high-resolution scanning tunneling microscopy, that asdeposited monolayers possess a large number of coexisting stable domains which upon annealing result in fewer, larger domains.49 It is here important to recall that despite close similarity with the domain structure obtained by Delamarche et al. and Anselmetti et al.,50 domain boundaries obtained in our study are different. The domains observed due to a slight height mismatch (below 0.9 Å) in the STM measurements49,50 were attributed to the differences in twists of the molecular chains about the molecular axis (rotation of the C-C-C backbone plane from the C-C-C axis), whereas the domains in the present study reveal grain boundaries that emerge due to the mismatch in the azimuthal rotation of the molecular rod (49) Delamarche, E.; Michel, B.; Gerber, Ch.; Anselmetti, D.; Gu¨ntherodt, H.-J.; Wolf, H.; Ringsdorf, H. Langmuir 1994, 10, 2869. (50) Anselmetti, D.; Baratoff, A.; Guntherodt, H. J.; Delamarche, E.; Michel, B.; Gerber, Ch.; Kang, H.; Wolf, H.; Ringsdorf, H. Europhys. Lett. 1994, 27, 365.

Langmuir, Vol. 19, No. 5, 2003 1479

from an arbitrarily designated initial structure. Also, topographic measurements in the STM cannot straightforwardly resolve the domain walls based on azimuthal rotational angle discrepancies measured here. A closer inspection of the images shows a number of attributes that can be ascribed to the level of the selfassembly process itself and captured in the rigid model. First, we note that the azimuthal and tilt structures obtained for the longest simulation time (left panel in Figure 5d) reveal large domains represented dominantly by a single yellow-green texture in the figure. This orientational state of the chain in our model corresponds to the molecular chain tilted at 28° and pointing approximately in the direction of the next-nearest neighbor. This observation is highly consistent with the previous experimental and theoretical studies that have revealed the preference for the NNN azimuthal direction for the molecular chain rotation.17,29 Second, comparing the azimuthal images for the four snapshots above, we note that the patches of differently oriented domains (e.g., a purple-coded rectangular domain in the center of the image) appear to form quickly and remain trapped in the otherwise annealed matrix. Thus, the formation of longlived patches of rotational defects appears to be intrinsic to the assembly dynamics. This is due to the intrinsic frustration of the (6-fold) local degeneracy of the preferred energetically identical φ directions (NNN at n60°) when the headgroups are hexagonally packed in the gold hollow sites. Third, we note distinct point defects or “holes” in the surface coverage. These holes are typically a couple of molecular spacings (1 nm) in range and also seem intrinsic to the packing behavior in the monolayer. The presence of similar holes has been observed in a number of recent high-resolution STM measurements.13,51,52 Fourth, comparing the images for the chain orientation in the left panels with those for S positional order confirms that domain formation of sulfur positioning in specific hollow sites takes place together with the domain formation of the molecular orientation. But the domain boundaries obtained for the S positional order do not appear to correlate strongly with the domain structure in the orientational picture. In fact, the sulfur positional order of domains, with sulfur occupying every third hollow site, gives rise to the very specific preferred directions of the molecular tilt (toward the next-nearest neighbor, φ ) n60°). This apparent discrepancy can be understood if one notes that the six (three fcc) competing patterns for S occupancy do not have direct impact on the azimuthal preference for the chains. In other words, differing domains of the S lattices result in the identical preferred azimuthal orientation for the molecular chains, since the three preferred sulfur domain types are simple translations of each other. As a result, the chain assembly intrinsically heals the defects created at the S sites due to occupancy differences. This observation is quite curious and can be corroborated by experimental studies jointly utilizing structure tools that are sensitive to microscopic descriptions of chain packing (e.g., vibrational spectroscopies) and the S-Au interface. To the best of our knowledge, such studies for the characterization of alkanethiols on Au have not yet been performed. To quantify the visual impressions gleaned from Figure 5 above, we have conducted several distinct simulations and gathered statistics on domain formation and defect content as a function of simulation time. Figure 6 shows (51) Stranick, S. J.; Parikh, A. N.; Allara, D. L.; Weiss, P. S. J. Phys. Chem. 1994, 43, 11136. (52) Poirier, G. E.; Fitts, W. P.; White, J. M. Langmuir 2001, 17, 1176.

1480



Figure 5. Simulated time evolution of a 10 044 molecule sample at optimal packing density, at kBT ) 0.6 kcal/mol, and for periodic boundary conditions in the two surface directions. The initial condition is a randomly positioned rectangular lattice of thiols with θ ) 0. Snapshots of azimuthal orientation (top left), headgroup position (top right), and azimuthal angle density (bottom) are shown for (a) time ) 25τ0, (b) time ) 125τ0, (c) time ) 1000τ0, and (d) time ) 5000τ0. The color map for the top views of the sample is shown in Figure 4.

the time evolution of average domain sizes for both orientational and S positional properties, measured as the number of participating molecules. The shown data are composed from eight simulations of 10 044 molecules. We observe a rapid saturation of average domains at about half the system size for the tilt direction and at about 103 molecules for the headgroup ordering. It is important to note here that the evolution periods for the chain orientational and S positional properties are comparable,

suggesting that both order parameters evolve concomitantly during this phase and then saturate to a stable and presumably optimal (or long-lived “glassy”) size. Since the domains of azimuthal angle seem to saturate at a size of the same order as the simulated system size, we have conducted additional simulations of 25 000, 50 000, and 100 000 rods to determine if the domain characteristics depend on the number of degrees of freedom. So far, we have determined that finite size effects may be somewhat


Figure 6. Domain formation as a function of simulated time for eight samples of 10 044 molecules at optimal packing density and at kBT ) 0.6 kcal/mol. Average size domains of azimuthal orientation and headgroup positions are shown. Domain sizes are measured in number of participating molecules in a domain.

important when simulating 10 044 rods, indicating that we need to conduct long-time simulations of at least 50 000 rods in order to adequately account for the tilt domain evolution. However, the sulfur position domains and statistics show little or no sensitivity to system size for simulations larger than 10 044 rods, suggesting that the headgroup domains saturate at about 10 nm size. We notice that even if the domains of average tilt have relatively strong finite size signatures, the defect domain patches of sizes of ∼5 nm seem to prevail in all system size simulations, indicating that these systems have an intrinsic propensity for generating sharp edge domains at the nanometer scale in the range of ∼5-10 nm in both tilt direction and headgroup positions. The significant difference in size between tilt domains and headgroup positions is not surprising given the uncorrelated domain boundaries of the two types of domains. One may rationalize the relatively large size of the tilt domains in terms of the relatively large van der Waals penalty of a tilt domain wall compared to the small energy penalty of a simple headgroup translation, which is the signature of a headgroup domain wall. B. Formation Mechanism and Emergent Structures at Submonolayer Coverages. This section describes results of molecular simulations for the dynamics of alkanethiols in a system of surface area of ∼1743 nm2, which can accommodate M ) 8064 molecules at optimal packing density. At room temperature, we simulate various molecular surface densities in order to construct a simple model phase diagram, in the spirit of, for example, refs 36 and 37, validated by a number of recent experimental observations30,52-54 for monolayer growth. To capture the essential mechanistic pathway proposed for the formation of alkanethiol monolayers on gold,30 we extended our simulations of rigid rods to include variable suboptimal molecular packing densities. We conducted simulations of a surface area of A ∼ 1743 nm2 with a varying number of molecules to assess the model properties as a function of molecular density. For all packing densities, we chose the initial configuration comparable to that used for the optimal coverage above: evenly distributed pseudo-molecular rods in a rectangular lattice with vertical orientation, or zero tilt (θ ) 0). The color scheme used to depict the molecular orientation and S positional order was also the same as that used for optimal coverage simulations. Figure 7 shows the orientational and S positional snapshots together with the distribution (53) Poirier, G. E. Langmuir 1999, 15, 1167. (54) Fitts, W. P.; White, J. M.; Poirier, G. E. Langmuir 2002, 18, 1561.

Langmuir, Vol. 19, No. 5, 2003 1481

of tilt, θ, for four different packing densities for a normalized simulation time of t ) 4000, which we estimate to be a reasonable evolution period. The packing densities used for Figure 7a-f were chosen between the optimal coverage of one thiol per A0 ) 21.6 Å2 (coverage ) 1) and dilute coverage of one thiol per 8A0 (coverage ) 0.125). These simulations provide several noteworthy features that characterize the trends in evolution of orientational and S positional preference as a function of molecular coverage. At very low molecular densities (Figure 7f), the molecular rods align along the substrate surface. The image shown for 1/8 of the optimal coverage clearly shows the presence of flat-lying molecular rods distributed on the (111) surface of gold. Even at these low molecular densities, patterns of closely packed molecules are evident as also seen in ref 37. The patterns can be understood in terms of the preference of molecular chains to pack due to intermolecular van der Waals interactions. Images of the S positional order accordingly show an alignment of S along a linear dimension consistent with the orientational picture. As the coverage increased, this scenario persisted until the entire surface of gold was covered with surface-aligned molecules at about 1/5 (Figure 7e) of the optimal coverage leading to a “monolayer” characterized by a uniform 90° tilt but revealing no significant correlations in the rotational azimuthal angle or the S occupancy at the Au surface. Beyond this packing, at the images for 1 /2 of optimal coverage (Figure 7d), a new phase is observed to nucleate as islands. S positional order at this coverage also begins to show 2D domains in contrast to 1D ordering that is seen for the surface-aligned phase. This nucleating phase shows a number of key characteristic properties. They represent a number of small islands characterized by molecules in single orientational preference. The latter preference is invariably a uniform tilt of 28-30° and varies only in the azimuthal rotational angles. Molecules within a contiguous domain almost always display an identical azimuthal angle. This suggests an abrupt formation of a uniformly tilted phase displaying features reminiscent of those of optimal monolayers. Note also a direct correspondence between the sizes of S positional order and the orientational order in the nucleating phase. Domains of molecules with lifted tailgroups also display a single habit for S occupancy. At even higher coverages of 5/8 and 3 /4 of optimal density (Figure 7b,c), we begin to note a greater fraction of the molecules in the 30° tilted phase with larger tilted phase domains, but the registry between the S positional order and the orientational order domains is beginning to weaken, suggesting that the S positional order reaches its critical in-plane correlation length at a much lower value than the tilt correlation. As discussed earlier, the loss of S translational order does not significantly penalize the coherence in the tilt domains (see above) and this lack of correspondence is not inconsistent with the general domain growth behavior captured here. At this time, we do not fully understand what limits the S positional translational order, but it appears to be a consequence of high nucleation rate for the tilted phase and the isoenergetic nature of the competing S positions in the 3-fold hollow patterns. Further studies of the suboptimal packing shown here are being conducted to determine the statistics of the pattern at much longer time scales. The above observations lend themselves to a very simple binary model. It is clear from the images that both the optimally covered system and the dilute one have rather uniform molecular tilts throughout the systems. To analyze the intermediate configurations, we have accompanied the snapshots in Figure 7 with statistical

1482


averages of the probability, f(θ), of finding molecules with tilt θ. These data are acquired as an average over 500 time steps preceding time ) 4000τ0. What is emerging from this density of tilt, f(θ), is a rough confirmation of the binary state theory, where molecules, for coverages between optimal and dilute, form coexisting domains of optimal packing (θ ≈ 30°) with domains of maximum surface contact (θ ≈ 90°). The molecular content in these domains can then be estimated by assuming that the free energy is minimized for maximum molecular contact with


the surface. Elaborating on and quantifying the simple binary model, we write the average tilt of a sample in the following manner:

〈θ〉 )

1

M

∑ θm ) csθs + ctθt Mm)1

) θs + ct(θt - θs)

(13) (14)

where M is the total number of molecules in the sample,


Langmuir, Vol. 19, No. 5, 2003 1483

Figure 7. Simulated samples of area A ) 1741.824 nm2 with different molecular coverages. The initial conditions are uniform rectangular lattices of θ ) 0 molecules, randomly positioned in the two surface directions. The simulated temperature is kBT ) 0.6 kcal/mol. Snapshots are shown for a simulated time of 4000τ0. Snapshots of azimuthal orientation (top left), headgroup position (top right), and tilt angle density (bottom) are shown for the following number, M, of molecules: (a) M ) 8064 (optimal packing, A0/A ) 1); (b) M ) 6048 (suboptimal packing, A0/A ) 3/4); (c) M ) 5040 (suboptimal packing, A0/A ) 5/8); (d) M ) 4032 (suboptimal packing, A0/A ) 1/2); (e) M ) 1612 (dilute packing, A0/A ≈ 1/5); (f) M ) 1008 (dilute packing, A0/A ≈ 1/8). The color map for the top views of the sample is shown in Figure 4.

and cs and ct are the fractions of molecules with tilts θs and θt, respectively, such that cs + ct ) 1. We can choose θs ≈ 30° and θt ≈ 90°. The standard deviation of the tilt in such an ensemble is given by

δθ ) x〈θ2〉 - 〈θ〉2 ) xct(1 - ct)|θs - θt|

(15) (16)

We will assume that the relative content, ct, of molecules at tilt θt is determined by a requirement of maximizing the coverage of the gold surface. The area (per molecule) covered by the two states is

A ) csA0 + ctAf 1ct )

A A0

Af -1 A0

(17)

(18)

where A0 ) 21.6 Å2 is the area covered by the θs ≈ 30° molecules, and Af is the surface contact area of a molecule at tilt θt ≈ 90°. Working with this model, we will assume A0 e A e Af. Since the characteristic surface coverage area, Af, of a θt ≈ 90° molecule is roughly Af ≈ (l12 + σ)σ, the above eqs 14, 16, and 18 provide an average tilt, as well as the spatial fluctuations in the domains, as a function of the surface coverage, A0/A. Figure 8 shows the statistics of the tilts, 〈θ〉 and δθ, as a function of surface coverage, A/A0. The solid lines in the interval A0 e A e Af show the result of the above binary

Figure 8. Average molecular sample tilt, 〈θ〉, and corresponding spatial standard deviation, δθ, as a function of molecular surface density. Solid lines represent the binary theory of eqs 14, 16, and 18. Simulated samples have total area A ) 1741.824 nm2, temperature kBT ) 0.6 kcal/mol, and data are taken at simulation time 4000τ0. The initial conditions are as described in Figure 7. Four characteristic phases are identified. (I) An unrealistically dense phase with average tilt 〈θ〉 ) 0. (II) A high-density phase, transitioning from 〈θ〉 ) 0 to optimal coverage 〈θ〉 ) θs ≈ 30°. (III) The mixed domain, or coexistence phase, of θs ≈ 30° and θt ≈ 90°. And finally a phase (IV) in which all molecules are flat, θt ≈ 90°, and the gold surface cannot be completely covered by the alkanethiols.

theory, while the dots show the corresponding simulated data. The figure provides four distinct phases: (I) Ultrahigh-density phase with θ ) 0 (0 < A < Av). This phase is unrealistically dense. The characteristic area per thiol, Av, is strongly dependent on model (molecule) details. (II)

1484


Superoptimal density (Av < A < A0). This phase is characterized by 0 < 〈θ〉(

Rigid Molecular Model for the Assembly Characteristics and Optimal

Recommend Documents