#### Abstract

Applying the reptation algorithm to a simplified perfluoropolyether Z off-lattice polymer model an NVT Monte Carlo simulation has been performed. Bulk condition has been simulated first to compare the average radius of gyration with the bulk experimental results. Then the model is tested for its ability to describe dynamics. After this, it is applied to observe the replenishment of nanoscale ultrathin liquid films on solid flat carbon surfaces. The replenishment rate for trenches of different widths (8, 12, and 16 nms for several molecular weights) between two films of perfluoropolyether Z from the Monte Carlo simulation is compared to that obtained solving the diffusion equation using the experimental diffusion coefficients of Ma et al. (1999), with room condition in both cases. Replenishment per Monte Carlo cycle seems to be a constant multiple of replenishment per second at least up to 2 nm replenished film thickness of the trenches over the carbon surface. Considerable good agreement has been achieved here between the experimental results and the dynamics of molecules using reptation moves in the ultrathin liquid films on solid surfaces.

#### 1. Introduction

The diffusion of molecules on surfaces is fundamental to phenomena such as heterogeneous catalysis, wetting, self-assembly, and nanofabrication [1]. Recent studies have focused on the behavior of mechanical properties [2, 3], glass transitions [4–11], diffusion in thin films [12–15] and colloidal particles [16, 17], and phase separation [18–20]. Obtaining detailed microscopic-level information on the interfacial structure and mobility of a polymer, and establishing links between this information and overall thermodynamic and mechanical performance, could greatly aid in the improvement of existing polymeric coatings and the development of new ones.

Spreading of molecularly thin films is key to many physically and industrially important processes that rely upon thin film adhesion and lubrication [21–23]. One specific industrial application is the use of submonolayer perfluoropolyether (PFPE) films to lubricate magnetic recording media (computer hard disks) within magnetic storage devices (hard disk drives). Reflow of the PFPE lubricant on the hard disk is vital to the reliability of disk drive [24]. The importance of lubricant mobility in hard disk drives motivated the initial studies investigating the spreading of PFPEs on silicon surfaces. These studies were done on PFPE Z with nonpolar end groups [25] and the spreading profile exhibits a smooth, monotonically increasing thickness profile. A rapid spreading is observed, and the profile is asymmetric about the position of the original step [26]. The shape of the spreading profiles for the PFPE and polydimethylsiloxane (PDMS) systems is derived from the film thickness dependence of the diffusion coefficient through the disjoining pressure dependence on film thickness [25, 27–29]. The diffusion coefficient derived from the dispersive force is a monotonically decreasing function of film thickness [27]. This component of the diffusion coefficient can account for the rapid spreading front and asymmetry in the spreading profile [27, 28].

Experimental diffusion coefficient has been calculated recently from the spreading profiles of PFPE Z on amorphous carbon surfaces [30] using the Matano interface method [31]. The thickness-dependent diffusion coefficient was found to be independent of time and initial film thickness, indicating a time dependence of the diffusion process. In this study a simplified, coarse grained, freely bending, and bead-spring model of PFPE Z is introduced to simulate replenishment of ultrathin liquid films on solid carbon surfaces. This is a novel approach to simulate the dynamics in the ultrathin liquid films from the molecular perspective and to exclude the necessity of expensive experiments by using direct Monte Carlo (MC) simulations.

#### 2. Model and Monte Carlo Simulation Method

The two most important simulation methods are Monte Carlo and molecular dynamics. Both methods can be used for the calculation of average values in equilibrium thermodynamic ensembles [32, 33]. In MC simulations a chain of configurations is generated (a Markov chain) in a stochastic way, by generating certain trial moves and accepting them with a certain probability. The choice of trial moves depends on the nature of the system that is simulated. To obtain maximum efficiency the correlation between successive configurations is made as small as possible.

MC simulations of nonequilibrium or relaxation processes are based on the interpretation of the Markov chain as the evolution of the system in time [34, 35]. The correlations between the successive configurations in the Markov chain are understood in terms of dynamic correlation functions with stochastic kinetics. Therefore the configuration changes have to correspond to real events in the stochastic system. Every type of move in the system has its own frequency reflecting the number of events of that type performed per unit time. This approach does not provide an absolute measure of physical time, only a proportionality.

When particle (here, bead) displacements are the only moves, one usually expresses time as the number of sweeps performed. Only by measuring a quantity during the simulation that is proportional to the real time, a better estimation of the physical time passed can be made. Huitema and van der Eerden [36] recently have shown that MC simulations can be used for simulating dynamical processes. They concluded that processes with rate-determining steps at time scales larger than about 1 ps and length scales larger than 0.40 pm can be simulated with MC resulting in similar dynamics as compared to MD. According to them these values apply to liquid Lennard-Jones systems with the usual argon fitted values of the interaction parameters and masses. In general the appropriate length scale will be on the order of the mean free path per molecule (about 0.1 for Lennard-Jones [37]) and the time scale will be roughly this length divided by the thermal velocity .

