First-order phase transitions in lattice bilayers of Janus-like particles: Monte Carlo simulations

The first-order phase transitions in the lattice model of Janus-like particles confined in slit-like pores are studied. We assume a cubic lattice with molecules that can freely change their orientation on a lattice site. Moreover, the molecules can interact with the pore walls with orientation-dependent forces. The performed calculations are limited to the cases of bilayers. Our emphasis is on the competition between the fluid-wall and fluid-fluid interactions. The oriented structures formed in the systems in which the fluid-wall interactions acting contrary to the fluid-fluid interactions differ from those appearing in the systems with neutral walls or with walls attracting the repulsive parts of fluid molecules.


Introduction
A growing interest has been recently observed regarding the self-assembly of Janus and patchy colloids into differently ordered structures. This is also due to the progress in the development of new experimental methods for obtaining particles that are capable of building blocks of the requested symmetry and structure [1][2][3][4][5][6][7]. One way to prepare spherically symmetric particles with anisotropic interactions relies on the creation of attractive patches on their surfaces [4,8]. Advanced experimental methods allow for the production of patchy colloids with different sizes, shapes and chemical compositions [9][10][11][12][13][14]. In parallel with experiments, remarkable theoretical and computer simulation efforts have been undertaken to predict the physical properties of these systems [15].
An important direction of the research is connected with investigations of phase transitions in such systems. The directional dependence of interactions leads to a very rich phase behavior and to the development of various mesoscale structures [16][17][18][19][20][21][22][23][24][25][26][27]. A review of unusual structures predicted for several types of colloids can be found in [2]. In particular, in the case of patchy spheres interacting via the Kern-Frenkel potential [28], diverse aggregates, such as micelles, vesicles, and bilayers, have been found [21][22][23]. Also, the coexistence of a dilute phase of micelles and denser phase of larger clusters has been observed. More recently [23] simulations of the formation of different crystalline structures have been carried out. Another interesting kind of behavior has been obtained within the framework of Wertheim theory [29] for the patchy colloids with multiple sites, which, in some cases, makes multiple bonding possible between sites [30][31][32][33][34][35].
Lattice models have also been used to study the self-assembly. For example, Dawson and Kurtović [36] introduced a lattice model of amphiphilic self-assembly in mixtures containing amphiphiles and water. Davis and Panagiotopoulos [37] examined mixtures of flexible chains comprising neutral and solvophobic monomers and big amphiphiles with multiple chains connected to them at common attachment points. They found that amphiphile geometry plays a key role in determining whether the micelles are spheres or flat bilayers. The behavior of Janus particles near solid surfaces was also investigated using density functional theory and computer simulations [38,39].
Recently, two of us reported the results of Monte Carlo simulations of the phase behavior of Janus disks on a square lattice [40]. The particles were composed of two different parts and interactions between neighboring particles depended on their orientations. It was found that the systems exhibited the first-order transitions between colloidal-rich and colloidal-poor phases and continuous order-disorder transitions to differently ordered structures.
Calculations performed in [40] assumed a restricted number of possible orientations of the particles. In this work we report Monte Carlo studies of bilayers of Janus-like molecules on a cubic lattice confined in slits, but in contrast to the previous work [40], the molecules can freely change their orientations on the sites. Moreover, the slit walls can also interact with particles with orientation-dependent potential. Our principal aim is to study the first-order transitions in the system and the dependence of the phase behavior on the competition between interparticle and the particle-wall interactions. The principal orderparameter is density. However, we neither carried out systematic studies of the finite size effects [41,42] nor the studies of continuous transitions in the system.

Simulation
We consider a cubic lattice of the size N x × N x × N z , N z = 2 in X , Y and Z directions, respectively. The lattice constant is assumed to be unit of length. Periodic boundary conditions are used in the X and Y directions. In the Z direction, however, the system is closed by planar walls, located at i z = 0 (the bottom wall) and at i z = 3 (the top wall). We consider spherical particles consisting of two hemispheres, one being solvophobic (A) and the other one being solvophilic (R). The solvent is involved implicitly.
We introduce the occupation number for each lattice site, n i , i = (i x , i y , i z ) that equals 0 for an empty site and 1 for an occupied site. The state of a molecule occupying the site i is characterized by a unit vector (sin θ i cos φ i , sin θ i sin φ i , cos θ i ), where θ i is the tilting angle measured with respect to the vector perpendicular to the bottom wall and φ i is the azimuthal angle, measured with respect to the positive X axis.
Due to the existing walls, the energy experienced by a molecule located on the site i is given by To evaluate the interparticle potential energy, we take into account only the nearest-neighbor interactions. If two nearest-neighbor sites, i and j, are occupied, the pair potential energy is where α i and α j are the angles measured with respect to the unit vector pointing from the site j to the site i. A similar potential of anisotropic interactions was used for amphiphilic molecules in [38]. We assume that ε > 0. In the absence of an external field, the particles prefer to couple in pairs whose A-patches face each other. On the contrary, the RR-contacts are energetically penalized.
Our approach can be treated as an extension of [40], where only four possible orientations of particles were considered. In this work, the orientation vector of any molecule can change continuously, 0 < θ < π and 0 < φ < 2π. The Hamiltonian can be written as follows: where µ is the chemical potential. In the above, the first summation is over all lattice sites. The subscript nn means summation over all pairs of the nearest neighbors.
The density is given by where V = N 2 N z . The potential energy per one lattice site, 〈U 〉, and the contributions due to external field, 〈U s 〉, and due to particle-particle interactions, 〈U p 〉 can be expressed as follows:

