Self-assembled monolayer formation of pentamers-like molecules onto FCC(111) surfaces: the case of curcuminoids onto Au(111) surface

The adsorption of rigid straight electrically polarized pentamers over a FCC(111) surface is studied. The model was inspired by the deposition of 2-thiophene molecules over the Au(111) surface, which was previously characterized by experimental techniques and simulated under the frame of the density functional theory. We now obtain and report the charge distribution of the molecule which allows to propose a deposition model followed by Monte Carlo simulations over an ad-hoc lattice gas model. We show that for a certain value of the chemical potential there exists an isotropic-nematic phase transition which can explain the formation of a self-assembled monolayer like the one observed in the transmission electron microscopy images. An order parameter is defined to characterize the transition which presents a step-like behavior at a critical chemical potential value. The possible nature of the nematic transition in conjunction with an ergodicity breakdown is discussed as future work by means of statistical physics techniques.


Introduction
The coating of materials surfaces is a standard technique widely used either to protect them or to tune specific properties or functionalities. In recent years an improvement of this technique consisting of the controlled deposition of specific molecules to get surface functionalization has been developed [1,2]. This fact has opened the path for applications in several fields, and encouraged an intensive collaboration among different disciplines like materials science, organic chemistry, biology and condensed matter physics [3][4][5][6]. Thus, the deposition of specifically designed molecules over inorganic surfaces has been proposed for new technologies in bio-medicine [7], nano-electronics [8,9], lubrification control [10], and energy storage [11,12].
Furthermore, when some molecules are adsorbed over inorganic surfaces they can self-organize, i.e. they present a spontaneous ordering [13][14][15]. The proper understanding of this behavior is crucial to extend the coating to the field of nanotechnology [16,17]. In particular, self-assembled monolayers (SAMs) of specific molecules are formed onto metallic surfaces, which suggests that they can be employed to adjust specific electronic and transport properties by tunning the competition between the molecule-substrate and intermolecular interactions.
Among SAMs, the better known are those formed by thiols and dithiols on different oxide-free metals [18][19][20] and on semiconductors [21,22]. They have been studied by many different and complementary techniques, where scanning probes microscopies are the most used evidences [23,24]. Moreover, these techniques are usually complemented by different electron spectroscopies. In this context, gold surfaces covered by thiol and dithiol SAMs have attracted considerable attention, mainly due to both their relatively easy preparation and high stability [25][26][27]. Additionally, the experimental data suggests that the self-assembly Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. process presents different stages, which are apparently governed by the balance between intermolecular interactions and molecule-substrate interactions.
In particular, the adsorption of 2-thiophene curcuminoid molecules on the Au(111) surface has been evidenced by x-ray photoelectron spectroscopy (XPS). According to the literature the obtained XPS signals can be attributed to molecules in a configuration 'lying down' onto metallic surfaces [20]. Furthermore, scanning tunneling microscopy (STM) images of these samples show atomic flat areas on top of the Au(111) surface which are compatible with a self-assembled monolayer [27,28]. From a theoretical point of view, some attempts have been performed under the frame of molecular dynamics [23,29,30], however this approach has two important limitations: an empirically adjusted force field is needed for each interatomic interaction, and computer simulations cannot be performed at low and intermediate molecule concentrations [31]. However, we should point out that very recently a new approach has been proposed which could reduce the effects of this limitation [32].
These observations have motivated us to study the adsorption process by an alternative method, the lattice gas model combined with Monte Carlo simulations, in the aim of elucidate the possibility and/or the conditions for a SAM formation.
Lattice-gas models have been extensively investigated in the last decades because they provide a theoretical framework for the description of many physical, chemical, and biological systems. The adsorption thermodynamics and the understanding of surface phenomena have been greatly benefited from the development of these models [33][34][35]. Particularly, the two-dimensional (2D) lattice-gas model with repulsive interactions between the adparticles has received considerable theoretical and experimental interest because it provides the theoretical framework for studies of surface phase transitions occurring in many adsorbed monolayer films [36][37][38][39][40][41][42][43][44][45][46][47].
To simulate self-assembled adsorption monolayers, the lattice-gas model is also often used. Self-assembly of molecules with complicated chemical structure on a solid surface in most cases is determined by directional intermolecular interactions such as hydrogen bonding, dipole-dipole interactions, and coordination interactions. There are several papers on modeling of such systems using the lattice-gas model [48][49][50][51].
In the case of self-assembly and orientational phase transitions, which is the topic of this paper, an interesting lattice model was introduced by Tavares et al [52], where effectively attractive patches induce the reversible selfassembly of particles into chains. From the work by Tavares et al [52], several papers exploring the selfassembled rigid rods model have been published [53][54][55][56][57][58][59]. These studies demonstrated that the application of lattice-gas models, combined with Monte Carlo simulations, can be very useful to describe self-assembled monolayers in the presence of orientational phase transitions.
We propose here a lattice gas model where the basic input data are the interaction energies: molecule-surface and molecule-molecule, which were obtained in the frame of density functional theory (DFT). The Au surface corresponds to a triangular lattice, indeed to the FCC(111) surface as a general case, while each molecule can be modeled as a linear k-mer. A linear k-mer is defined as an entity that occupies k consecutive empty lattice sites when it is adsorbed. For simplicity we assume that the k-mer is commensurate with the lattice.
In this paper we adopt an ad-hoc lattice gas model for linear pentamers adsorbed onto a triangular lattice, including intermolecule interactions. By simulating the deposition process under a Monte Carlo scheme, we identified an isotropic-nematic phase transition which indicates the formation of a self-assembled monolayer. An order parameter is proposed to characterize the nematic transition.