Recently reptation quantum Monte Carlo which can exploit the dynamical properties of the classical diffusion process—rather than retaining the asymptotic distribution alone—and map them onto the (imaginary-time) dynamical properties of the quantum system of interest is proposed as a method for unbiased ground-state averages and imaginary-time correlations by Baroni and Moroni [38]. They concluded that the reptation quantum Monte Carlo can be applied to clusters, films, and superfluids in restricted geometries and the algorithm is free from the mixed estimate and the population control bias. Furthermore, because it samples an explicit expression for imaginary time evolution, the algorithm gives access to quantities obtained by differentiation.

In this study, reptation algorithm [39] which is suitable for a chain with arbitrary intermolecular and intramolecular potentials in a continuum fluid is modified and applied. The method exploits the Metropolis solution to the transition matrix to asymptotically sample the Boltzmann distribution.

Chemical formula of PFPE Z considered is F-CF_{2}-[(OCF_{2}-CF_{2})_{m}-(OCF_{2})_{n}]-OCF_{2}-F where . The variation of this formula for molecular weights (MWs) of 3840, 2500, and 1700 g/mol is shown in Table 1. In this table, a polymer molecule is assumed by a chain of several beads and each bead is made of several monomer units of OCF_{2} and CF_{2}. The background assumptions on the selection of simulation beads are as follows.(a)CF_{2} is the basic backbone unit and O is distributed on the backbone uniformly.(b)CF_{3} is same as CF_{2} at the two ends.

According to the chemical structure of PFPE Z bonds are more rotatable around the O atoms on the chain. Therefore, the point of view of the rotatable units still persists with these assumptions on the selection of simulation beads. Moreover, it is thought that the simulation becomes more simplified with these assumptions which are thought to be justifiable at least in case of PFPE Z. Cross-sectional diameter of all the polymers is considered to be 0.7 nm and the lengths of the polymer molecules of MWs of 3840, 2500, and 1700 g/mol are taken to be about 13.44, 8.75, and 5.96 nm, respectively. Effective bond lengths (EBLs) are calculated on the basis of division of the total polymer lengths by the approximate divisor of Table 1 and are 1.12, 1.25, and 1.2 nm for MWs of 3840, 2500, and 1700 g/mol, respectively, which are about 1.5 times the cross-sectional diameter and are consistent with the ability of coil formation with the bulk radius of gyration as illustrated in Figure 1 for MW of 3840 g/mol.

##### 2.1. Interaction Potentials

Three different potentials are considered in this model as a polymer molecule is assumed by a chain of several beads. All beads have the dispersive interaction described by the truncated Lennard-Jones () potential, . Each polymer chain consists of a sequence of beads bonded by the potential energy of the limited harmonic potential, . For the dispersive bead-wall interaction, a 3–9 potential, , has been used.

Lennard-Jones () pair potential which operates between all pairs of beads in the system including the bonded ones is where since , which is the interaction cutoff distance [32], characterizes the strength of the interaction, and is a characteristic length of the interaction potential which is taken to be equal to the EBLs. This is an easy way to include excluded volume effects by putting spheres centered at each connection point on the chain [40].

For MC simulations, a slightly different choice of potential which links the adjacent beads in the same chain was recently proposed [41, 42]:
where since which is the attempted displacement of the head or the tail bead and is a constant. Also in the simulation , , , m^{−2} ( is the Boltzmann’s constant), and are considered because with this choice of parameters chain intersection is practically forbidden even if one uses large attempted displacements of the head or the tail (in each update, the head or the tail is attempted to be displaced into a randomly chosen direction and the step width is chosen randomly from a cube with linear dimensions equal to , the old position being in the center of the cube) [40].

The substrate is modeled by a flat continuum material and thought to be homogeneous. The total substrate-bead interaction potential which is obtained by integrating the potential over the half space [43] is where , , , nm (carbon atomic diameter), and and (constant) = 3000 (approximately) are considered. It can be shown that, for the substrate-fluid interaction, the depth of the potential is comparable to and the width is comparable to [44].

##### 2.2. Initialization

The length of the box is considered to be side , side nm, and side to disable the beads to sense the symmetry due to the periodic boundary condition in the direction for the fluid of beads in the box [32]. All the beads in a polymer are randomly inserted and the insertion of the total polymer is rejected if the randomly selected attempted insertion position of any of its beads is preoccupied. Periodic boundary condition at the directions, neutral wall condition at the directions, and solid boundary condition at the bottom of the box are then activated. The top surface of the box restricted any move of the polymers above it and the bottom surface acted on them with potential of (3). Moreover, all the interaction potentials as mentioned above are then activated. MC thermalization is then done until the system reaches a near to equilibrium condition which has been detected by very small changes in the polymer conformation and distribution by giving further runs. Though gravitational potential has been considered in this simulation, the effect of this potential is observed to be very negligible in the results.