33602-2
We have carried out Monte Carlo (MC) simulations in the grand canonical ensemble using the hyperparallel tempering technique [43][44][45][46]. The size N x of the system was equal to 48. Besides the systems exchange between replicas at different µ and T , a Monte Carlo step consisted in an attempt to insert a randomly oriented particle into the system at a randomly chosen position, an attempt to remove the existing particle and an attempt to change the orientation of a selected particle. A number of MC steps necessary for obtaining solid results depended considerably on the assumed values of the energy parameters. In the majority of simulations, we used 10 9 MC steps (per site) for equilibration and 10 10 for production runs.
However, for certain sets of energy parameters, the simulations had to be longer. We have calculated the average density and average potential energy and corresponding fluctuations, the compressibility and, the heat capacity, Besides thermodynamic quantities listed above, we have evaluated the histograms of ρ and the histograms of cos θ i z in each layer. The first-order phase transitions were investigated by analysing the density histograms, and the coexistence of different phases was located using the equal peak-weight criterion for the density histograms [44]. Moreover, we have calculated the angular-averaged local densities 〈ρ i z 〉, the average values of the cosine of the tilting angle in each layer, 〈cos θ i z 〉, and the angular-dependent local densities ρ i z (θ, φ).
The energy parameter ε was taken as unit of energy and the reduced temperature is T * = k B T /ε. Similarly, the reduced chemical potential is µ * = µ/ε.

Results and discussion
The aim of this study is to investigate the phase transitions in the systems. The order-parameter was the average density. We changed interactions of the walls with the R-sides of Janus molecules from attractive, through neutral, to repulsive, i.e., t = ε s /ε = 5, 1, 0 and −1 and, in order to distinguish the systems under study, we abbreviated them by the label Et , where t = 5, 1, 0 or −1.
We start the discussion with the case of Janus particles in relatively strong external field (the system E5). The results are presented in figures 1, 2 and 4-6. At low values of µ, the molecules expose their attractive parts towards the pore walls and, consequently, the molecules in bottom, i z = 1, layer expose their repulsive parts towards molecules in the second (i z = 2) layer. The external field acts against an "intrinsic" tendency to form an oriented bilayer due to attractive AA-interactions. In Three jumps on the isotherm in figure 1 (a) are seen. They are marked as I, II and III, respectively. The values of compressibility, κ are discontinuous at these jumps, cf. figure 1 (b). At low values of the chemical potential, µ * , the molecules are oriented almost perpendicularly to the surface (θ 1 ≈ π for the first and θ 2 ≈ 0 for the second layer molecules). Of course, this orientation results from the strong external field. For µ * < −1.6, the formation of clusters composed of four molecules at each surface is observed (see the snapshot in figure 2). The sites over the cluster formed within the bottom layer, as well as the sites below the clusters within the upper layer are empty. In-plane orientation within each cluster is such that 〈| sin(φ)|〉 ≈ 0.7, i.e., the values of the azimuthal angle in each layer for the cluster are, respectively, close to π/4, 3π/4, 5π/4 and 7π, starting from the "low-left" corner of the cluster.
With an increase of the chemical potential, the clusters start to flow together and to form ordered structures. These structures can be considered as patterns built of the occupied and empty lattice sites. In such patterns, orientations of particles can change. There is translational order without a more pronounced orientational order. We have found three translationally ordered phases for ρ = 1/2, ρ = 2/3 and ρ = 5/6. They are schematically depicted in figure 3. When a half of the lattice sites is occupied, the zigzag structures on each wall are observed. There are sloping steplike strips of occupied sites on the 33602-3 We return now to the discussion of phase transitions in the system E5. At the first plateau of the isotherm [figure 1 (a)], the average density is close to 〈ρ〉 ≈ 0.5. An example of configuration at the plateau of the isotherm (µ * = −1.5) is shown in figure 4 (a). Indeed, one sees here a well pronounced zigzag structure. Within this region, the average tilting angles of the top (bottom) layers become close to θ 1 ≈ 5π/6 (θ 2 ≈ π/6), cf. figure 1 (c). At T * = 0.06, the development of the zigzag structure is a continuous transition.

