On the phase behavior of mixed Ar-Xe submonolayer films on graphite

Using Monte Carlo simulation methods in the canonical and grand canonical ensembles, we discuss the melting and the formation of ordered structures of mixed Ar-Xe submonolayer films on graphite. The calculations have been performed using two- as well as three-dimensional models of the systems studied. It is demonstrated that out-of plane motion does not affect the properties of the adsorbed films as long as the total density is not close to the monolayer completion. On the other hand, close to the monolayer completion, the promotion of particles to the second layer considerably affects the properties of mixed films. It has been shown that the mixture exhibits complete mixing in the liquid phase and freezes into solid phases of the structure depending upon the film composition. For submonolayer densities, the melting temperature exhibits non-monotonous changes with the film composition. In particular, the melting temperature initially increases when the xenon concentration increases up to about 20%, then it decreases and reaches minimum for the xenon concentration of about 40%. For still higher xenon concentrations, the melting point gradually increases to the temperature corresponding to pure xenon film. It has been also demonstrated that the topology of phase diagrams of mixed films is sensitive to the composition of adsorbed layers.

A vast majority of theoretical studies of mixed adsorbed layers has been based on lattice gas models [28][29][30][31][32][33], which do not constitute a good basis for the discussion of the incommensurate-commensurate transitions in the films of rare gases on graphite. Such models cannot properly describe the incommensurate floating solid. However, there have also been some attempts to construct theoretical models for the commensurate-incommensurate transitions in monolayer films of rare gas mixtures on graphite [34][35][36]. The primary aim of the models proposed by Marti et al. [34,36], was to explain the anomaly in the phase behavior of Ar-Xe and Kr-Xe submonolayer films on graphite. Namely, experiments have demonstrated [34] that less krypton than argon is needed to induce the formation of commensurate phase. On the other hand, a rather general mean-field model of Villain and Moreira [35] requires the introduction of several approximations in order to make the resulting equations numerically tractable. Consequently, the agreement with experimental data is rather poor. Nevertheless, these authors have derived qualitative phase diagrams for the Ar-Xe mixture adsorbed on graphite. The theory predicts the existence of incommensurate, Ar-like and Xe-like structures at low temperatures, and the formation of the commensurate structure at higher temperatures over a rather limited range of the Xe mole fraction.
In this work, we present and discuss the results of rather extensive Monte Carlo simulations of mixed Ar-Xe films on graphite. Our main goal has been to investigate the structure of low temperature solid phases as well as to determine the changes of the melting temperature with the mixture composition. However, we also discuss the evolution of phase diagrams resulting from the changes in the film composition.
The paper is organized as follows. In the next section we present the model used and describe the Monte Carlo method used to determine the properties of mixed submonolayer films. Then, in section 3 we briefly discuss the behavior of pure Ar and Xe films. The last section 4 is devoted to the presentation of the results for the mixed Ar-Xe films.