Molecular interaction
The deposition of 2-thiophene molecules on a Au(111) surface was reported in previous works [27,28]. From a theoretical point of view, the adsorption energy (E ads ) and the corresponding final disposition of a single molecule onto the gold surface were obtained by performing density functional theory (DFT) calculations using the QUANTUM ESPRESSO package [60] with ultrasoft Perdew-Burke-Ernzerhof pseudopotentials [61] and the Car-Parinello dynamics [62] (more details are given in [28]). In this way and after the relaxation of the atomic positions, it was shown that a single 2-thiophene molecule is adsorbed laying down along one of three main equivalent directions imposed by a triangular lattice as those corresponding to the Au(111) surface (see figure 1(a)).
For an isolated 2-thiophene molecule adsorbed along one of the possible easy lying down directions, it is illustrative to obtain the density charge difference ρ diff induced by the molecule-surface interaction. This quantity corresponds to where r molec surf is the charge density for the molecule-surface system in equilibrium, while ρ surf and r molec are the charge densities obtained independently for the surface and for the molecule, respectively. In other words, ρ diff indicates how the molecular electronic charge rearranges after deposition. In figure 1(b), displayed with the help of the software XCrySDen [63], a positive (negative) charge difference is depicted in red (blue). This modified charge density explains the repulsive interaction between two molecules adsorbed over the gold surface. Even more, the interactions are not homogeneous along the molecules so configurational dependent interactions can be expected. As a first approach to the enormous resulting configurational space, 16 relative orientations of two adsorbed molecules were considered and their corresponding energy calculated in a previous work [28].
Although such set represents a small fraction of the configurations space, it allows to identify three predominant ways of interaction. For clarity in this presentation, in figure 2 we have reproduced with the help of VESTA code [64] some representative configurations for these molecules with their corresponding interaction energies E int . Indeed, by defining w HH =4.5eV as the energy corresponding to a head-to-head interaction, we observe in figure 2(a) that there exist two head-to-head interactions and the corresponding energy is about = E w 2.0 int HH . However, if the molecules are disposed as is shown in figure 2(b), the interaction is mainly due to the oxygen atoms in each molecule and it can be approximated by w 1.5 HH , meanwhile it can be attributed   Having in mind these observations, in the present work we attempt an explanation for the SAM formation after the molecule depositions taking into consideration the intermolecular interactions.

