The nature of the ordered phase of the confined self-assembled rigid rod model

# The nature of the ordered phase of the confined self-assembled rigid rod model

## Abstract

We investigate the nature of the ordered phase and the orientational correlations between adjacent layers of the confined three-dimensional self-assembled rigid rod model, on the cubic lattice. We find that the ordered phase at finite temperatures becomes uniaxial in the thermodynamic limit, by contrast to the ground state (partial) order where the orientation of the uncorrelated layers is perpendicular to one of the three lattice directions. The increase of the orientational correlation between layers as the number of layers increases suggests that the unconfined model may also exhibit uniaxial ordering at finite temperatures.

## I Introduction

State of the art techniques for synthesizing colloids monodisperse in shape and size allow their collective behavior to be investigated Pawar (). The new particles may become the ‘molecules’ of new materials if they can be tailored and assembled into useful structures Glotzer (). In fact, the possibility of particle decoration (through, e.g., glancing angle deposition, templating, or lithography) produces colloids with pre-determined surface patterns (patches). Patches yield new features such as anisotropic interactions, control of the valence, and the formation of permanent electrical dipoles, paving the way for the rational development of novel self-assembled materials (e.g., super-molecules) with highly tunable mechanical, optical, and thermal properties Pawar (); Sciortino ().

Self-assembly has been exploited theoretically for a primitive model of patchy colloids and state of the art simulation studies revealed how the number, type, and distribution of the patches determine the self-assembled structures. In systems with two bonding sites per particle, only (polydisperse) linear chains form and there is no liquid-vapor phase transition Sciortino (). If the linear chains are stiff they will undergo an ordering transition, at fixed concentration, as the temperature decreases. The minimal model of this transition considers the effects of the equilibrium polydispersity and the polymerization process of the rods. In this context, we proposed a model of self-assembled rigid rods (SARR), composed of monomers with two bonding sites that polymerize reversibly into polydisperse chains Tavares2009a () and carried out extensive Monte Carlo simulations to investigate the nature of the ordering transition on the square and triangular lattices Almarza2010 (); Almarza2011 (). The polydisperse rods undergo a continuous ordering transition that was found to be in the two-dimensional (2D) Potts q=2 (Ising) and q=3 universality classes, respectively, as in similar models where the rods are monodisperse Matoz2008b (). These findings refute previous claims, based on Canonical Monte Carlo simulations, that equilibrium polydispersity and the statistical ensemble change the criticality of these models to random percolation Lopez (); Comment (); Reply ().

The nature of the ordering transition of the three-dimensional (3D) SARR model on the simple cubic lattice is much more difficult to establish. The model consists of particles with two patches aligned along , where represents one of the three lattice directions (,,). Particles on nearest-neighbor (NN) lattice sites and interact attractively with energy if their patches are aligned along . Monte Carlo Simulations using efficient algorithms suggest that the ordered phase (below the transition temperature) exhibits a bias towards uniaxial behavior (i.e. the system exhibits a tendency to align different layers, by contrast to the ground state partial order). This tendency is observed only when the system size, defined by with the number of sites considered in the simulation, is sufficiently large at temperatures that are not too low, . Despite the use of efficient cluster algorithms we have not been able to establish the nature of the ordered phase as the system sizes required to observe uniaxial behavior increase rapidly as the temperature decreases, as discussed below.

Here we consider the confined 3D SARR model as a first step towards elucidating the nature of the ordered transition on the cubic lattice. We investigate the nature of the ordered phase and the orientational correlations between adjacent layers of the confined model, on the cubic lattice, and find that the ordered phase at finite temperatures becomes uniaxial in the thermodynamic limit, by contrast to the ground state order where the orientation of the uncorrelated layers is perpendicular to one of the three lattice directions. In addition, we find that the orientational correlation between layers increases as the number of layers increases from two to three suggesting that the unconfined model may also exhibit uniaxial order at finite temperatures.

The paper is arranged as follows: In section II we describe the ground state and the simulation methods used to analyze the 3D SARR model at finite temperatures while in section III we present the simulation results for the ordering transition and the order parameters. In section IV we introduce the confined SARR model. In section V we present the simulation results for the order parameters and the correlations between adjacent layers, for models with two and three layers. We conclude in section VI with a discussion of the results.

## Ii Three-dimensional SARR model