The model and Monte Carlo methods
The interaction between adsorbate atoms is assumed to be represented by the (12,6) Lennard-Jones potential u i,j (r i j ) = 4ε i,j σ i,j /r i j 12 − σ i,j /r i j 6 , (2.1) where r i j is the distance between a pair of atoms and i and j mark the species Ar and Xe. The values of the parameters ε i,i and σ i,i used in this work are given in table 1. The corresponding parameters The potential (2.1) has been cut at the distance 3σ i,j .
We are aware of some drawbacks that the assumption of the LJ potential has got, and that other authors have used different potentials to reproduce experimental data for adsorbed films of pure rare gases on graphite [37]. Also, we have not taken into account the surface mediated interactions, which are known to affect the strength of adsorbate-adsorbate interaction in the vicinity of solid substrates [7].
The interaction of rare gas atoms with the graphite basal plane can be represented by the potential proposed by Steele [38] v In the above, the first term in the square brackets is the fluid-solid potential averaged over the entire surface, while the second term represents the corrugation part of the fluid-solid potential. Assuming that the interaction between an adsorbate atom and the carbon atom of the graphite substrate is also represented by the Lennard-Jones potential, the explicit expressions for v 0,i (z), the Fourier components v k,i (z) and the functions f k (x, y) are given by the following equations: with the sum running over all graphite reciprocal lattice vectors of the length q k . In the above equations A i = σ i,C /a 1 , where a 1 = 2.46 Å is the graphite lattice constant, i.e., the distance between the centers of adjacent carbon hexagons, the values of σ i,C and ε gs,i (i = Ar, Xe) are given in table 2, ∆z = 3.4 Å is the spacing between graphite planes, a s = 5.24 Å 2 is the area of the graphite unit cell, K 2 and K 5 are the modified Bessel functions of the second kind and of the second and fifth order respectively, and q k 's are the lengths of the graphite basal plane reciprocal lattice vectors. In the case of only partially filled monolayer films and at sufficiently low temperatures, the promotion of the second layer is likely to be negligibly small. This allows us to consider a simple strictly twodimensional model with the external field of the form where the parameter V b,i (i = Ar, Kr) determines the amplitude of the corrugation potential. The magnitudes of V b,i (V b,Ar = 0.07 and V b,Xe = 0.08) have been adjusted in such a way that the results for each component are more or less consistent with the full 3D calculations.
Simulations have been performed using the Monte Carlo method in the canonical and grand canonical ensembles [39,40]. In the case of two-dimensional model, the rectangular simulation cell of the size La 1 × L 3a 1 /2, with L = 60 and with the standard periodic boundary conditions has been used. In three dimensional calculations, the simulation cell has been the rectangular parallelepiped of the size 60a 1 × 60a 1 3/2×10a 1 , with the periodic boundary conditions applied in the directions parallel to the substrate surface and with the reflecting hard wall located at z = 10a 1 .

23601-3
The quantities recorded included the average potential energy, 〈e〉, the contributions to the potential energy due to the fluid-fluid interaction, 〈e gg 〉 and the contributions due to the fluid-solid interaction for each component 〈e gs,i 〉 and the heat capacity obtained from the fluctuations of the potential energy, (2.8) In order to monitor the structure of solid phases we have used radial distribution functions, g i j (r ), for different pairs of species i and j , and appropriate order parameters. The formation of hexagonally ordered phases has been monitored using the bond-orientational order parameters [41,42] measured separately for each adsorbate (i = Ar or Xe). In the above, the first sum runs over all atoms of the i -th component, the second sum runs over all nearest neighbors of the same type, φ m,n is the angle between the bond joining the atoms m and n and an arbitrary reference axis, chosen here to be the xaxis of the simulation cell, and N b,i is the number of bonds between pairs of the like atoms. Also, we have calculated the total bond-orientational order parameter where the first sum runs over all atoms in the system and the second over all nearest neighbors. The above defined bond-orientational order parameters make it possible to detect the hexagonally ordered structures, but are not suitable to distinguish between the commensurate and incommensurate phases. In the commensurate phase, the atoms are localized over the centers of carbon hexagons, and the appropriate order parameter allowing to monitor such localized structures can be defined as [43] exp iq n r m,i . (2.11) The first sum is taken over all atoms of the i -th component, while the second sum runs over the six reciprocal lattice vectors q n of the graphite substrate and r m,i is the position of the m-th atom of component i .
The above defined order parameters have been supplemented by the corresponding susceptibilities where 'op' stands for any of the above given order parameters. When grand canonical simulations have been carried out, we have also recorded the adsorptiondesorption isotherms.
Throughout this paper, we use reduced quantities, assuming that the graphite lattice constant a 1 is the unit of length, and the Lennard-Jones parameter ε Ar,Ar is the unit of energy. We have decided, however, to give the temperature in Kelvins, as it allows for easier comparison of our simulation results with experimental data. All the densities are expressed in commensurate monolayers.