Lattice-gas model
In order to simulate the adsorption of 2-thiophene molecules over the Au(111) surface, we have employed a rather typical lattice-gas model. Under this framework, the Au(111) surface is modelled as a 2D rhombusshaped triangular lattice of size M=L×L (L is the linear size of the lattice) with periodic boundary conditions (PBC).
The linear molecules under study are not homogeneous and, as the DFT calculations show, they have segments which interact with different intensity. This statement is explained by the charge rearrangement along the molecule (see figure 1(b)), which is induced by the interaction with the surface of sulfur atoms at both ends, as well as by the carboxyl group at the middle of each molecule. In this line, the 2-thiophene molecule is modeled as a rigid rod that occupies 5 adsorption sites when it is adsorbed, i.e. a linear k-mer of 5 consecutive segments where the distance between adjacent segments equals the lattice constant. Thechoice of a k-mer with k=5 (or pentamer) to model the thiol molecules responds to two facts: (1) It is the smallest and simplest structure that manages to closely reproduce the number of adsorption sites occupied by a 2-thiophene curcuminoid molecule when deposited over a triangular lattice (like the one formed by the Au atoms on the 111 surface) according to the previous DFT calculations. And (2), and more importantly, an entity with 5 segments allows us to define different interaction centers in order to capture the high-polarity nature of these 2-thiophene curcuminoid molecules, as we will detail below. This is a simplification encouraged by the charge distribution presented in figures 1 and 2, where there is one center, two similar ends and two similar connecting sectors, thus a total of five distinct sectors along the k-mer. It is also well known that as the longitude of the adsorbates grows, reaching an equilibrium state is not a trivial task in this kind of studies and the computational time is severely affected by this property even in the absence of repulsive lateral interactions. Taking this last point into account, a pentamer is the smallest k-mer possible of adequately reproducing the properties of the adsorbate without compromising too much the computational cost.
In addition, we consider reversible monolayer adsorption, i.e. the molecules are adsorbed and desorbed of the surface as a whole (no dissociation, reorganization or distortion is allowed) until some equilibrium state is reached, and a lattice site can host only one segment of a molecule.
In order to write the Hamiltonian of the system, the occupation variable c i is introduced, where c i =0 if adsorption site i is empty and c i =1 if it is occupied.
Here, w (expressed in k BT units, where k B is the Boltzmann constant) is the corresponding interaction energy, á ñ i j , stands for the nearest-neighbor sites and N is the total number of k-mers adsorbed. The term ( ) -N k w 1 is subtracted since the first summation overestimates the total number of interactions by adding the internal bonds of a k-mer. The summation of the last term runs on every one of the M lattice sites and therefore, it represents the adsorption energy of the whole lattice, being ε 0 the adsorption energy of each site. Since the lattices considered in this work are homogeneous, we set ε 0 =0 without any loss of generality.
Regarding the interactions between k-mers, we have modelled them in the following way: the pentamers are considered to have 3 interaction centers, both ends ('heads'), and the central segment. The heads are of the same kind, while the segment at the center is of a different kind. Such distribution leads to 3 different basic interactions, head-head (HH), head-center (HC) and center-center (CC). A schematic representation of this model is shown in figure 3. Only nearest-neighbor interactions are considered ( º w w w w , , Before leaving this model section it is necessary to clarify that, the triangular lattice proposed to model the real system within the lattice gas model is aimed to adequately reproduce the preferred relative orientations of the adsorbed molecules, since this is the most relevant feature in our model (see figure 4). The adsorption sites of the modeled triangular lattice do not necessarily coincide with the centers of the Au atoms or with the empty spaces between them as conmensuration is not perfect.

Monte Carlo simulation scheme
The systems investigated, modelled as described in the previous section, were studied under the framework of Grand Canonical Monte Carlo simulations. We implemented a typical adsorption-desorption algorithm following the Metropolis scheme [65], using the Parallel Tempering simulation method [66][67][68]. A given number N rep of replicas of the system are generated, each one at a different value of chemical potential between μ 0 and μ f . The difference in chemical potential between a given replica of the system μ i and its immediate neighboring replica m + , where ΔN and Δμ represent the difference in the number of adsorbed molecules and the difference in chemical potential between the two interchanging replicas respectively. This algorithm allows the system to reach equilibrium in a considerably shorter time than standard sequential algorithms. It is particularly helpful to unblock 'freezing states' where the system is occasionally trapped into local energy minimums.  The general outline of the algorithm is resumed next: Initially, a random configuration for every one of the N rep replicas is generated. The first step of the algorithm consists of randomly selecting one of the N rep replicas, and then a linear k-uple (a set of k consecutive lattice sites) belonging to that replica is selected at random. If the selected kuple is empty (no k-mer adsorbed on any of those lattice sites) a k-mer is adsorbed with probability where ΔH represents the difference between the Hamiltonians of the final and initial states. If, on the contrary, the selected k-uple is occupied (a k-mer is already adsorbed on those lattice sites) the k-mer is desorbed with the same probability P. After repeating M=L×L attempts of adsorption/desorption, a replica exchange step is followed and the configurations of two adjacent replicas are swapped with probability In order to measure the ordering of the adsorbed phase and characterize the Isotropic-Nematic phase transition, we introduce an order parameter defined as, where N i is the number of adsorbed k-mers in the direction i of the triangular lattice and  v i is an unit vector along the direction i.

Results
We begin by considering a triangular lattice having M=50×50 adsorption sites. The corresponding adsorption isotherm and nematic order parameter as a function of the chemical potential μ are shown in Figure 5. Adsorption isotherm and Nematic order parameter as a function of the chemical potential μ for a triangular lattice of 50×50 lattice sites.
figure 5. As expected, the surface coverage rises as the chemical potential increases, but there are two changes of slope at μ≈8 and μ≈13 where the coverage changes abruptly. The first step is clearly related to the isotropic-nematic phase transition where the adsorbed molecules spontaneously align along one of the three possible preferential directions. This orientational phase transition is evidenced by the order parameter which jumps from δ≈0.15 to d  1.0 at μ≈8. The second step of the isotherm at around μ≈13, is due to a reorganization of the already aligned molecules in order to optimize the packing, space is produced to accommodate the incoming molecules reaching very high coverage values θ≈1 at the cost of slowing down the dynamics. Eventually, this reorganization can be accompanied by a switch of ergodic valley in the configuration space, namely, the orientation axis can change as it will be discussed later. This reordering is not evidenced by the order parameter since the system is already in the nematic phase and δ is dynamically normalized with respect to N, the number of deposited k-mers.
In order to get more insight about the evolution and properties of the adsorbed layer, we have obtained snapshots of the system directly from the MC simulations. Figure 6 shows snapshots of three different states of the system along with the adsorption isotherm as a reference. At intermediate coverage the system is isotropic i.e. the molecules are adsorbed along any of the three directions of the lattice with no preferential orientation. This regime is evidenced by the low values of the order parameter for that chemical potential. Increasing the chemical potential, the system goes trough the isotropic-nematic phase transition and the adsorbed molecules spontaneously align along a given direction, as can be seen from the snapshots at μ=11. For higher values of μ, as the snapshot for μ=17 shows, a collective reordering can occur allowing the occupation of almost all the empty sites. This snapshot shows the state of the system at this high concentrations where it can be seen that the molecules are adsorbed forming a highly ordered structure, resembling those of self-assembled monolayers similar to the ones observed in the experiment [28]. This nematic transition is triggered by an ergodicity breakdown which tends to select only one of the three equivalent deposition directions.
We now extend the analysis to larger systems. Thus, in figure 7, we present the results obtained from simulations corresponding to triangular lattices of different sizes (50×50, 100×100 and 150×150 lattice sites) run over 20×10 6 MCs (10×10 6 MCs to equilibrate plus extra 10×10 6 MCs to collect averages). Figure 7(a) concentrates on the order parameter v/s the chemical potential showing the strong jump just over μ=7. Regarding the isotherms ( figure 7 (b)), we can see that the recognition of the nematic transition found for the system M=50 is well reproduced for the larger systems at almost exactly the same values of the chemical potential (just over 7). This feature has been previously observed [69] and it is a reaffirmation of the nematic phase transition observed for the initial system.
At higher chemical potential values differences of two kinds are found. First, δ(μ) decreases its value to almost 0.65 near μ=13 for L=150, while it remains at 1.0 for the smaller lattice sizes. This is due to slow dynamics effects which traps the systems at intermediate metastable states which are in between two or three ergodic valleys. We will come back to this point later.
Second, as already discussed in the case of figure 5 order parameter δ is not appropriate to discuss the second phase transition so we turn our attention now to figure 8 where μ is plotted v/s θ for L=50, 100, and 150; we also include δ(θ) for L=100 just as a reference. As it can be seen the three curves for μ(θ) coincide exactly up to μ≈10 showing the change of slope as δ jumps. Then, the second change of slope is clearly shown at μ≈13 for L=50 and to a smaller degree for L=150; the curve for L=100 does not present any change of slope. This is again a consequence of the slow dynamics for large lattices at these high coverage values. Moreover, this is reinforced by the fact that neither μ(θ) for L=100 or L=150 reaches the saturation value in spite δ does it already.
Let us go back to the behavior of the order parameter for L=150 near the second transition. Figure 9 presents the evolution of parameter δ as μ increases near and over 13. Clearly, patches showing alignment along two or three different directions are observed. Eventually larger computer times would be necessary to unblock these mixed metastable states in the case of large enough systems.
At the moment, it is not clear if this behaviour near μ=13 corresponds to a second 'entropic' phase transition (partially detected by the nematic order parameter for L=150) or if this behaviour is merely the result of a non-equilibrated system at high coverage values. In order to properly detect this second phase transition it is necessary to use different ways of characterizing and measuring it. We refer to finite-size scaling  techniques and information theory tools, which is beyond the scope of the present paper and will be studied and presented separately.

Conclusions
In the aim of understanding the self-assembled monolayer formation of 2-thiophene molecules over Au(111) surfaces, we have studied this system from a theoretical point of view. By extending previous density functional theory calculations, we define three equivalent easy lying directions for the adsorption of any molecule over the gold surface. Moreover, we identify three kinds of basic repulsive interactions among neighboring molecules. In this way we can define a lattice gas model consisting of polar linear pentamers adsorbed over a 2D rhombusshaped triangular lattice.
By performing Monte Carlo simulations we have identified an isotropic-nematic phase transition occurring at lattice coverage of around θ≈0.8. When the system goes through this phase transition, the adsorbed molecules spontaneously align over one of the three main directions of the lattice. Such phase transition indicates that this may be a mechanism for the 2-thiophene molecules to form a self-assembled monolayer over the Au(111) surface in agreement with previously reported experimental results already discussed in the Introduction.
Additionally, we have found some evidence favoring a second phase transition at high concentrations, which may be entropic in essence, and it is related to a redistribution and increase in the surface coverage. Some initial studies of us have indicated that ergodic separation plays an important role in the separation of these phases. Such, second phase transition will require another way of measuring it, eventually another order parameter, but still, it will be hard to detect since equilibrium problems dominate at high values of lattice coverage. For larger systems these equilibrium problems are even more important bringing in huge computational costs. At this moment we point out this unusual phenomenon but it is left as an open question.
The whole methodology exposed in this work can be extended to understand and/or follow the selfassembly process for other kind of molecules deposited over different substrates. For example, thiol-based molecules adsorbed on gold produce a rich structural variety of SAMs which strongly depend on the intermolecular interactions [70]. In a similar way, N-heterocyclic carbenes deposited over gold surfaces can form ultra stable SAMs, which have promissory technological applications [71]. On the other hand, surfaces as metal oxides can also support SAMs, where phosphonic acids and perylenediimide are two cases having interesting projections for electronics [72] and for energy storage [73] industries. Another interesting system, which deserves a theoretical study, corresponds to silicon oxide micro-spheres forming SAMs on a sode-lime glass, with promising applications in cooling down devices, useful for green technologies [74].
Finally, the simplified lattice-gas model presented in this work has shown that it reproduces the phenomenological behavior of the adsorption process of the 2-thiophene molecules, particularly, the formation mechanism of a SAMs. In this way, it has proven to be an adequate tool to characterize this kind of adsorption phenomena where a high number of parameters needs to be taken into account. Future efforts will be done following two directions: first, to extend this lattice-gas treatment to other experimental systems where SAMs have been observed and second, to perform a more rigorous study of these phase transitions in the framework of statistical-mechanics, finite-size scaling and information theory.