33602-4
Bilayers of Janus-like particles A further increase of the chemical potential leads to the rearrangement of the molecules and to the formation of well-developed "pillars" between the walls. The configuration displayed in figure 4 (a) is stable and a significant increase of the chemical potential is needed in order to enforce more molecules to enter the pore. The first isotherm jump is connected with the filling up of the empty zigzag structure. This filling leads to the location of the particles "one atop the other" and both such particles expose their repulsive parts one towards the other. In order to reduce the repulsive potential energy, the titling angles for the molecules located within the first, θ 1 , and the second layer, θ 2 , deviate more from the orientations perpendicular to the walls. The first step on the isotherm (I) corresponds to the transition from the zigzag structure to the structure containing single strips of pillars (SP).
A typical snapshot of the configuration before the second transition [marked as II in figure 1 (a)] is shown in figure 4 (b). After this transition, the second plateau on the adsorption isotherm develops. The isotherm plateau (the average density is ≈ 0.67) is shorter than the first one. A representative snapshot of the configuration of molecules is displayed in figure 4 (c). The second transition is also connected with an increase of the value of cos θ 1 to its maximum value. At T * = 0.06, the maximum of cos θ 1 is close to −0.83. The step II is associated with the transition from the SP-phase to the phase in which double strips of pillars begin to appear. At the end of the second plateau, the third step on the adsorption isotherm appears. We think that the last step (III) is a transition from the ordered DP-phase to a phase with the filled empty sites that remained. This step is connected with a decrease of cos θ 1 . A further increase of the chemical potential leads to a gradual increase of the density inside the slit.

a b c
At high values of the chemical potential, the pore is completely filled. At the third final plateau, cos θ 1 ≈ −0.94. Instantaneously, the azimuthal orientation of molecules within the filled first and second layers is similar to that observed at very low values of the chemical potential, when disconnected "islands" composed of four molecules have been formed, i.e., 〈| sin(φ)|〉 ≈ 0.7. was difficult and would require extremely long runs to obtain reliable histograms. We also did not attempt a precise evaluation of the critical temperatures. However, one sees that the critical temperatures of the phase transition increase in the order T I < T II < T III [the indices I, II and III abbreviate the consecutive transitions from figure 1 (a)]. The changes of the density during phase transition are small and this is a rather striking behavior. It is also worth stressing that the changes of the chemical potential with the temperature are non-monotonous. We now consider the system E1. Similarly to the previous case, the external field acts against the "intrinsic" tendency to form an oriented bilayer with facing R-sides of the neighboring Janus particles, but now this effect is weaker. In figure 7, we show the adsorption isotherm, the plot of 〈cos θ 1 〉 (T * = 0.06) and the phase diagrams for the system E1.
One sees three jumps on the adsorption isotherm and at very low chemical potentials, the structure of the confined fluid is similar to the system E5. The molecules within each layer form squares. However, since the external field is weaker, the number of contacts between the "square clusters" in bottom and in top layers is larger than for the system E5, cf. the snapshot in figure 8

33602-7
between fluid-fluid and fluid-wall interactions at low chemical potential, the tilting angle in the system E1 significantly differs from the tilting angle in the system E5. Now 〈cos θ 1 〉 ≈ −0.74, i.e., θ 1 is close to 3π/4. The in-plane orientation, however, is similar to the system E5, but the average value 〈| sin φ|〉 is lower, 〈| sin φ|〉 ≈ 0.64. This is the consequence of the existence of contacts between bottom and top layers.
With an increasing chemical potential, a small plateau on the adsorption isotherm is reached at µ * ≈ −0.54 [ figure 7 (a)]. Before the transition, a the zigzag structure is observed. However, as we have already stressed, the tilting angle is different, 〈cos θ 1 〉 ≈ 0.64 [cf. figure 7 (b)]. The transition I almost does not influence the tilting angle. After the transition, the average density of the confined fluid is close to 0.6. The translational ordering of the fluid is developed from the zigzag structure in such a way that the stripes at each wall are "compressed" and they alternately contact the filled and empty stripes at the other wall, cf. figure 8 (b).
The second transition [marked as II in figure 7 (a)] leads to the structure of the density 〈ρ〉 ≈ 0.67 that is shown in figure 8 (c). This transition is connected with a large change of the tilting angle [ figure 7 (b)]. At the isotherm plateau that appears for the chemical potential 0.32 < µ < 0.45, the tilting angle is 〈cos θ 1 〉 ≈ 0.18. This means that the molecules are re-oriented from the configuration with molecules facing their attractive parts towards the walls (before the transition) to the configuration with the attractive parts of the molecules partially directed towards the pore interior (after the transition). A small density jump at the transition III is not connected with any significant change of the tilting angle. For a completely filled pore, the tilting angle θ 1 is close to π/4.

Figures 7 (c) and 7 (d)
show the phase diagram in the density-temperature and the chemical potentialtemperature planes. However, contrary to the system E5, for the system E1 the critical temperature T II is the highest and the critical temperature T I is the lowest. Moreover, the critical temperatures for the system E1 are in general higher than for the system E5, which is the consequence of the total potential energy being considerably lower in the case of the system E5.
We next discuss the system E0, where the pore walls are just hard walls and orientation effects result entirely from the fluid-fluid interactions. Figure 9 shows the isotherm [part (a)], tilting angles [part (b)] and the phase diagrams [parts (c) and (d)]. Additionally, in part (b) we have also displayed the values of 〈| sin φ|〉. As in the previous cases, three first-order phase transitions appear within the investigated range of temperatures. However, we observe the formation of considerably different ordered structures in the system. Before the first transition, the confined molecules form "pillars" built of 8 molecules, four at each wall, cf. figure 10 (a). We also observe the appearance of "incomplete pillars", but the analysis After this transition, the density is close to 〈ρ〉 ≈ 0.67, the tilting angle is close to θ 1 ≈ 70 • and 〈| sin φ|〉 ≈ 0.67. Also, during the third transition, the changes of 〈cos θ 1 〉 and of 〈| sin φ|〉 are very small. After the
In figures 9 (c) and 9 (d) we show the phase diagrams for the system E0. The first transition appears only at low temperatures and its critical temperature is lower than for the system E1. The critical temperatures of the second and third transitions, however, are slightly higher than for the system E1. Figure 11 summarizes the results obtained for the system E-1 at temperature T * = 0.075. In contrast to the previous cases, only two first-order transitions occur within the investigated range of temperatures. The "survived" transitions correspond to the transitions II and III observed in the systems investigated previously. Indeed, the the first jump on the isotherm is connected with the formation of the stripes.
After the transition, 〈ρ〉 ≈ 0.67 and the translational structure of the system is the same as displayed in

Summary
We have carried out Monte Carlo simulations of Janus-like particles on a cubic lattice, confined in very narrow slits. The particles consisted of two hemispheres that mutually interacted via short-range directional potential. Directional interactions were also assumed between particles and pore walls. We formulated a simple model describing these interactions.
A great majority of the calculations have been performed for bilayers. In this case, we have studied the effects of the competition between fluid-fluid and fluid-pore walls interactions. The principal aim of our calculations was to investigate the phase transitions in the system. The transitions were located by analyzing the histograms of the density. The behavior of the histograms and of other thermodynamic quantities (the heat capacities and compressibilities) suggests the first-order character of the transitions.
However, our calculations were carried out for one (though rather large) size of the system in X and Y directions, and we did not perform the analysis of the dependence of the calculated characteristics on the systems size. To be precise, one should stress that only the analysis of the finite size effects [41,42] could conclusively prove the first-order character of the observed transitions. Nevertheless, we believe that our simulations correctly reveal the number and the range of temperatures, densities and chemical potential of the first-order transitions in the systems under study.
The oriented structures found for the systems in which the fluid-wall interactions are strong and act contrary to the fluid-fluid interactions differ from those appearing in the systems with neutral walls or with the walls attracting the repulsive parts of fluid molecules. In the former case, at low values of the chemical potential, the structures (squares or zigzags) formed at one wall connect the empty sites at the other wall. The extortion of overlaps of the structure at both walls requires a significant increase of the chemical potential. This is particularly visible for the system E5. In the latter case of fluid-wall interactions, the tendency to form overlapping structures is enforced by reinforcement of the orientational effects resulting from fluid-fluid interactions by fluid-wall forces. In all the studied bilayers, the single and double strips of pillars were formed at higher values of the chemical potential.
One can expect that depending on the relations between parameters characterizing Janus-like molecules confined in wider pores, various complicated ordered structures can appear. This opens up new directions in the modelling of modern materials. We hope that the presented results provide a guidance for a further study of concrete systems within the framework of more realistic models.