The results for pure Ar and Xe films
The phase diagram of xenon monolayer on graphite exhibits the vapor-liquid critical point and the triple point (see parts (a) and (b) of figure 1). The triple point temperature, equal to T tr ≈ 89 K, agrees quite well with some experimental data [44], yet it is lower than the value of about 99 K stemming from other experiments [22,45] and from the recent Monte Carlo results of Przydrozny and Kuchta [37]. The underestimation of the triple point temperature is associated with our choice of the interaction potential and its parameters. Przydrozny and Kuchta applied a semi-empirical potential proposed by Aziz and Slaman [46], while we have used a simple Lennard-Jones potential. We should mention that in the earlier molecular dynamics simulation studies of the melting transition of xenon on graphite [47,48], also based on the Lennard-Jones potential, but with slightly different values of the parameters ε Xe,Xe and σ Xe , the triple point temperature was found to be located at T tr ≈ 0.4ε Xe,Xe /k, while our result is T tr ≈ 0.402ε Xe,Xe /k. Also, the critical point temperature T cr ≈ 109 K is lower than the experimental value by about 18 K [45]. Nevertheless, the qualitative agreement with the available experimental data is good enough to assume that the results for the mixed films are also qualitatively correct.
In the case of argon, the phase diagram derived from our grand canonical simulation (see parts (c) and (d) of figure 1) agrees very well with experiment. In particular, the triple point temperature, equal to 49.5±0.5 K, is practically the same as the experimental value of 49.7 K [11]. Also, the critical temperature agrees very well with experiment [49].
The freezing of submonolayer xenon and argon films leads to the formation of incommensurate structures. The xenon incommensurate solid attains the density of about 0.85 at the monolayer completion and is expended with respect to the commensurate ( 3 × 3)R30 structure. Figure 2 shows the temperature

23601-5
changes of the order parameters Ψ 6,Xe and Φ Xe for xenon films of different total densities and one sees that the order parameter Φ Xe is quite low even at very low temperatures, indicating the lack of localization of adatoms over the minima of the graphite lattice. On the other hand, the bond-orientational order parameter Ψ 6 demonstrates the formation of hexagonally ordered phase below the freezing point. The argon incommensurate solid also shows a well developed hexagonal symmetry, but it is contracted with respect to the commensurate ( 3 × 3)R30 structure and attains the density of about 1.25.
At sufficiently low temperatures, the film exhibits epitaxial rotation of about 3 degrees (see figure 3). Again, this result agrees very well with earlier theoretical [18,19] and computer simulation [16,27] results.
The results of Monte Carlo simulation for pure argon and xenon films given above will serve as reference data for the study of mixed films.

The results for mixed films
We begin with the presentation of canonical ensemble Monte Carlo simulation results aiming at the determination of the melting temperature and the structure of solid phases in submonolayer mixed films. Since the solid phases (commensurate and incommensurate) exhibit hexagonal symmetry, the location of the melting point can be estimated using the bond-orientational order parameter, Ψ 6 , and its susceptibility, χ Ψ 6 . Of course, one also expects that the melting transition is manifested by sudden changes of the potential energy and the heat capacity anomalies. Figure 4 gives an example of our results, obtained for the submonolayer film of the total density ρ c = 0.4 and the xenon mole fraction equal to x Xe = 0.1667. Part (a) of figure 4 shows the heat capacity curve and one sees a sharp peak at the melting point at T ≈ 54 K. At the same temperature, the total potential energy u and the contributions to the potential energy due to Ar-graphite and Xe-graphite interactions exhibit sudden drops (see part (b) of figure 4). Finally, part (c) of figure 4, which shows the temperature changes of the bond-orientational order parameter Ψ 6 , and its susceptibility χ Ψ 6 demonstrates that the melting transition is accompanied by the loss of hexagonal ordering. It should be emphasized that a large increase of the Ar-graphite and Xe-graphite interaction energies accompanying the freezing transition marks a sudden increase of localization of the adsorbed argon and xenon in the solid phase. Upon a decrease of temperature, the localization of xenon gradually increases, while argon exhibits a decrease of localization at temperatures below T ≈ 30 K. This behavior can be attributed to the transition between the commensurate phase, stable at T > 30 K, and the incommensurate phase, stable at T < 30 K. Note that the transition is not accompanied by any changes in the

