Langmuir 1986,2, 219-227
219
Computer Simulation Study of the Multilayer Adsorption of Fluid N2 on Graphite Alexei V. Vernovt and William A. Steele* Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802 Received August 20, 1985
-
Computer simulations of nitrogen adsorbed on the graphite basal plane at 73.6 K are reported for coverages ranging from 1 monolayer to -2.5 layers. Quantities evaluated include average potential energies and the component molecule-solid and molecule-molecule parts of these energies, translational and orientational order parameters of Nz molecules at the surface, and molecular pair correlation functions projected onto planes parallel to the surface. It is shown that both the orientations and the density of molecules in the first layer change as the total coverage changes, that there is a definite preference for the commensurate (or ( d 3 X d 3 ) R 3 0 ° ) translational ordered structure for first-layer molecules, and that the presence of a dense second layer can actually cause a slight decrease in the number of molecules held in the first layer. Comparisons with available experimental data are good, indicating that the simulations are based on a realistic model of the actual N2-graphite system. gassolid forces, so that molecules transfer from one layer 1. Introduction to another. Thus one eventually obtains values for the Recently, a considerable amount of new information has equilibrium layer densities as well as for the fluctuations been produced concerning the thermodynamic and about the equilibrium values. Although there is almost structural properties of nitrogen film formed on the excomplete orientational disorder for these molecules in the posed basal plane of graphite.l-19 The techniques used plane of the surface, calculations of the angles between the in these investigations include thermodynamic measuremolecular axes and the surface plane show considerable ments of heats and isotherms at temperatures around the orientational order in this degree of freedom. Furthermore, boiling point (78 K)'r9J0 and of heat capacities a t lower order, as measured by the positions of the temperatures,2J1J2 diffraction of X-rays,lg n e u t r o n ~ , ~ ~ ~ translational J~ molecular centers-of-mass relative to the carbon atom and electron^^^^-^ by the ordered structures that form in hexagons, was characterized and found to be significant, the monolayer and multilayer films at temperatures below at least in the first layer of adsorbed molecules. We note -50 K, and computer simulations of model N2-graphite that the interaction energies and effective size for a N2 It is now clear that the monolayer film can molecule rotating freely in the plane of the surface are occur in one of several different solid structures, starting similar to those for krypton molecules, so that our results with the ( ~ ' 3 x 4 3 commensurate ) system wherein the N2 and conclusions are of interest not only for the N2-graphite molecules adsorb in a regular array over the centers of system but have implications also for the complicated one-third of the carbon hexagons of the graphite basal phase behavior observed for the analogous krypton sysplane and, progressing as density increases, through at least tem.2@22 two other ordered structures. A t moderate coverage, the commensurate film melts directly to a supercritical fluid phase a t T N 47 K. If two-dimensional liquid exists at (1) Larher. Y.J. Chem. Phvs. 1978.68. 2257. all in the monolayer, it is present only in an extremely (2j Chung,'T. T.; Dash, J. G. Surf. hi.' 1977,66, 559. (3)Fain, S. C., Jr.; Toney, M. F.; Diehl, R. D. In "Proceedings of the limited region of the p , T phase diagram for this system. Ninth International Vacuum Congress and Fifth International ConferAt temperatures less than 80 K the fluid N2 phase which ence on Solid Surfaces": de Seaovia, . J. L.,. Ed.:.~ImDrenta Moderna: forms at low coverage appears to transform into an ordered Madrid. 1983: DD 129-137. (4)Kjems, J: K.;Passell, L.; Taub, H.; Dash, J. G.; Novaco, A. D. Phys. phase over a narrow range of density near monolayer Reu. B 1976,13, 1446. completion, but the interpretation of the experimental data (5) Eckert, J.; Ellenson, W. D.; Hastings, J. B.; Passell, L. Phys. Rev. is complicated by the onset of second-layer adsorption in Lett. 1979,43, 1329. (6)Diehl, R. D.; Toney, J. F.; Fain, S. C., Jr. Phys. Rev. Lett. 1982,48, this region. Thus, it appears that the detailed molecular 177. information obtainable from computer simulation of a (7) Diehl, R. D.; Fain, S. C., Jr. J . Chem. Phys. 1982, 77,5065. realistic model of this system might be very helpful in (8) Diehl, R. D., Fain, S. C., Jr., Phys. Reu. B 1982,26,4785. generating information concerning the molecular transla(9)Piper, J.; Morrison, J. A.; Peters, C.; Ozaki, Y. J . Chem. SOC. Faraday Trans. 1983, 79,2863. tional and orientational order. Indeed, simulations at 73.6 (10)Chan, M. H.W.; Migone, A. D.; Miner, K. D., Jr.; Li, Z. R. Phys. K and a t densities ranging up to the commensurate value Rev. E 1984,30,2681. have already been reported.16 In the present paper, sim(11)Miner, K. D., Jr.; Chan, M. H. W.; Migone, A. D. Phys. Reu. Lett. 1983,51,1965. ulations at this temperature are reported in which the (12)Zhang, Q.M.; Kim, H. K.; Chan, M. H. W. Phys. Reu. B 1985,32, coverage of nitrogen on the surface is extended from 1985. commensurate up to 2.5 commensurate layers. (13) Peters, C.; Klein, M. L. Mol. Phys. 1985,54, 895. (14)Talbot, J.; Tildesley, D. J.; Steele, W. A.Mol. Phys. 1984,51,1331. It is found that the molecules can easily be assigned to (15)Talbot, J.;Tildesley, D. J.; Steele, W. A. Surf. Sci., in press. Joshi, specific layers, so that one can evaluate thermodynamic J. P.; Tildesley, D. J. Mol. Phys. 1985,55,999. and structural properties for the first, second, and third (16)Talbot, J.: Tildesley, D. J.: Steele, W. A. Faraday Discuss. Chem. SOC.,in press. layers in these systems. The molecular dynamics tech(17) Migone, A. D.; Kim, H. K.; Chan, M. H. W.; Talbot, J.; Tildesley, nique employed puts no restrictions on the molecular D. J.: Steele. W. A. Phvs. Reu. Lett. 1983.51. 192. motion other than those due to the intermolecular and (18) Nielsen, M.; KGr, K.; Bohr, J. J. hkctron Spectrosc. Rel. Phe-
-
Permanent address: Department of Physical and Colloid Chemistry, People's Friendship University, 117302 Moscow, USSR.
nom. 1983,30, 111. (19)Morishige, K.; Mowforth, C.; Thomas, R. K. Surf. Sci., in press. (20) Koch, S. W.; Rudge, W. E.; Abraham, F. F. Surf. Sci. 1984,145, 329.
0743-7463/86/2402-0219$01.50/0 0 1986 American Chemical Society
220 Langmuir, Vol. 2, No. 2, 1986
Vernov and Steele
namics algorithm25(canonical ensemble, to order 1 / P ) ; to our knowledge, this is the first usage of an algorithm of this kind in simulations of systems with large inhomogeneities due, in this case, to the rapidly varying gas-solid potential energy function. Results obtained with the constant temperature algorithm will be compared with those simulated a t const.-it energy and it will be argued that the constant-temy . r u e algorithm is an important improvement in the mr ~culardynamics approach to these systems. (Of course, the previous Monte Carlo multilayer simulations were performed in a canonical ensemblegrand canonical, to be specific-which is normal for that technique.)
0
0
'I
co
I
r
CD
A portion of the currently accepted phase diagram for N2 films on graphite in the coverage-temperature plane is shown in Figure 1. This figure shows the coverage values chosen for study in the present work. It is evident that the results reported in this paper have the potential for producing detailed new information concerning the technically and theoretically interesting question of how simple nonspherical molecules behave as one progresses from monolayer- to multilayer-adsorbed films on a wellcharacterized solid substrate. Previous computer work in this range of coverage has been limited to a few Monte Carlo studies of (spherical) rare gas atom^.^^,^^ Although these pioneering studies were of great utility in revealing the nature of the layer-by-layer growth in multilayer adsorption near the bulk boiling points of the adsorbates, they suffer from two defects relative to the present simulations: it is impossible to extract information from a Monte Carlo simulation concerning the dynamics of molecular translation and rotation in the adsorbed layer, and, of course, questions concerning molecular orientational changes with coverage could not be addressed. Note that one of the points of particular interest in the present work is translational dynamics as related to the rate of transfer from one adsorbed layer to another and the influence of such density fluctuations upon other thermodynamic properties. In section 2, the model system taken for study will be described, together with simulation techniques used. Because of the large potential energy changes associated with interlayer transfer, the usual microcanonical (or constant energy) ensemble was discarded in favor of a constant-temperature molecular dy-
2. Potential Model and Simulation Algorithm In order to maintain close contact with the simulations of monolayer adsorption of N2on the graphite basal plane a t T 73.6 K, the potential energies of interaction between pairs of N2 molecules27and between an Nz and the surface were taken to be identical with those used in the previous study.16 To recapitulate briefly, the N2-N2energy is assumed to contain Lennard-Jones (W) and electrostatic terms, each of which is represented as a sum over discrete sites in the molecule. Two LJ sites were located at the N atoms, so that the N2-N2 LJ interaction was given by the sum of the four site-site inverse 12-6 potential functions with well depth t"/k = 36.4 K and size U" = 0.332 nm. For calculating the electrostatic N2-Nz energy, charges of -6.49 X C were placed a t the atomic sites together with a charge of 12.98 X lo-%C at the center of symmetry of the molecule. The N2solid energy2sis also assumed to be a sum of site-site energies: each N atom interacts with each C atom in the solid via a LJ 12-6 function with tNC/k = 31.92 K and aNc/k = 0.336 nm. (For computational speed, it is convenient to express each of the two sums over the solid atoms as a Fourier series in the reciprocal lattice vectors of the surface lattice-details of this procedure are now standard29and need not be repeated here.) The usual molecular dynamics algorithm for an ensemble of nonspherical molecules consists in numerical solution of Newton's equations of motion for a sample of a few hundred molecules. The sample is surrounded by replicas of itself, so that molecules can pass through the walls of the computer "box" while conserving number N , volume V, and total energy E of the sample. In physisorption simulations, the gas-solid interaction enters as a (timeindependent) external potential. Periodic replication of the sample is done in the x,y directions parallel to the surface, and in the present work, possible loss of molecules by translation away from the surface is handled by reflection at a conveniently distant plane from the surface. Equations of motion for orientational degrees of freedom are often expressed in terms of quaternions rather than Euler angles in order to avoid numerical difficulties associated with finite arithmetic in Euler-angle space.30 This procedure was followed here. A timestep of length 2.022 x ps gave adequate energy conservation (to a few parts in lo5) in the solutions to the equations of motion for samples of molecules ranging from 96 to 240, depending upon coverage. The fifth-order predictor-corrector al-
(21) Thomy, A.; Duval, X.; Regnier, J. Surf. Sci. Rep. 1981, 1 , 1. Vilches, 0. E. Annu. Reu. Phys. Chem. 1980, 31, 463. (22) McTague, J. P.; Als-Niesen, J.; Bohr, J.; Nielsen, M. Phys. Reu. B 1982,25, 7765. Butler, D. M.; Litzinger, J. A.; Stewart, G. A,; Griffiths, R. B. Phys. Reu. Lett. 1979, 42, 1289. (23) Whitehouse, J. S.; Nicholson, D.; Parsonage, N. G. Mol. Phys. 1983, 49, 829. (24) Rowley, L. A.; Nicholson, D.; Parsonage, N. G. Mol. Phys. 1976, 31, 365, 389.
(25) Evans, D. J.; Morris, G. P. Computer Phys. Rep. 1984, I , 297; Chem. Phys. Lett. 1983,98A,433. (26) Nos&,S.J. Chem. Phys. 1984,81,511. (27) Murthy, C. S.;Singer, K.; Klein, M. L.; McDonald, I. R. Mol. Phys. 1980,41, 1387. (28) Steele, W. A. J. Phys. (Paris) 1977, C4, 61. (29) Steele, W. A. Surf. Sci. 1973, 36, 317. (30) Evans, D. J. Mol. Phys. 1977, 34, 317. Evans, D. J.; Murad, S. Mol. Phys. 1977, 34, 327.
i I 0
1
1 1 1--I ' 1 1
1
1
IO 20 30 40 50 60 70 80 Temperature
Figure 1. Phase diagram for nitrogen on graphite, based on current experimental information. Monolayer phases: TID, which is triangular, incommensurate, orientationally disordered (inplane); UIO, uniaxial, incommensurate, ordered; CO and CD, the commensurate, orientationally ordered and disordered phases. T h e positions of some of the lines drawn are still speculative, especially a t high coverage and high temperature. Densities and the temperature investigated in this simulation are indicated by the points.
Langmuir, Vol. 2, No. 2, 1986 221
Computer Study of Multilayer Adsorption
Table I. Simulations of Multilayers of N2on Graphite at 73.6 K no. of molecules in run simulation box 1 96 2 100 3 106 4 120 5 130 6 144 7 168 8 192 9 240
size of simulation box p*
0.385 0.4 0.425 0.481 0.521 0.575 0.674 0.77 0.962
N,* * L,* 1 1.04 1.1 1.25 1.35 1.5 1.75 2.0 2.5
18 18 18 16 16 18.9 16 18 18
L,*/31/2 8 8 8 9 9 7.65 9 8 8
starting config herringbone herringbone herringbone triangular triangular triangular triangular herringbone last config from run 8 plus 48 molecules added
no. of timesteps in no. of timesteps additional no. microcanonical ensemble in NVT ensemble of timesteps in after equilibration after equilibration NVT ensemble' 8800 9775 14990 15925 17075 7875 5000 10925 19475 7675 3800 7775 18125 16100 6375 8700 20200 5550 23575
a p * = density in reduced units = (molecules/A2)(2.46 A)2. *N* = number of molecules/number in a commensurate layer. cUsed for calculation of mean values.
gorithm employed is identical with that used previously and is described in the earlier papers. Numbers of molecules, dimensions of the simulation box, and other parametersof interest are listed in Table I for the simulation runs reported here. When the conditions were chosen such that significant multilayer formation occurred (essentially, in all but the first two runs listed), it was soon noticed that transfer of a molecule from the first layer to any other and vice versa gave rise to difficulties in a constant energy (microcanonical) simulation. The problem arises because such molecules experience a large change in gas-solid potential; conservation in total energy requires an accompanying change (in the opposite direction) in average kinetic energy and, consequently, a significant jump in the instantaneous temperature, which is obtained from the equation ( K E ) = (5/2)NkT (KE = kinetic energy). Since interlayer transfers occur relatively infrequently (see below), one is attempting to evaluate thermodynamic averages in a system exhibiting long-lived, large temperature fluctuations. In order to eliminate this, a constant-temperature algorithm due to Evans et al.25was adapted to this system of nonspherical molecules and utilized for much of the simulation work. In this algorithm, the equations of motion are solved subject to the constraint that the total kinetic energy of the sample remain constant at'all times. The method of Gaussian least curvature yields modified equations of motion which were numerically solved, giving kinetic energy (and thus, temperature) which fluctuates by a few parts in lo5 from timestep to timestep. From a statistical mechanical point of view, has shown that such a system yields equilibrium averages that are nearly the same as those for a canonical ensemble (differences of order 1/N) and Evans31 has shown that the dynamical properties of this constant-temperature system are also nearly identical with those expected for a canonical ensemble. Thus, the procedure followed was to select a reasonable starting configwation for the N2 molecules near the surface (briefly indicated in Table I) and a random choice of initial velocities such that the total kinetic energy was close to the desired value. The system was then allowed to equilibrate at the desired temperature for 3000 timesteps. The subsequent time evolution of the molecular coordinates and velocities was used to calculate averages of thermodynamic and local order parameters. The number (31) Evans, D. J.; Morriss, G . P.Chem. Phys. 1984, 87, 451.
of timesteps in each of these runs is indicated in Table I, where it can be seen that both the constant total energy and the constant kinetic energy algorithms were employed for various time intervals at the different surface coverages chosen for study. The starting configurations were taken either to be herringbone structures which contained 96 molecules with axes parallel to the surface in the first layer at a height of 0.332 nm and extra molecules with the same coordinates and orientation in the second layer at the height 0.6765 nm or a two-dimensional triangular lattice with 120 molecules with N-N bonds perpendicular to the surface in the first layer at a center of mass height of 0.344 nm and extra molecules in the second layer at the height 0.672 nm with the same perpendicular orientation to the surface. In Table I and later we will use reduced units of length, energy, temperature and time, defined as follows: length, I* = l / a ; energy, E* = E/cNC; temperature, T* = kT/cNc; time, t* = t(cNC/Ma2)'J2, where M is the molecular mass of nitrogen and a is a surface lattice spacing parameter equal to 0.246 nm. In all cases except the calculation for 144 molecules a simulation box with area 15.09 nm2 was utilized. The surface lattice for such a cell has the correct periodicity to ensure continuous periodicity in the underlying graphite. One of the interesting and unexpected findings of these simulations was that careful inspection of time-dependent average quantities obtained after the nominal equilibration period of 3000 timesteps revealed that not only did these averages exhibit the expected fluctuations but that long time drifts were present as well. The two primary quantities observed in this connection were the average energy per molecule and the total number of molecules in a given layer. Of course, this total energy fluctuated and drifted only in the canonical or constant NVT algorithm, but the numbers of timesteps listed in Table I indicate that these fluctuations and drifts were followed for long time intervals, especially for the highest coverages where the drifts were particularly noticeable and long-lived. Figure 2 shows dependence of the system energy upon timestep number in the NVT system with 240 N, molecules over a period of 15000 timesteps; the accompanying plot of the number of molecules held in the first layer indicates that the energy drift in the first 10000 timesteps is probably due to the slow increase in first-layer molecules. The second part of the plot indicates that the statistical averages obtained over the final 5550 timesteps indicated in Table I are reliable measures of the equilibrium statistical properties in this case. Similar analyses were made for the other systems considered, with results indicated in Table I. In all cases, molecular configurations in the NVT ensemble were utilized for statistical averaging. At low
-
222 Langmuir, Vol. 2, No. 2, 1986
-310
5000
l0000
Vernov and Steele
,
-310
5000
6
5000 Timesteps N1 IOOC
9%-0 Timesteps
2 r
2 r
Figure 2. Time dependence of the energy U* (in reduced units) and the number of molecules in layer 1 are indicated here for a simulation in the NVT ensemble at Nt = 240 molecules. Two time sequences are shown, with the second pair of figures referring to a later time than the first and, specifically, to a time interval during which the system appeared to be fluctuating around equilibrium. coverages, long-term drifts were less of a problem, but nevertheless equilibrations were run for the initial 3000 timesteps plus a variable number not less than 3800 timesteps before statistical averaging was begun. 3. Results: Singlet Densities and Energies As mentioned above, these simulations are used to extract considerable microscopic information not available by other techniques. As a first example, the molecular center of mass (com) density was evaluated as a function of position. Although the use of a gas-solid potential depending upon all three coordinates means that the density is a function of all three, we first present the surface-averaged density p ( z ) , which gives an indication of the distribution of molecules in layers in these systems. Figure 3 shows such curves, evaluated for all nine systems of Table I. It is obvious that these simulations include a range of coverages from nearly 1 monolayer up to one in which partial filling in both the second and third layer is present. (Such partial filling has been observed previously in the simulations of multilayers of argon on graphite at 90 and 115 K.24) It is also clear that division of these systems into first, second, third layers is unambiguous. The most probable first-layer gas-solid separation distance hardly varies with increasing total coverage and is equal to -1.45 (reduced units) or 0.357 nm-a number which is only slightly larger than uNCand, thus, which indicates that the molecules are predominantly lying parallel to the surface a t this temperature. (The small coverage-dependent changes in the density curves near their first-layer maxima are real and will be discussed in more detail below.) This point is reinforced by plots of the z dependence of the nitrogen atom densities shown in Figure 4,where sharp maxima are observed a t a distance corresponding to the minimum in the nitrogen site-solid potential function. The tailing off in these curves at slightly larger distances is due to configurations where one end of the molecule is at the bottom of the potential well, but thermal agitation has caused the other to move away to a larger distance. Note finally that the broad N-atom density peak indicates relatively little preferential orientation of the molecules in the second layer-however, the degree of orientational ordering can be seen more clearly in other calculations to be discussed below.
Figure 3. Orientationally and surface averaged density of molecules is shown here as a function of distance from the surface for all the systems simulated. Numbers are run numbers, as indicated in Table I.
7
L*
Figure 4. Dependence of nitrogen atomic densities upon molecule-solid distance are shown here for two simulations-these show small but significant differences from the center-of-mass densities of Figure 3. If one now defines layer separation planes at z given by the minima in p ( z ) , the fluctuating number of molecules in a layer is readily calculated from the simulation data. The time dependence of the number of the molecules in the first layer is shown in Figure 5 for all systems except run 1,which was completely and continuously in a single layer. Evidently, single molecule transfers between the first and second layers occur on a time scale of -500 timesteps or -1 ps, at least for coverages where significant second-layer formation has occurred. Since the system contains 100 first-layer molecules, one finds a lifetime of 100 ps for the first second layer transfer process. A quantity of considerable interest is the “monolayer coverage” of N2 on a solid surface. For the graphite basal plane, this quantity is readily calculated, either by integration of the first peaks in Figure 4 or by time averaging
- -
-
Langmuir, Vol. 2, No. 2, 1986 223
Computer Study of Multilayer Adsorption A
A Expt.
N I 98 97
:"i ,
N I 102 101
looo
5
,
IO ti mesteps/1000
o Total U*
ifF 100 "0
5
IO iimesteps/1000
I02 IOIO
-
g1
5
IO
timesteps/1000
Figure 5. Time-dependenceof the number of molecules in the first layer is shown for the various total coverages simulated-see Table I. Note that 1000 time steps 2.2 ps. (Molecules in run 1never left the first layer.) The arrow for run 5 denotes the time of the end of the simulation.
I.o
-250
2.0
N;^
Figure 7. Average values of the reduced potential energy per molecule of nitrogen adsorbed on graphite are plotted as a function of total coverage in reduced units. These results are compared with experimental points obtained from the integral heats of Also shown are Ums,i*,the molecule-solid contriads~rption.~ butions to the total energy due to molecules in layer i (calculated per molecule). Note the vertical scale change for the second-layer energies. Figure 6. Dependence of the number of molecules in the first layer upon the total number of molecules adsorbed. Reduced units are defined by dividing by the number in a commensurate layer. the plots of Figure 5. We define reduced coverages by evaluating Ni*= Ni/N,,,, where N,, is the number of molecules in a commensurate monolayer. Figure 6 shows N,*, the first layer coverage as a function of Nt*,the total coverage. Several interesting features emerge: First, the commensurate coverage is a preferred one. (Thermodyn a m i ~ ' and , ~ neutron scattering experiments4i5 both indicate that this is an expected result.) As the total molecules in the system (and, thus, in the second layer) is increased above Nt* = 1, more molecules are forced into the first layer. Ultimately, additional increases in total molecules cause this first layer density to decrease again to a value close to but slightly higher than commensurate. These density changes are not insignificant, amounting to -lo%, and it is safe to assume that they are temperature-dependent. An interesting question arises concerning the role of changes in molecule-surface orientation in connection with this phenomenon; we defer discussion of this point to follow a presentation of the average energies of the adsorbed N2 molecules. In classical statistical mechanics, energies of adsorption are due only to changes in the potential energies of the molecules involved. We can generate a detailed picture of this process by splitting the total potential energy of the adsorbed film into its component parts, which are Umsj*, the average N2-solid interaction for molecules in the j'th layer, and Uji*,the average N2-N2energy for pairs of N2 molecules, one of which is in layer i and one in layer j . (This "lateral interaction" could be further subdivided into electrostatic and nonelectrostatic parts, but we omit this detail here.) Figures 7 and 8 show average potential energies obtained from the simulations. In Figure 7, the total energy U* at coverages equal to or greater than commensurate is com-
I
I
I
I
I
Figure 8. Contributions to the total average potential energy due to N2-N2 interactions. Vi.*denotes the energy for the molecule in layer i interacting witk those in layer j and calculated per molecule in layer j . Points at coverages Nt* < 1were taken from the previous study.16 Values of U,,* were calculated only for Nt*corresponding to significant second-layer coverage. pared with the previous simulations at lower coverages and with values obtained from the experimental integral heats of ads~rption.~ The large discrepancy between simulation and experiment at low coverage is due almost entirely to surface heterogeneity in the experimental graphite substrate; when a reasonable correction is made for this, excellent agreement is obtained, as was shown in the previous work.16 Thus, the agreement between experiment and simulation encourages the belief that the interaction en-
224 Langmuir, Vol. 2, No. 2, 1986
Vernov and Steele
ergies assumed in this model may be close to those for the actual Nz-homogeneous graphite system. The general behavior of the curve of total energy per molecule vs. coverage lends support to the generally accepted descriptions of physisorption on a uniform surface, especially when the overall energy is supplemented by calculations of the component parts. The rise in (negative) energy across most of the monolayer regime is due primarily to increasing lateral interactions within the layer; near the commensurate density, a sharp decrease in the energy is observed which is primarily due to the initiation of second-layer formation with its smaller Um,,2*value. Nevertheless, it is interesting to see that this picture does not entirely account for the observations; in fact, there is a significant coverage variation in the average moleculesolid energy U,,,,* which, for example, amounts to a change of - 3 (reduced units) compared to total energy change of -6.5 in going from zero to commensurate coverage. This fact, which actually derives from the previous study,16must be taken into account when making estimates of lateral interaction parameters in monolayers of nonspherical molecules. It is caused primarily by changes in the molecule-surface orientation with first-layer density. Here, it can be seen that the changes in Ums,l*,although much reduced, continue as second-layer filling proceeds. Calculations of the density-induced changes in orientational distributions (which will be discussed below) show that the second layer adsorption does indeed affect the first-layer distributions and, thus, the first-layer molecule-surface energy. As expected, Ums,p*, the average gas-solid energy for molecules in the second layer, is nearly independent of total coverage and thus of second-layer coverage. It is also amounts to of Ums,l*,as is often found for physisorption on uniform surfaces. The N2-N2 energies are shown in Figure 8 and also exhibit the expected behavior for a high-temperature fluid. The linear increase in (negative) energy with increasing density in the layer is due to the lateral N2-N2 interaction in a van der Waals-like fluid; the leveling off in Ull* reflects primarily the leveling off in Nl* with increasing Nt*, as shown in Figure 6. The interlayer interaction U12* (calculated per molecule in layer 2) is constant for most of the filling range of layer 2, as expected for the interaction of a second-layer molecule with the dense first layer. It is likely that the change in this energy at high coverage is more apparent than real and is due to some third-layer occupancy which was not properly excluded from the calculation of second-layer coverage.
0
I
o
0
Figure 9. Order parameter OP,defined in eq 1 is plotted as a function of Nt*,the total coverage in reduced units.
mensurate lattice and T~ denotes the x,y coordinates of the com of molecule i (in the first layer). A problem with this order parameter is encountered in this work which was absent from the previous low-temperature calculations. This arises from the fact that although OP, = 1for molecules in a perfect commensurate layer on one of the three sublattices offered by the graphite basal plane, it is equal to if the commensurate phase forms on either of the other two. If, as is likely, the commensurate phase forms on more than one sublattice or if it changes sublattices as a function of time (due to random interlayer transfers of N2 molecules, for example), OP, may become quite small and will clearly no longer be a useful measure of the degree of commensurability. For this reason, a "local" order parameter was also evaluated. This calculation is based on a suggestion by Whitehouse et al.23and was utilized in the previous simulation of N,/graphite at 75 K.16 One defines F(x*,y*) =
2a(x* - y*/fi) cos 2a(x* y*/&)
COS
+
+ + cos ( 4 a y * / f i )
(2)
Ordering The observation of sharp features in the isosteric heats of adsorption and in the neutron diffraction of N2 on graphite at 73.6 K gives convincing evidence for the formation of solidlike commensurate monolayer phase under the conditions mimiced in these simulations. It is difficult to simulate the thermodynamic energy of adsorption in this region with sufficient accuracy to confirm the experiments, but one can calculate various translational order parameters which should shed light on this problem. The two possibilities investigated in this work are both based on in-plane com density distribution calculated relative to the periodicity of the surface. First, a structure factor denoted by OP, is defined as14
where x*,y* define the com position of the N2 molecule, with origin at the center of a hexagon (and with substrate C atoms at x*,y* = 0, &l/v'3; kl/*, f1/2v'3). Thus, for a given molecule, F ranges from -1.5 when the N2 is directly over a C atom to +3.0 when it is directly over the center of the hexagon (defined to be an adsorption site). Evidently, a perfect commensurate layer will have F = 3 for all molecules regardless of the sublattice occupied. A useful measure of "localization" in this sense will be to evaluate the probability that a molecule has a given value of F , since sitewise adsorption produces a sharp spike in this function a t F = 3 whereas completely mobile adsorption would yield a constant probability. Note that this behavior would be observed regardless of which sublattice is occupied by adsorbed molecules or, for that matter, regardless of whether the order is long ranged or not. Simulated values of OP1 are shown in Figure 9 for various total surface coverages. It is evident that a significant degree of commensurate behavior is present at coverages very close to those for a commensurate layer. However, at higher coverages where interlayer transfers do occur (see Figure 3) the monolayer appears not to be commensurate. It is interesting to note the slight rise in OP, at the highest coverages, where the data indicate partial return to the commensurate value for first-layer coverage. Nevertheless, it was felt that calculations of site distributions might throw additional light on this problem and, indeed, the curves in Figure 10 show that this is the case. The function p ( F ) plotted there is defined asz3
OP, = ( I / N ) ( ~ ( e L k + l ~eLky71)) T~
dF) = n ( F ) / a ( F )
4. Results: Translation a n d Orientational
(1)
I
where k,, k, are the reciprocal lattice vectors of a com-
(3)
where n(F) is the fraction of the total molecules with a value of F in the range F - AF/2 to F + AF/2 and a(F)
"
Langmuir, Vol. 2, No. 2, 1986 225
Computer Study of Multilayer Adsorption
c
c
c
piF1 2 0.50
00 .-
PiFi
;p:k, , -, ?I50
000
150
300
!I50
000
150
-i.w-nso
:k, 300
9150 0 0 0
150
0.00 n5o 1.00
c
8
2
0.50 l
300
0 0 . 0
0.00-i.m-aso
t
1.00
o
o
0
h
n
-I.W-Q5o 0.00 0.50 100
I50
3iO
00-9
is the corresponding fraction of the area of a basal plane hexagon (which is equivalent to the fraction of molecules in the same interval of F for a uniform distribution of molecules). (Values of a(F) for A F = 0.225 have been calculated previously by the Monte Carlo technique for this system.)16 One sees significant localization in the first layer over the adsorption sites in systems 1 and 2 (see Table I), almost completely random probability distributions for the coverages of systems 3-7, and a return to partially localized behavior for 8 and 9. (Both p ( F ) and OP, indicate this, but p ( F ) seems the more reliable measure.) At sufficiently high second-layer coverage, a few molecules are extracted from the first layer and the tendency toward commensurate behavior reappears. We conclude that the degree of localization (into a commensurate layer) as measured by p ( F ) is extremely sensitive to the layer density. For example, Figure 5 indicates that the number of molecules in layer 1 fluctuates between 96 and 97 in run 2, and between 96 and 98 in run 3. However, the curves of p ( F ) are quite different in the two cases, and indicate that this slight change in N1 is sufficient to push the molecules away from the site centers and, thus, to destroy commensurability. We have also evaluated p ( F ) for molecules in the second layers of the various simulations and find, as expected, no tendency for preferential adsorption over points on the graphitic sites. An important aspect of this work is of course the calculation of molecular orientations and their dependence upon coverage. It is convenient to split this problem into in-plane and out-of-plane parts, particularly since the inplane ordering of these molecules was found to be negligible for all coverages at 73.6 K. Thus, we consider only the statistical properties of cos /3, where j3 is the angle between the N2 molecular axis and the normal to the graphite surface. Figure 11 shows the calculated probability distributions for the first-layer molecules in all the simulation systems, and Figure 12 shows the same functions for molecules in the second layers. The central peaks in Figure 11 show a preference for coplanar molecular configurations for all first-layer densities. However, significant variations in these curves with changing density are also seen. In particular, the numbers of molecules that are almost perpendicular to the surface (cos j3 = f l ) is nonnegligible and increases as the layer becomes more dense. When the first layer density is significantly larger than commensurate, some of the molecules are squeezed into vertical orientations. One consequence of this is the change in Ums,l*shown in Figure 7; it is also likely that
0
o
o
-100-0.50 0.00
am
1.00
h
0
u
o.00-100-0.50 0.00 0.50 1.00
-1.W-0.50 0.00 0.50 1.00
c
5
0
v
c
l"o,h/c
0.50
W - l?
F
Figure 10. Probability densities for site localization p ( F ) defined in eq 3 are plotted as a function of the position variable F for the first-layer molecules in all systems simulated in this work. (See ref 16 for similar curves evaluated at lower N2 coverages.)
5
00 .-
100
5
.
c
9/50 O b 0
0.00 a50 1.00
n
o o k -i.w-nm 0.w 0.50 1.00 cos
0
p
.
0
o
k
- 1 w - 0 ~ 0ooo 050 100 cos p
O
O
O
aw
-IW-050
cos
a50
p
L
ioo
Figure 11. Probability distributions for the molecule axissurface plane orientation @ are shown here for the first-layer molecules in all simulations here. (Similar curves are given in ref 16 for lower coverages). The probability of a given value of cos p is shown. For homonuclear molecules, these curves should be symmetric about cos @ = 0; any differences shown are due to statistical uncertainty. -
I O O ~
9
. 000-
100
0001 ' -100
I
1
'
000
I
1
I
1
IO0
1.00
1
1
,
0.00 -1.00
6
1.00
1.00 0
000
2
; ; 100
6
-100
7
5
0
0.00 COlp
b
,
050+,
1.00
o~o-4.00
000
1.00
cos B
Figure 12. Same curves as in Figure 11, but for the molecules in the second layers of the simulation systems.
Ull*, the average N2-N2 lateral interaction, will be affected but the strong first-order density dependence of this energy a t low coverage makes it impossible to use it to estimate the effects of small orientational changes upon the overall interaction. It should also be noted that the slight decrease in first-layer density a t the highest total coverages is accompanied by a small but noticeable decrease in the fraction of molecules that are nearly perpendicular to the surface. The curves of Figure 12 for the out-of-plane orientational probabilities of second-layer molecules indicates that some preference for coplanar orientations persists, especially for the lowest second-layer densities. Of course, the small number of such molecules means relatively poor statistics under these conditions. Nevertheless, the loss of this preference at higher density is obvious and not easy to understand, especially since low-temperature (35 K) simulations of N2 bilayers show strong coplanar ordering in both layers. 5. Results: Pair Correlations Functions The physical picture derived from these simulations of N2at high coverage on graphite at 73.6 K is one of a dense,
226 Langmuir, Vol. 2, No. 2, 1986
Vernov and Steele 5
I
2
-i tl
I I
Idllrlllllll
l , i l , l l l l l l l
'0
2
4 ,*6
8
10
'0
2
4 r+6
8
IO
Figure 14. Same as Figure 13, but for molecules in the second layer. Only the peaks of the atom-atom correlations are shown (by solid lines). r*
r*
r*
Figure 13. Two-dimensionalprojections of the center-of-mass and atom-atom pair correlation functions are shown for the first-layer molecules in the three systems indicated (see Table I for coverages). The peaks in the com correlations occur at the distances expected for a commensurate lattice (shown by the vertical lines) with heights that are more or less in agreement with the expected values, after allowing for the considerable thermal broadening present in these systems. The center panel shows the correlation functions for a layer which exhibits separation distances slightly smaller than for a commensurate lattice. The atom-atom correlations are shown by dotted lines and are lower and broader than the com peaks, due primarily to the averaging effect of nearly free in-plane rotation of the atoms about the com position. (Similar curves for lower coverages are given in ref 16.) solidlike monolayer in which the molecules minimize their gas-solid potential energies by being predominantly parallel tQ the surface. Changes of the total coverage produce slight changes in density of this underlayer which are associated with partial reorientations of the N2 molecules. (To make room for additional first-layer molecules, for example.) Another probe of the structure of these N2 layers is offered by the pair correlation functions.32 Thus, both the intralayer com and the nitrogen site-site pair correlation functions were evaluated from the simulations. Intralayer here signifies the two-dimensional correlations projected into a plane parallel to the surface for the molecules in either the first or the second layer. Correlation functions for first-layer molecules are shown in Figure 13 for selected densities. Both the com and the site-site correlations demonstrate the damped long-range order expected for a solidlike structure at high temperature (and thus, with significant thermal motion). It is of interest to compare the relative heights and positions of the peaks with the expected values for the 0 K commensurate lattice. These are indicated by vertical lines of height proportional to the number of neighborsf radial distance. After taking account of thermal broadening, the results for N,* = 1.0 are in satisfactory agreement with those for the commensurate lattice and, interestingly, indicate that each of the peaks for second-, third-, and fourth-neighbor shells actually should be associated with pairs of overlapping, thermally broadened peaks obtained for a regular, triangular lattice. At the intermediate coverage N,* = 1.35, the basic structure is unchanged, but the peak positions are shifted slightly inward, indicating a solidlike, regular triangular lattice which is not commensurate, in agreement with conclusions drawn previously from the simulated one-particle properties. At the highest density, the peaks shift outward again and are thus closer to the commensurate lattice spacings than those for N,* = 1.35, which is also in agreement with previous conclusions. Turning now to the site-site correlations, these curves are completely consistent with the picture of an ordered (32) Nicholson, D.; Parsonage, N. G. 'Computer Simulation and the Statistical Mechanics of Adsorption"; Academic Press: London, 1982.
com lattice and free in-plane rotation of the nitrogen molecule sites about the com positions. Since the N2 site-com separation is 0.22 (reduced units), one expects the peaks to be centered on the com distances but smeared out by the rotational averaging, just as observed. The character of the second-layer pair correlations shown in Figure 14 is quite different from those for the first, at least for the temperature and coverages considered here. Not only is there more (three-dimensional) orientational smearing in the site-site functions, there is also clearly fluidlike behavior for the com distributions, which show strongly damped peaks for only a few shells of neighbors. 6. Conclusions Taken together with the previous simulation,I6this study of nitrogen adsorbed on graphite a t 73.6 K in the monolayer-bilayer regime provides us with a detailed picture of the molecular behavior of this system. While the level of agreement between the simulation results and experiment is gratifying, an interesting aspect of this work is the information generated concerning molecular orientation relative to the surface and its role in determining the heats of adsorption (and, ultimately, the adsorption isotherm-to be discussed elsewhere). The magnitude of the periodicity of the molecule-solid potential is not large compared to kT-a molecule adsorbed over the center of a hexagon has energy -0.43 reduced units or -14 K lower than if it were randomly distributed (or if the periodic part of the potential was absent) at a distance corresponding to the minimum in the z dependence of the potential, and the calculation shows that this value drops off very steeply as the molecular com moves outward from the potential minimum. It appears that the effect of this periodic term is reinforced by the inherent tendency of these molecules to pack with a spacing close to that for a commensurate lattice. Thus, the observed solidlike behavior of the first layer of nitrogen at high coverage must be a result of the combination of the substrate periodicity and a fortunate matching with the molecular size parameters. To illustrate this, the contribution of the average of the periodic energy Uperto the total can be estimated with the aid of the p ( F ) curves. Assuming that the z-dependent average can be separated from the average over area, one has
uper = -2E,(z)
s
hexagon
Fp(F) dA
(4)
If we convert the integral to a finite sum over LW (remember that dAfdF = a ( F ) ) ,the calculation is straightforward. The results are shown in Figure 15, where we have used p ( F ) from the previous as well as this work to generate estimatesof this energy over the entire coverage region. If we take E1(z)to be the value for z corresponding
Langmuir, Vol. 2, No. 2, 1986 227
Computer Study of Multilayer Adsorption 0.8
I
I
I
I
I
I
T
R
0.61
(and, possibly, for 20-30" below this), the facts that the distribution of orientation relative to the surface plane are not random and depend upon both coverage and temperature means that an accurate description of the thermodynamics of such systems will necessitate taking accountof the orientational changes in average moleculesurface and in molecule-molecule energy. Evidently, such complications are absent in the Kr-graphite system. An additional feature of the orientational distributions obtained for the N2 graphite system is the coupling between orientation and maximum first-layer coverage. Changes in this first-layer coverage are associated with changes in the orientational distributions, which are in turn associated with changes in the density of molecules in the second layer. These concepts have been suggested in various previous papers, but the present simulation study provides a detailed picture not available heretofore. Finally, it should be emphasized that the simulations reported here should not be construed as giving a theory for the monolayer-bilayer adsorption of Nzon graphite. However, they have generated a great deal of information which should serve as a guide and a test of future theoretical arguments. As an example, the phase diagram shown in Figure 1would appear to need some adjustments in the region of coverage higher than commensurate. Certainly the sequence fluid CD TID at 74 K seems to be followed in the simulations. However, the model system does not appear to show sharp phase transitions. In addition, the simulated CD phase is not unambiguously solid. No evidence for first-order phase changes into or out of the CD phase was obtained, and, furthermore, the time-dependent properties (to be published elsewhere) indicate considerable molecular mobility parallel to the surface in all systems simulated here. More importantly, the phase diagram of Figure 1 should indicate the formation of a mixed TID-bilayer adsorbed phase at total coverages higher than commensurate. Finally, there seems to be no way to adapt such a phase diagram to show the changes in first-layer density in the mixed TID-bilayer system that occur with increasing total coverage.
0'4[Jr I
I1
I
0
0.2
I 0
I
0.4
I
0.8
I
I
1.2
1.6
I
2.0
I
2.0
I 2.8
N:
F i g u r e 15. Averages of the periodic function F for first-layer molecules, evaluated using the curves of Figure 10 plus the results given in ref 16 (filled points) are shown as a function of total coverage. These values can be used to estimate the contribution of the periodic part of the gas-solid energy to the total potential energy of a molecule adsorbed in the first layer and exhibit maxima a t the points where the layer is closest to commensurate packing.
to the minimum in the potential, we should obtain a reasonably realistic value for Upr. Since this amounts to only 0.27 reduced units just a t the commensurate coverage where it is a maximum, one concludes that the periodic energy is a nonnegligible but small part of the total and that changes in this energy as the layer becomes noncommensurate are also a small but nonnegligible part of the total energy change. The shift toward perpendicular molecular orientation as the density increases is not unexpected. The magnitudes of the energy changes accompanying these shifts are a significant part of the overall change and show clearly that comparison between experimental heats of adsorption and theoretical estimates of the lateral interaction energy for nonspherical molecules in general and for N2 in particular should take account of this reorientation effect. Perhaps the most surprising finding of this work is the observation of a small decrease in first-layer density as the second-layer density is increased past completion. The reasons for this are not entirely clear, but one might guess that as the second-layer energy becomes more negative a t moderately large density, the energy separation between first and second layers becomes less, and any resulting molecular redistribution would be in the direction observed. The argument that N2/graphite and Kr/graphite have similar phase diagrams at high temperatures where the N2 rotation is free or nearly free is now seen to require some caution in its applications. Although the in-plane orientational averaging is certainly complete for Nz a t 70 K
- -
Acknowledgment. Support from a grant from the Materials Research Division of the National Science Foundation is gratefully acknowledged. This work was done while A.V.V. was an exchange visitor at Penn State under the auspices of the International Research and Exchange Board and the Ministry of Higher Education of the U.S.S.R. Helpful discussions with Moses Chan and Dominic Tildesley are also acknowledged. Registry No. N B ,7727-37-9; graphite, 1182-42-5.