### ii.1 Ground state

In the full lattice limit every site is occupied by one particle aligned in one of the three lattice directions. In the ground state, the SARR model exhibits partial order: One lattice direction (say, ) is suppressed, with the uncoupled layers aligned in one of the remaining lattice directions ( or ). The ground state potential energy is then , where is the number of lattice sites, with degeneracy:

 ωGS=3×(2L−1). (1)

The entropy per site vanishes in the thermodynamic limit.

### ii.2 Simulation procedures

It has been shown that the full-lattice 2D SARR model on the square lattice can be mapped on the 2D Ising model Almarza2010 (). This mapping allows the use of cluster algorithms developed for Potts modelsSwendsen (); Landau-Binder () to enhance the efficiency of the Monte Carlo simulations. In what follows we describe how the Swendsen-Wang algorithm may be adapted to the SARR model on lattices where the mapping does not exist. We recall that the 3D SARR model on the cubic lattice or the 2D SARR model on the triangular lattice cannot be mapped on Ising or Potts models Almarza2011 (). We can, however, develop cluster algorithms based on the layer structure of ground state of the 3D SARR model and the mapping of the 2D SARR model on the square lattice.

#### Cluster sampling

The cluster algorithm samples at each MC step a subset of all sites as described next. One of the three lattice directions is chosen at random, say . Then the sites oriented along are blocked, i.e. their state is frozen during the MC step. Sites with orientations or are active, and their states may change during the MC step. Given the NN character of the Hamiltonian the procedure may be (and it is) applied to all active sites in one MC step. For simplicity, however, we consider one layer, i.e. all sites, , with . The procedure starts by checking the links between pairs of NN active sites (not to be confused with bonds of the original model). Two NN active sites with the same orientation are linked with probability , where . Links cannot be formed between pairs of NN active sites with different orientations. We define clusters of active sites based on the links generated in the previous step. The new configuration is obtained by choosing, independently, for each cluster a new in-plane orientation. The probability of the new cluster orientation, ( and ), is given by:

 Ak(^α)∝exp[−βϵ2nk(^α)]; (2)

where is the number of patches of the cluster that point to blocked sites when the active sites are oriented along .

The procedure is validated using the plaquette formalismAlmarza2010 () that maps the model with blocked sites orientated along to a Potts model Wu () in an external field. For a system with layers () the intralayer potential energy can be written as:

 U=−ϵ∑δ(^αi,^αj)δ(^αi,^rij); (3)

where runs over the NN active sites in one layer, is the orientation of site , and is one if , and zero otherwise. This intralayer Hamiltonian can be mapped to a Potts modelWu () with blocked sites on the square lattice. Using the plaquette formalism Almarza2010 () we find:

 U=−K∑δ(^αi,^αj)−∑

where the subindex runs over pairs of NN sites on the layer with an active site and a passive one. In Eq. (4) is the coupling constant, while and describe the interactions between active sites and the blocked ones (this may be viewed as the interaction of an external field with the active sites). , , and are given in terms of the energy of the patchy model by: , , and . As expected the interaction energy of the active sites pointing to blocked ones is unfavorable. It is now straightforward to implement the cluster algorithm described above. The Potts coupling defines linking criterion between active sites according to the Swendsen-Wang rules Swendsen () while the single-particle interactions are taken into account by considering the effect of an external field Landau-Binder ().

#### Sublattice sampling

In addition to the cluster moves we implemented sublattice single-particle moves. We consider systems with even, and divide the sites into two sublattices: those with odd (sublattice 1) and those with even (sublattice 2). Notice that two NN sites belong to different sublattices. In one sublattice sampling move, we choose one of the sublattices at random and then update the state of each site by computing the interaction with its NNs, , for the three orientations: . The new configuration is obtained by choosing, for each particle, a new orientation with probability: .

In order to check the cluster algorithm and its implementation we have run pairs of simulations using either cluster moves or sublattices moves only. The results were found to be the same within error bars. The cluster algorithm is much more efficient than the sublattice algorithm and the relative efficiency increases as the system size increases. Nevertheless, as we will discuss later, its performance is far from optimal for very large systems and temperatures slightly below the order-disorder transition.

The Monte Carlo simulations of the 3D model are run in cycles. We choose at random, with equal probability, one of the five cycles to be run, namely (two) sublattice and (three) cluster samplings and then proceed as described above.