23601-6
Ar-Xe mixed film on graphite behavior of the bond-orientational order parameter, but produces a well seen heat capacity anomaly.
The inspection snapshots of configurations recorded during the simulation runs have shown that the commensurate phase is mixed, while the incommensurate phase consists of argon only. In the case of small xenon mole fraction, as in the system considered now, we expect to observe only a partially developed commensurate phase. Indeed, the snapshot given in figure 5 (a) shows that the film peripheries are predominantly occupied by argon atoms, which also show a rather high degree of incommensuration. At the temperature below commensurate-incommensurate transition, we find coexisting domains of mixed commensurate and argon-like incommensurate phases (see figure 5 (b)). In the snapshots given in figure 5, we have assigned the atoms to commensurate and incommensurate positions using the following order parameter [50]: φ(r) = cos(q 1 r) + cos(q 2 r) + cos([q 1 − q 2 ]r) , (4.1) and assuming that the atom is commensurate (incommensurate) when φ > 0 (φ 0).One should note that even in a rather small system used, consisting of only 480 atoms, the argon-like incommensurate domain exhibits epitaxial rotation, just the same as observed for pure argon films.
In order to determine the locations of the commensurate-incommensurate transition, we have monitored the behavior of the order parameters Φ Ar and Φ Xe , defined by the equation (2.11). Figure 6 shows the changes of these two order parameters with the xenon mole fraction at two different temperatures in  submonolayer films of the total density ρ c = 0.4. Quite similar results have been obtained for the films of different total densities and using two-as well as three-dimensional models. From the observed changes of the order parameters Φ Ar and Φ Xe , it follows that the increase of the xenon concentration leads to a sequence of changes in the film structure. For small x Xe we find that xenon is highly localized, while the degree of localization of argon increases with x Xe . In this region, the film consists of two coexisting phases: one being the incommensurate argon-like solid and the second being the mixed commensurate solid. A gradual increase of the xenon mole fraction causes the commensurate domain to become larger and the size of incommensurate domain to shrink gradually. Then, there is a region of xenon concentration over which both adsorbates are highly localized. This corresponds to the presence of pure mixed commensurate phase and terminates at the xenon mole fraction close to about 0.4. Then, both order parameters gradually decrease when x Xe increases up to about 0.75. In this region, the mixed commensurate phase coexists with the demixed xenon-like incommensurate phase. Finally, for x Xe exceeding about 0.75, the film consists of xenon-like incommensurate phase with the argon atoms located at its peripheries. This has been confirmed by the inspection of snapshots and radial distribution functions.
The central result of the canonical ensemble Monte Carlo study is given in figure 7, which contains the phase diagrams showing the locations of the melting transition and the regions of stability of different solid phases in the films of different total densities. Part (a) of figure 7 gives the results for submonolayer films of different total densities, equal to 0.4, 0.667 and 0.8. The results for ρ c = 0.4 and 0.667 have been obtained using a two-dimensional model, while those for ρ c = 0.8 have been obtained within a more realistic three-dimensional model. The locations of the melting point are more or less the same over a wide range of x Xe between 0 and about 0.8. The independence of the melting temperature of the total density indicates that the melting occurs at the triple point temperature. For the xenon concentration higher than 0.8 the triple point melting occurs for ρ c = 0.4 and 0.667, but not for ρ c = 0.8. This suggests that the density ρ c = 0.8 is higher than the liquid density at the triple point. A rather sharp increase of the melting temperature with the xenon mole fraction, for x Xe above 0.8, results from the fact that the increase of xenon concentration brings the film closer to the monolayer completion. We should emphasize that even for x Xe close to unity there is no trace of the promotion of adsorbed argon and xenon to the second layer, even at the temperatures above the melting point. Thus, the films remain practically two-dimensional.
Although we have not performed any simulation at ρ c = 0.8 using the two-dimensional model, it can be anticipated that the results should be quite the same as those obtained with the three-dimensional model. The two-dimensional approximation is expected to fail when the film density starts to exceed the monolayer capacity. We have carried out the canonical ensemble calculations assuming that ρ c = 1.0. This density is lower than monolayer capacity of pure argon film, but it is well above the monolayer capacity of pure xenon film. Part (b) of figure 7 shows the phase diagram obtained. In this case, the two-dimensional model works well only in the region of the xenon mole fraction not greater than about 0.2, and starts to overestimate the stability of the solid phase for higher xenon concentrations. Three-dimensional calculations have shown that argon is partially promoted to the second layer when the temperature becomes high enough. The temperature at which the second layer promotion begins depends upon the film composition and it is higher than the melting temperature only for x Xe not exceeding about 0.2 and lies below the melting temperature for higher xenon mole fractions. It is therefore not surprising that the two-dimensional One sees that the xenon mole fraction range over which the commensurate phase is stable is considerably wider than in the previously discussed films of lower total density (cf. figure 7 (a)). Moreover, we find that over a certain range of xenon concentrations, between about 0.48 and 0.83, the mixed incommensurate phase appears at the temperatures just below the freezing transition, whereas no trace of such a phase has been found in submonolayer films. Upon the decrease of temperature, this phase transforms either into the a commensurate phase, when x Xe is lower than 0.6, or into the coexisting commensurate and xenon-like incommensurate phase, when x Xe is higher than 0.6. Figure 8 shows the changes of the heat capacity (part (a)) and the order parameters (part (b)) in the case of the film with x Xe = 0.7. The heat capacity has two pronounced anomalies. The first one, at T ≈ 108.5 K, is the signature of freezing transition, accompanied by the development of hexagonal order in the film (see the behavior of the bond-orientational order parameter Ψ 6 in figure 8 (b)). At the temperatures between the freezing point and the second heat capacity anomaly at T ≈ 70 K, the solid phase is incommensurate and the order parameters Φ Ar and Φ Xe remain very small, as expected for the incommensurate phase. The inspection 23601-9 of snapshots has shown that the incommensurate phase is mixed and the argon is partially promoted to the second layer. The second heat capacity anomaly, at T ≈ 70 K, is due to the onset of the transition accompanied by a rather large increase of the order parameters Φ Ar and Φ Xe , indicating the increase of localization of Ar and Xe upon the lowering of temperature. Figure 9 shows the snapshots recorded at 72 K and 12 K, which demonstrate that the transition observed leads to the formation of domains consisting of the mixed commensurate and demixed Xe-like incommensurate phases. The snapshot recorded at 12 K also demonstrates that argon is partially promoted to the second layer, and forms a compact island of a solid-like phase. It is also noteworthy that the solid-like patch of argon in the second layer is located over the demixed xenon domain rather than over the domain formed by the mixed commensurate phase. This can be readily understood by taking into account the magnitudes of Ar-Ar and Ar-Xe interaction energies, measured by the Lennard-Jones potential parameters ε Ar,Ar and ε Ar,Xe , and of course ε Ar,Xe is considerably larger than ε Ar,Ar (cf. table 1). When the argon atoms from the second layer are located over the pure xenon patch, each of them has three xenon atoms from the first layer as nearest neighbors. On the other hand, if the argon patch were located over the mixed commensurate patch then some nearest neighbors from the first layer would be argon atoms, and this situation is energetically less favorable.
When the xenon concentration exceeds about 0.83, the first layer consists only of xenon, while all argon atoms are promoted to the second layer. Of course, when the amount of xenon in the film exceeds the monolayer capacity of pure xenon film, then the excess of xenon is also located in the second layer.  We now proceed to the discussion of the changes in the phase behavior resulting from grand canonical Monte Carlo simulation. The calculations have been performed under the condition of the fixed chemical potential of xenon, so that the xenon concentration in the film was not conserved. Along the adsorption isotherms obtained by changing the chemical potential of argon, the amounts of xenon change as well. Figure 10 shows two examples of adsorption isotherms, both recorded at T = 60 K, but with different values of the xenon chemical potential. It is quite evident that the gas-liquid transition is accompanied by a sudden increase of the xenon density and that this effect is much stronger when the chemical potential of xenon is higher. We have determined the phase diagrams for a series of systems with different values of the xenon chemical potential, µ Xe = −20.0, −19.5, −19.0 and −18.6.
In the case of µ Xe = −20.0, the amounts of xenon in the film are very low, with x Xe < 0.01, over the entire range of temperatures and film densities studied. It is, therefore, not surprising that the phase diagram obtained is very similar to that of pure argon film (see figure 11). In particular, the melting transition appears to be continuous and the solid phase is an incommensurate argon-like phase. However, we find that even very small amounts of xenon shift the locations of the triple and critical points towards higher temperatures. An increase of µ Xe to −19.5 leads to some qualitative changes in the phase diagram topology (see parts (a) and (b) of figure 12). In particular, at the temperatures below about 50.5 K, the two-dimensional gas condenses into the commensurate krypton-like structure, which undergoes a transition into the incommensurate phase when the argon chemical potential becomes high enough. This commensurate- The vertical dotted lines mark the locations of subsequent maxima in a perfectly ordered commensurate phase. Part (a) also shows the argon-xenon radial distribution function obtained for µ Ar = −9.9.

23601-12
incommensurate transition is continuous. One should note that the xenon concentration in the commensurate phase is rather low (x Xe < 0.2), and becomes still lower in the incommensurate phase. For still higher value of µ Xe equal to −19.0 the phase diagram topology remains the same and only the stability region of the commensurate phase becomes wider and extents up to T ≈ 57 K. Also, the commensurateincommensurate transition occurs at higher values of the argon chemical potential. At the temperature between the upper limit of the commensurate phase, i.e., T ≈ 57 K, and the critical point, the gas phase condenses into the liquid phase. The picture presented above is very well confirmed by the behavior of radial distribution functions. Part (a) of figure 13 gives the argon-argon distribution functions recorded at the temperature of 51 K and for the argon chemical potentials below and above the commensurateincommensurate transition. It is evident that the maxima, apart from the first one, coincide very nicely with the locations of subsequent neighbors in the commensurate phase. The stability of commensurate phase is due to the presence of argon-xenon nearest neighbors, and the argon-xenon distribution function (also shown in figure 13 as a dashed line) exhibits the first maximum very close to 3a 1 , as expected for the commensurate phase.
At the temperature of 58 K, i.e, above the upper limit of the commensurate phase stability, the argonargon distribution function recorded at µ Ar = −9.0 demonstrates the presence of a liquid phase, while at µ Ar = −8.8 it corresponds to an incommensurate solid phase. The liquid is of course partially ordered due to the effects of periodic corrugation potential.  Figure 14. The phase diagram for the system with µ Xe = −18.6. Parts (a) and (b) show the temperatureargon chemical potential and the temperaturetotal density projections, respectively. The abbreviations G, L, C and IC stand for the gas, liquid, commensurate solid and incommensurate solid, respectively. The vertical dash-dotted line marks the temperature above which the commensurate phase does not appear.
One should also note a gradual increase of the critical temperature resulting from the increase of xenon concentration. The phase behavior changes when the chemical potential of xenon is increased to −18.6. Figure 14 shows that the gas condenses into a liquid phase of rather high xenon concentration, ranging between x Xe ≈ 0.67 at T = 56 K and x Xe ≈ 0.44 at T = 66 K. When the argon chemical potential increases, we observe the transition between the liquid and commensurate phases. This transition, quite well illustrated by the change in the behavior of the argon-argon radial distribution function given in figure 15, occurs only at the temperatures lower than 61 K. A further increase of the argon chemical potential does not lead to the transition between the commensurate and incommensurate solid phases, as in the previously considered cases, but rather again to the liquid phase. The liquid undergoes a transition into the incommensurate solid-like phase at still higher values of the argon chemical potential (cf. figure 15). This re-entrant behavior can be understood by taking into account that the upper limit of the film density in the commensurate phase is equal to 1.0, while the transition into the incommensurate solid in the argon rich film and at the temperatures used occurs at the densities well above unity. An increasing chemical potential leads to a gradual removal of xenon, so that the dense film becomes more and more argon-like. Note that the liquid-incommensurate solid transition in pure argon film at T = 56 K occurs at the density of about 1.07 (cf. figure 1 (b)), and at still higher densities at higher temperatures.
Concluding, we would like to present the comparison of our results with the available experimental data for submonolayer films of the total density equal to ρ c = 0.4. One readily notes a qualitative agreement between Monte Carlo and experimental results. However, the present simulation predicts a considerably narrower range of xenon concentrations corresponding to the stability region of the com- mensurate phase. Unfortunately, we cannot propose any reasonable explanation for the underestimation of the commensurate phase stability by computer simulation. One can speculate that our model based on Lennard-Jones potential and standard mixing rules overestimates the tendency towards demixing in submonolayer films. The commensurate phase is mixed, while the incommensurate xenon-like phase is demixed. It is also possible, however, that x-ray diffraction data overestimate the range of xenon mole fractions corresponding to the commensurate phase. Note that within the region of coexisting commensurate (C) and xenon-like incommensurate (IX Xe ) phases the paches of incommensurate phase may be quite small and hence escape detection. We recall that Villain and Moreira [35] have also questioned the reliability of experimental results given in reference [23] and suggested that the results were affected by metastability effects.