##### 2.3. Simulation Method

In every MC cycle all the polymers are interrogated. Each polymer is considered randomly and again the head or the tail is chosen randomly to be displaced. Suppose that the initial coordinates of the beads in the th polymer are . Here and are the head and the tail, respectively, and one of them is randomly chosen. A new position of the selected end bead (suppose, ) is then selected such that where the direction of is chosen at random on the surface of a sphere by the acceptance-rejection technique of von Neumann [45] and the magnitude of is chosen according to the probability distribution , where and is given by (2). In this way, the model includes the bond angle vibrational and bond torsional motion of the polymers randomly in the simulation. The polymer chain has a new trial configuration after the bead positions are shifted along the chain contour. The change in the total interaction potential is calculated as follows:

The trial configuration is accepted with probability equal to

If the move is accepted the polymer beads take the new configuration and, on the contrary, the polymer retains its old configuration. In any case the next randomly selected polymer is then considered for interrogation and the selection processes are repeated.

After thermalization of 6000 MC cycle (MCC) at least two thermalized films are placed side by side in the direction keeping trenches of different widths (8, 12, and 16 nm for molecular weight (MW) of 3840 g/mol and 12 nm for MWs of 2500 and 1700 g/mol) between them separately as shown in Figure 2 and the middle mutual facing two neutral walls of the two films are then taken away to allow the films to spread and replenish the trenches. Dynamic equilibrium condition of the replenished trenches is traced by averaging over number () of different dynamic equilibrium distributions obtained in sequence and separated by 100 MC cycles (one MC cycle consists of the attempted MC moves equal to the number of polymers taken) from the initial static equilibrium condition. The number of polymers taken in the previously mentioned box is 133, 220, and 332 for MWs of 3840, 2500, and 1700 g/mol, respectively (which should be multiplied by the number of beads per polymer, (13, 8, and 6 for MWs of 3840, 2500, and 1700 g/mol, resp.) to get the similarity with the number of atoms) and taken is 4. Higher has also been taken to verify any change in the averaged results but no significant difference is observed.

#### 3. Results and Discussion

##### 3.1. Bulk Radius of Gyration

At first simulation has been performed in bulk situation where periodic boundary condition has been used in all the three coordinate directions of the simulation box. The lengths of the box in all the three directions are considered to be side = side = side to disable the beads to sense the symmetry due to the periodic boundary conditions in all the coordinate directions for the fluid of beads in the box [32]. Only bead-bead Lennard-Jones () potential of (1) is used in the bulk situation since there is no substrate surface in this case. Equilibration has been performed in the same way as stated above in the case of ultrathin liquid films on solid surfaces. In the bulk simulation, it is observed that the average radius of gyration values with all the MW cases (i.e., 3840, 2500, and 1700 g/mol) is nearly equal to their bulk experimental values [26, 46]. According to [46], the experimental radii of gyration of PFPE Z for MWs of 3840, 2500, and 1700 gms/mol are 1.35, 1.07, and 0.87 nms, respectively. From the MC simulations the average radii of gyration of PFPE Z for MWs of 3840, 2500, and 1700 gms/mol are 1.38, 1.10, and 0.91 nms, respectively. The average components of* Rg* in the , , and directions,* Rgx*,* Rgy*, and* Rgz*, are almost equal to each other in the bulk denoting the average shape of a polymer as a sphere.

##### 3.2. Mean-Square Displacement of the Center of Mass of the Chain and the Self-Diffusion Constant

The self-diffusion constant can be determined from the mean-square displacement of , the center of mass of the chain,

Within the Rouse model for all times which can be observed in Figures 3 and 4, the diffusion constant is defined as where is the degree of polymerization.

Even more importantly the data are accurate enough to display the slowing down in the motion of the center of mass of the chains as shown in Figure 5. This phenomenon was also found in case of bond fluctuation MC [40]. After the initial effective slope of nearly 0.9 it seems to decay to approximately 0.5 and then it again increases to a value nearly to 1, as expected by the reptation model [40].

According to the simple Rouse model, we should have independent of time. This is true only for density, , while for even smaller densities (and of course larger densities like that of melt) one can see an initial decrease of with , until settles down on a kind of plateau, which hence defines the self-diffusion constant given by (7). Therefore, has been plotted with in Figure 6 for several MWs and approximate self-diffusion constants; have been obtained as shown in the same figure. The plots in this figure have fluctuations; those may be due to the selection of larger effective bond lengths (EBLs) to model the polymers in this simulation. The larger EBLs are responsible for the suppression of the details of the scale less than the EBLs. However, larger EBLs along with the artificial reptation moves decrease the computational time considerably to access the experimental results. Moreover, to model a relatively complicated molecule of PFPE we had to divide a molecule into beads as shown in Table 1 for different MWs and consider relatively larger EBLs to keep the chemical composition of each bead same. We have plotted the log-log plot of these values with respect to the corresponding MWs in Figure 7 and the data are compatible with a simple (a Rouse model prediction), where is the degree of polymerization that is equivalent to MW multiplied by a constant.