## Iii Simulation results for the 3D SARR model

The order of the transition can be inferred from the scaling with the system size, of the peak of the excess heat capacity: , where is the potential energy per site. At first order transitions the peak is expected to scale as: Landau-Binder ()

 cmaxv(L)=c0+(ΔU)24kBT2cL3, (5)

where is the transition temperature. In Figure 1 we plot the excess heat capacities in the transition region, and the scaling behavior of their peaks, . The scaling of with indicates that the transition is first order. A least-square fit yields the latent heat of the transition: .

We consider two order parameters to characterize the transition. The first: is based on the partial ground state order:

 O=1−3min[Nx,Ny,Nz]/N, (6)

where is the number of sites with orientation , and is the total number of sites. The second: measures the uniaxial order:

 S=12[3N2x+N2y+N2zN2−1]. (7)

In the ground state (one direction suppressed), while . Values of greater than below the order-disorder transition temperature signal the tendency for uniaxial order. In Figure 2 we plot the two order parameters as a function of the temperature for several system sizes.

Both sets of curves and show, as increases, an abrupt change at (or close to) the temperature where the heat capacity peaks. In small systems, varies monotonically with the temperature. However, for the largest systems peaks at temperatures slightly below the transition temperature, and saturates at values larger than the ground state value at low temperatures. This finding suggests a surprising tendency for uniaxial order, i.e. at finite subcritical temperatures and large system sizes the particles prefer to align in one direction rather than aligning in two directions as expected from the ground state analysis. Notice that as increases the results for below have large error bars. This is a signature of the lost of efficiency of the cluster algorithm to sample slightly below as the system size grows.

With the current algorithms and computational resources, however, we cannot establish the nature of the ordered phase at finite temperatures, in the thermodynamic limit.

## Iv Confined SARR model

In order to investigate the mechanism that may drive uniaxial order in the 3D SARR model, and to quantify it, we have considered a simpler model where the system consists of a number of layers () with sites. These layers are taken perpendicular to direction. Periodic boundary conditions are only considered in directions and . The ground state of the confined SARR model is -degenerate with all the particles in a given layer aligned along the or the direction. We note that geometrical confinement is also used to assist the self-assembly process of 3D systems and thus the study of confinement is also of some practical relevance Kretzschmar ().

The confined systems are simulated using the algorithms described in section II.2, with minor adaptations: namely, the cluster moves are carried out only for layers perpendicular to .

The confined SARR model exhibits an order-disorder transition as the models in 2D and 3D. Note that the limit of confinement (single layer) of the 3D model is not equivalent to the 2D SARR model as the patches can be aligned in three distinct directions. Given the nature of the ground state, where each layer is ordered in an arbitrary direction ( or ) we consider the single-layer order parameters , defined as:

 Si=1L2[Nx(i)−Ny(i)], (8)