##### 3.3. Present Study on Replenishment

Diffusion coefficients of PFPE Z have been calculated recently by Matano interface method [31] from the experimental spreading profiles of ultrathin liquid films of the same on amorphous carbon surfaces. The diffusion equation considering the disjoining pressure as the motivating force has also been proposed [46] as the equation to simulate spreading and replenishment of the ultrathin liquid films on solid surfaces. In this study this diffusion equation has also been solved using the experimental diffusion coefficients of Ma et al. [30] to simulate replenishment of the trenches between the two previously mentioned films from the same initial profiles as used in MC simulations. The system conditions are (a) temperature which is 25°C (MC simulation) and 26°C [30], (b) relative humidity which is 0% in both the cases, and (c) carbon surface roughness which is flat (MC simulation), average surface roughness of 6 Å [30]. As an illustration, replenishment profiles with time solving the diffusion equation for the trench width (TW) of 8 nm with MW of 3840 g/mol are shown in Figure 8. On the other hand, replenishment profiles with number of MC cycles (MCC) using the MC simulation for the same TW and MW are shown in Figure 9 (only the trench cross-sectional view has been zoomed in the figures). In all these figures replenishment curves which are self-explanatory for meeting of the two film edges and replenishment up to several nanometers have been illustrated. In Figure 9 the film is divided into vertical layers of 5 width (slightly lower than the cross sectional diameter of PFPE Z) in the direction and the average number of beads per layer is assumed to be directly proportional to the average thickness of the film, 4 nm. Similarly, 5 Å cell width has been taken in solving the diffusion equation. The solution of the diffusion equation has also been verified with the solution of the same by Ma [46] and is found to be completely in agreement. Comparison of MC simulation results with the solutions of diffusion equation has been made on the basis of the crossing of the lower portions of the replenishment curves from both the methods some nanometer level thickness above the bottom surface. Figures 10 and 11 show curves of number of MC cycles versus time in microseconds needed for the meeting of the two film edges and replenishment up to several (0.5, 1, 1.5, and 2 nanometers) nanometers from the surface for different trench widths and MWs, respectively (TW = 8, 12, and 16 nm for MW of 3840 g/mol in Figure 10 and 12 nm for MWs of 3840, 2500, and 1700 g/mol in Figure 11). In Figure 9 the replenishment curves from the MC simulation show fluctuations; those make it difficult for the analyzer to compare between the number of MC cycles and the real time. Thus, we have included error bars for each point of comparison between the number of MC cycles and the time in microseconds in Figures 10 and 11 and the comparison between the two can be justified within these error bars. All the curves for different TWs and MWs show almost constant ratio of number of MC cycles to time in microseconds that is about 10^{4} considering the average values of the error bars. In a broader sense, this means that the simulation retains the time correlation during the diffusion process of the polymer molecules with respect to a characteristic length of the space (TW) and a characteristic length of the polymers (MW). Obviously, question arises as to whether the same phenomena are obvious to many other polymers or not? Or whether the constant ratio of MC cycle and real time with respect to the above characteristic lengths is universally observable in many other chain polymers?

**(a)**

**(b)**

**(c)**

**(d)**

The result signifies the capability of the newly introduced model of PFPE Z to trace the length and time scales (at least from the point of view of longer time scales) of spreading and replenishment of ultrathin liquid films of the same on solid surfaces by direct MC simulation technique. Furthermore, the model of PFPE Z introduced here introduces a possibility that the time scale of many other phenomena in the nanometer scale may be traced having these results in good agreement with the experiment for different TWs and MWs as a baseline of the time per MC cycle (MCC).

#### 4. Conclusions

The work presented here is on the development of an MC simulation program, which generates average conformation and longer time dynamics of PFPE Z in the ultrathin liquid films on solid surfaces. In the bulk simulation, it is observed that the average radius of gyration values with all the molecular weight cases are nearly equal to their bulk experimental values. Experimental results are compared with the molecular counterparts and the results seem to be consistent considering the average longer time dynamic motion of polymers near the solid surface. All the results for replenishment of ultrathin liquid films on solid surfaces with different trench widths and molecular weights show almost constant ratio of number of MC cycles to time in microseconds. The results reveal that a simple model of PFPE Z is capable of simulating replenishment and reptation MC simulation can retain the time correlation of the same and may be any other polymers.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

This work has been supported partially by the Japan Society for the Promotion of Science (JSPS).