where the index refers to one layer. and are the number of sites on layer with patches in the and directions, respectively. Due to the symmetry of the model, the average single-layer order parameters vanish, . In the thermodynamic limit, at low temperatures, we expect an ordering transition described by:

 ⟨S2i⟩{=0;T≥Tc,>0;T

We anticipate a discontinuous transition for systems with a large number of layers (as in 3D), and a continuous one for thin slabs (as in 2D). The continuous transition is expected to be in the 2D Ising class as in single layer systems.

## V Simulation Results for the confined systems

In order to proceed we consider the scaling behavior of the heat capacity, or the related quantity , with the potential energy per site. For systems in the 2D Ising universality class, the scaling behavior at criticality is:

 u′β(βc)∼lnL. (10)

In addition we investigate the scaling behavior of the ratios , related to the Binder cummulants Landau-Binder (). For 2D Ising critical behavior for different cross at the universal value Salas ().

The results for one layer, comply with the expected 2D Ising critical behavior. The critical temperature, , is estimated from the Binder cummulant following standard proceduresAlmarza2010 (). Considering system sizes in the range we find . This result is consistent with the behavior of the pseudocritical temperatures defined by the peaks of the heat capacity as a function of (results not shown). In addition, at the estimated the scaling of the average order parameter exhibits the expected Ising behavior: (where is the critical exponent for the magnetization). The results for and the scaling of are shown in Fig. 3.

For and , the curves for different system sizes () cross at a value of close to that of the Ising universality class (See Fig.4). Somewhat surprisingly, for , the crossing of of the inner layer occurs at a temperature slightly below that of the outer layers (plot not shown). This is likely to be a finite-size effect. A possible explanation is that, for these values of , the different layers are almost independent; within this assumption the pseudo-critical temperature of each layer depends on the density of defects (number of sites with orientation ). The inner layer is expected to have a larger number of defects since these sites oriented along can establish two bonds; by contrast, in the outer layers the sites oriented along can form at most one bond. The larger density of defects reduces the stability of the ordered layer and thus its local pseudo-critical temperature is lower.

At finite temperatures, some particles will be aligned in the direction. An interaction between adjacent layers results from bond formation between particles in different layers (z-bonds) and a correlation between and may appear. If present, these correlations may drive the bias to uniaxial behavior observed in the simulations of the 3D model. Let us define a global order parameter as:

 S=1hh∑i=1Si; (11)

where is the number of layers of the confined model. Again symmetry implies . The average value of may be written as:

 ⟨S2⟩=1h2[h∑i⟨S2i⟩+2h−1∑i=1h∑j=i+1⟨SiSj⟩], (12)

The correlation between two layers is defined as:

 cij=⟨SiSj⟩[⟨S2i⟩⟨S2j⟩]1/2. (13)

The correlation depends both on and , , and it is expected to vanish at low and high temperatures. Inspection of Eq.(12) reveals that behaves in the limit of low temperatures as:

 limT→0⟨S2(L,T)⟩=1h. (14)

We located the order-disorder transition of the model with in the full lattice limit by considering the behavior of . In Figure 4 we plot the results for different system sizes. In the region around the maximum we observe the scaling in line with 2D Ising criticality. In Figure 5 we plot the results for the single-layer, (L,T), and the global, , order parameters for the same model. The single-layer order parameter exhibits the usual dependence on and . The global order parameter, however, exhibits a different behavior: the curves for different system sizes cross around , and then merge as the temperature decreases.

The same qualitative behavior is observed for the confined model with in Figure 6.

While the scaling of and suggests a continuous transition, the crossing of the curves at criticality suggests a (weak) first-order transition, as the order parameter, , in the thermodynamic limit, could exhibit a discontinuity at the transition jumping from zero to a finite value .

In Figure 7a we plot the correlation function between the layers of the system. As expected the correlation decreases and appears to vanish at low and high temperatures. The most relevant feature, however, is that the correlation between layers increases markedly with the system size, . Figure 7b reveals that the correlation increases as (in the range of sizes considered). These results suggest that for , in the thermodynamic limit , the confined model becomes uniaxial (i.e., the layers will align along a unique direction).

Now, we consider the effect of the number of layers on the correlation between adjacent layers. In Figure 8 we plot the layer-layer correlation for and , for systems with . Note that the correlation functions and are equal (except for statistical errors) due to the symmetry of the model. The main conclusion from the results of Fig 8 is that for a fixed value of the correlation between adjacent layers increases with the number of layers. This suggests that the ordered phase of the three dimensional SARR model may become uniaxial, in the thermodynamic limit, at finite temperatures . Note, however, that as the temperature decreases the system size required to observe uniaxial ordering increases very rapidly.

It is clear that the simulation algorithms used in this work loose efficiency as increases at temperatures slightly below the critical temperature. In order to confirm the trend to uniaxiality suggested by the results presented so far, we return to the two layer system, and use an indirect method to compute the free energy difference , where the subscripts indicate configurations where the layers are oriented preferentially in the same () and in different directions. is computed for large using thermodynamic integration from low temperature (where the free energy of the two types of configurations is the same: , as ). The free energy difference is then:

 ΔA(T)T=ΔA(T0)T0+∫TT0ΔU(T′)d1T′ (15)

where . The potential energies are computed using Monte Carlo simulation of relatively large systems, , and , without cluster moves to avoid interconversion between the two types of configurations. The results are plotted in Figure 9. increases with temperature and is proportional to . Thus for large systems and moderate temperatures, the confined SARR model is expected to exhibit uniaxial order.

## Vi Discussion

These surprising results may be interpreted as follows: Consider a two-layer system at low temperature with most particles aligned along the or directions. At , however, a number of particles will align along . A -bond lowers the energy by with respect to two independent -sites, one in each layer, but isolated -bonds do not contribute to the orientational correlation between layers. Now, suppose that two -bonds occur in NN positions (for instance one between sites and , and the second between sites and ) (See Figure 10). At (low) subcritical temperatures this pair of NN -bonds promotes the alignment of both layers in the direction defined by the pair (the direction in this example as shown in Figure 10. In practice, configurations with a different number of pairs of NN bonds along the and directions favor the alignment of the layers.

A bond counting argument gives the probability of aligning one layer along the easy direction over the probability of aligning it along , when a single pair of NN -bonds, along , is present:

 pxpy=eβϵ; (16)

The ratio of the probabilities of aligning the layers over the probability of not doing so is then:

 pxx+pyypxy+pyx=cosh(βϵ). (17)

An estimate of the density of NN -bonds, is:

 α(β)≈e−4βϵ. (18)

Note that only configurations where the number of NN -bonds in the and directions are different contribute to the orientational correlation of the layers. Let us, however, consider the rough estimates given above. At , and . For a system with most configurations will not have NN -bonds, and about one in five (0.186) will have one. An estimate of is then:

 c12≃0.186×cosh(2.5)−1cosh(2.5)+1≃0.134. (19)

which is close to the value obtained from the simulation . This estimate supports the hypothesis that the uniaxial behavior results from the orientational correlation between adjacent layers driven by the presence of NN -bonds.

The characterization of the ordering transition of the confined SARR model with requires the development of more efficient cluster simulation algorithms and thus the behavior of the 3D SARR model cannot be investigated at present. The problem is related, but not identical, to the model for crystallization and vitrification of semiflexible living polymers investigated by Menon and co-workers in 2D and 3D MenonEPL (); MenonPRE ().

###### Acknowledgements.
We acknowledge M. Simões for stimulating discussions in various stages of this work. NGA gratefully acknowledges the support from the Dirección General de Investigación Científica y Técnica under Grant No. FIS2010-15502, and from the Dirección General de Universidades e Investigación de la Comunidad de Madrid under Grant No. S2009/ESP-1691 and Program MODELICO-CM. MMTG and JMT acknowledge financial support from the Portuguese Foundation for Science and Technology (FCT) under Contracts nos. PEst-OE/FIS/UI0618/2011 and PTDC/FIS/098254/2008.

### References

1. A. B. Pawar and I. Kretzschmar, Macromol. Rapid Commun. 31, 150 (2010).
2. S. C. Glotzer and M. J. Solomon, Nature Materials 6, 557-562 (2007).
3. F. Sciortino and E. Zaccarelli, Current Opinion in Solid State and Materials Science 15, 246 - 253 (2011).
4. J. M. Tavares, B. Holder and M. M. Telo da Gama, Phys. Rev E 79, 021505 (2009).
5. N. G. Almarza, J. M. Tavares and M. M. Telo da Gama, Phys. Rev. E 82, 061117 (2010).
6. N. G. Almarza, J. M. Tavares and M. M. Telo da Gama, J. Chem. Phys. 134, 071101 (2011)
7. D. A. Matoz-Fernandez, D. H. Linares, and A. J. Ramirez-Pastor, Europhysics Letters 82, 50007 (2008).
8. L.G. López, D.H. Linares, and A.J. Ramirez-Pastor, Phys. Rev E 80, 040105(R) (2009).
9. L.G. López, D.H. Linares, and A.J. Ramirez-Pastor, Phys. Rev E, 85, 053101 (2012).
10. N. G. Almarza, J. M. Tavares and M. M. Telo da Gama, Phys. Rev. E, 85, 053102 (2012).
11. R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58, 86 (1987).
12. D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, 2nd ed. (Cambridge University Press), Cambridge, 2005.
13. F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
14. I. Kretzschmar and J. H. K. Song, Curr. Opin. Colloid Interface Sci. 16, 84 (2011).
15. J. Salas and. A. D. Sokal, J. Stat. Phys., 98, 551 (2000).
16. G. I. Menon, R. Pandit and M. Barma, Europhysics Letters, 24, 253 (1993).
17. G. I. Menon and R. Pandit, Phys. Rev. E, 59, 787 (1999).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters