Papers in Physics, vol. 1, art. 010005 (2009)

Received: 9 July 2009, Accepted: 6 November 2009
Edited by: J. A. Bertolotto
Licence: Creative Commons Attribution 3.0
DOI: 10.4279/PIP.010005

www.papersinphysics.org

ISSN 1852-4249

The effect of the lateral interactions on the critical behavior of long
straight rigid rods on two-dimensional lattices

P. Longone,1 D. H. Linares,1 A. J. Ramirez-Pastor1∗

Using Monte Carlo simulations and finite-size scaling analysis, the critical behavior of
attractive rigid rods of length k (k-mers) on square lattices has been studied. An ordered
state, with the majority of k-mers being horizontally or vertically aligned, was found. This
ordered phase is separated from the disordered state by a continuous transition occurring at
a critical density θc, which increases linearly with the magnitude of the lateral interactions.

I. Introduction

The study of systems of hard non-spherical col-
loidal particles has, for many years, been attract-
ing a great deal of interest and the activity in this
field is still growing [1–14]. An early seminal con-
tribution to this subject was made by Onsager [1]
with his paper on the isotropic-nematic (I-N) phase
transition in liquid crystals. The Onsager’s theory
predicted that very long and thin rods interacting
with only excluded volume interaction can lead to
long-range orientational (nematic) order. Thus, at
low densities, the molecules are typically far from
each other and the resulting state is an isotropic
gas. However, at large densities, it is more favor-
able for the molecules to align spontaneously (there
are many more ways of placing nearly aligned rods
than randomly oriented ones), and a nematic phase
is present at equilibrium.

Interestingly, a number of papers have appeared
recently, in which the I-N transition was studied
in two dimensions [10–14]. In Ref. [10], the au-

∗E-mail: antorami@unsl.edu.ar

1 Departamento de F́ısica, Instituto de F́ısica Aplicada,
Universidad Nacional de San Luis-CONICET, Cha-
cabuco 917, D5700BWS San Luis, Argentina.

thors gathered strong numerical evidence to sug-
gest that a system of square geometry, with two
allowed orientations, shows nematic order at inter-
mediate densities for k ≥ 7 and provided a qualita-
tive description of a second phase transition (from
a nematic order to a non-nematic state) occurring
at a density close to 1. However, the authors were
not able to determine the critical quantities (crit-
ical point and critical exponents) characterizing
the I-N phase transition occurring in the system.
This problem was resolved in Refs. [11, 12], where
an accurate determination of the critical expo-
nents, along with the behavior of Binder cumulants,
showed that the transition from the low-density dis-
ordered phase to the intermediate-density ordered
phase belongs to the 2D Ising universality class for
square lattices and the three-state Potts univer-
sality class for honeycomb and triangular lattices.
Later, the I-N phase transition was analyzed by
combining Monte Carlo (MC) simulations and the-
oretical analysis [13, 14]. The study in Refs. [13, 14]
allowed (1) to obtain θc as a function of k for
square, triangular and honeycomb lattices, being
θc(k) ∝ k−1 (this dependence was already noted
in Ref. [10]); and (2) to determine the minimum
value of k (kmin), which allows the formation of a
nematic phase on triangular (kmin = 7) and hon-
eycomb (kmin = 11) lattices.

010005-1



Papers in Physics, vol. 1, art. 010005 (2009) / P. Longone et al.

In a recent paper, Fischer and Vink [15] indi-
cated that the transition studied in Refs. [10–14]
corresponds to a liquid-gas transition, rather than
I-N. This interpretation is consistent with the 2D-
Ising critical behavior observed for monodisperse
rigid rods on square lattices [11]. This point will
be discussed in more detail in Sec. III.

In contrast to the systems studied in Refs. [10–
14], many rod-like biological polymers are formed
by monomers reversibly self-assembling into chains
of arbitrary length. Consequently, these systems
exhibit a broad equilibrium distribution of fila-
ment lengths. A model of self-assembled rigid
rods has been recently considered by Tavares et
al. [16]. The authors focused on a system composed
of monomers with two attractive (sticky) poles that
polymerize reversibly into polydisperse chains and,
at the same time, undergo a continuous I-N phase
transition. The obtained results revealed that ne-
matic ordering enhances bonding. In addition, the
average rod length was described quantitatively in
both phases, while the location of the ordering tran-
sition, which was found to be continuous, was pre-
dicted semiquantitatively by the theory.

Beyond the differences between lattice geometry
and the characteristics of the rods (self-assembled
or not), one fundamental feature is preserved in
all the studies mentioned above. This is the as-
sumption that only excluded volume interactions
between the rods are considered (except in Ref. [16],
where monomers with two attractive bonding sites
polymerize into polydisperse rods). Moreover, one
often encounters phrases in the literature, such
as “This theory [Onsager’s theory] shows that re-
pulsive interactions [excluding volume interactions]
alone can lead to long-range orientational nematic
order, disproving the notion that attractive inter-
actions are a prerequisite” [17], which could be am-
biguous with respect to the role that attractive lat-
eral interactions between the rods should play in
reinforcing (or not) the nematic order.

In this context, it is of interest and of value to
inquire how the existence of lateral interactions be-
tween the rods influences the phase transition oc-
curring in the system. The objective of this paper
is to provide a thorough analysis in this direction.
For this purpose, an exhaustive study of the phase
transition occurring in a system of attractive rigid
rods deposited on square lattices was performed.
The results revealed that (i) the orientational or-

der survives in the presence of attractive lateral in-
teractions; (ii) the critical density shifts to higher
values as the magnitude of the lateral interactions
is increased; and (iii) the continuous transition be-
comes first order for interaction strength w > wc
(in absolute values).

The outline of the paper is as follows. In Sec.
II we describe the lattice-gas model and the sim-
ulation scheme. In Sec. III we present the MC
results. Finally, the general conclusions are given
in Sec. IV.

II. Lattice-gas model and Monte
Carlo simulation scheme

We address the general case of adsorbates assumed
to be linear rigid particles containing k identical
units (k-mers), with each one occupying a lattice
site. Small adsorbates would correspond to the
monomer limit (k = 1). The distance between k-
mer units is assumed to be equal to the lattice con-
stant; hence exactly k sites are occupied by a k-mer
when adsorbed (see Fig. 1). The surface is repre-
sented as an array of M = L × L adsorptive sites
in a square lattice arrangement, where L denotes
the linear size of the array. In order to describe
the system of N k-mers adsorbed on M sites at a
given temperature T , let us introduce the occupa-
tion variable ci which can take the values ci = 0
if the corresponding site is empty and ci = 1 if
the site is occupied. On the other hand, molecules
adsorb or desorb as one unit, neglecting any pos-
sible dissociation. Under these considerations, the
Hamiltonian of the system is given by

H = w
∑
〈i,j〉

cicj − N (k − 1)w + �o
∑

i

ci (1)

where w is the nearest-neighbor (NN) interaction
constant which is assumed to be attractive (nega-
tive), 〈i, j〉 represents pairs of NN sites and �o is
the energy of adsorption of one given surface site.
The term N (k − 1)w is subtracted in eq. (1) since
the summation over all the pairs of NN sites over-
estimates the total energy by including N (k − 1)
bonds belonging to the N adsorbed k-mers. Be-
cause the surface was assumed to be homogeneous,
the interaction energy between the adsorbed dimer
and the atoms of the substrate �o was neglected for
the sake of simplicity.

010005-2



Papers in Physics, vol. 1, art. 010005 (2009) / P. Longone et al.

Figure 1: Linear tetramers adsorbed on square lat-
tices. Full and empty circles represent tetramer
units and empty sites, respectively.

In order to characterize the phase transition, we
use the order parameter defined in Ref. [11], which
in this case can be written as

δ =
|nh − nv|
nh + nv

(2)

where nh(nv) is the number of rods aligned along
the horizontal (vertical) direction. When the sys-
tem is disordered (θ < θc), all orientations are
equivalents and δ is zero. As the density is in-
creased above θc, the k-mers align along one direc-
tion and δ is different from zero. Thus, δ appears
as a proper order parameter to elucidate the phase
transition.

The problem has been studied by grand canon-
ical MC simulations using a typical adsorption-
desorption algorithm. The procedure is as fol-
lows. Once the value of the chemical potential
µ is set, a linear k-uple of nearest-neighbor sites
is chosen at random and an attempt is made to
change its occupancy state with probability W =
min {1, exp (−∆H/kBT )}, where ∆H = Hf −Hi is
the difference between the Hamiltonians of the final
and initial states and kB is the Boltzmann constant.
In addition, displacement (diffusional relaxation) of
adparticles to nearest-neighbor positions, by either
jumps along the k-mer axis or reptation by rotation
around the k-mer end, must be allowed in order to
reach equilibrium in a reasonable time. A MC step

Figure 2: Adsorption isotherms (coverage versus
chemical potential) for k = 10, L = 100 different
w/kBT ’s as indicated. Inset: Adsorption phase di-
agram of attractive 10-mers on square lattices.

(MCs) is achieved when M k-uples of sites have
been tested to change its occupancy state. Typi-
cally, the equilibrium state can be well reproduced
after discarding the first r′ = 107 MCs. Then, the
next r = 2×107 MCs are used to compute averages.

In our MC simulations, we varied the chemical
potential µ and monitored the density θ and the
order parameter δ, which can be calculated as sim-
ple averages. The reduced fourth-order cumulant
UL introduced by Binder [18] was calculated as:

UL = 1 −
〈δ4〉

3〈δ2〉2
, (3)

where 〈· · ·〉 means the average over the MC simu-
lation runs. All calculations were carried out using
the BACO parallel cluster (composed by 60 PCs
each with a 3.0 GHz Pentium-4 processor and 90
PCs each with a 2.4 GHz Core 2 Quad processor)
located at Instituto de F́ısica Aplicada, Universi-
dad Nacional de San Luis-CONICET, San Luis,
Argentina.

III. Results

The calculations were developed for linear 10-mers
(k = 10). With this value of k and for non-

010005-3



Papers in Physics, vol. 1, art. 010005 (2009) / P. Longone et al.

interacting rods, it is expected the existence of a
nematic phase at intermediate densities [10]. The
surface was represented as an array of adsorptive
sites in a square lattice arrangement with conven-
tional periodic boundary conditions. The effect of
finite size was investigated by examining lattices
with L = 50, 100, 150, 200.

In order to understand the basic phenomenology,
we consider, in the first place, the behavior of the
adsorption isotherms in presence of attractive lat-
eral interactions between the k-mers.

Fig. 2 shows typical adsorption isotherms (cover-
age versus µ/kBT ) for linear 10-mers with different
values of the lateral interaction (the solid circles
represent the Langmuir case, w/kBT = 0).

The isotherms shift to lower values of chemical
potential, and their slopes increase as the ratio
w/kBT increases (in absolute value). For inter-
action strength above the critical value (w > wc,
in absolute values) the system undergoes a first-
order phase transition, which is observed in the
clear discontinuity in the adsorption isotherms1. In
the case studied, this critical value is approximately
wc/kBT ≈ −0.80 (or kBTc/w ≈ −1.25). The be-
havior of the adsorption isotherms also allows us to
calculate the phase diagram of the adsorbed mono-
layer in “temperature-coverage” coordinates. In
fact, once obtained the real value of the chemical
potential (or critical chemical potential µc) in the
two-phase region, the corresponding critical densi-
ties can be easily calculated. By repeating this pro-
cedure for different temperatures ranging between
0 and Tc, the coexistence curve can be built [20].
A typical phase diagram, obtained in this case for
attractive 10-mers, is shown in the inset of Fig. 2.

On the basis of the study in Fig. 2, our next
objective is to obtain evidence for the existence of
nematic order in the range −0.80 ≤ w/kBT < 0 of
attractive interactions. For this purpose, the be-
havior of the order parameter δ as a function of
coverage was analyzed for k = 10, L = 100 and dif-

1In this situation, which has been observed experimen-
tally in numerous systems, the only phase which one ex-
pects is a lattice-gas phase at low coverage, separated by a
two-phase coexistence region from a “lattice-fluid” phase at
higher coverage. This condensation of a two-dimensional gas
to a two-dimensional liquid is similar to that of a lattice-gas
of attractive monomers. However, the symmetry particle-
vacancy (valid for monoatomic particles) is broken for k-
mers and the isotherms are asymmetric with respect to
θ = 0.5.

Figure 3: Surface coverage dependence of the ne-
matic order parameter for k = 10, L = 100 different
w/kBT ’s as indicated.

ferent values of the lateral interaction. The results
are shown in Fig. 3, revealing that (i) the orien-
tational order survives in the presence of attrac-
tive lateral interactions and (ii) the critical density
shifts to higher values as the magnitude of the lat-
eral interactions is increased.

In order to corroborate the results obtained in
the last figure, we now study the dependence of
θc on w/kBT . In the case of the standard the-
ory of FSS [18, 19], when the phase transition is
temperature driven, the technique allows for vari-
ous efficient routes to estimate Tc from MC data.
One of these methods, which will be used in this
case, is from the temperature dependence of UL(T ),
which is independent of the system size for T = Tc.
In other words, Tc is found from the intersection
of the curve UL(T ) for different values of L, since
UL(Tc) =const. In our study, we modified the con-
ventional FSS analysis by replacing temperature by
density [11]. Under this condition, the critical den-
sity has been estimated from the plots of the re-
duced four-order cumulants UL(θ) plotted versus
θ for several lattice sizes. As an example, Fig. 4
shows the results for w/kBT = −0.125. In this
case, the value obtained was θc = 0.542(2). In
the inset, the data are plotted over a wider range
of temperatures, exhibiting the typical behavior of
the cumulants in the presence of a continuous phase

010005-4



Papers in Physics, vol. 1, art. 010005 (2009) / P. Longone et al.

Figure 4: Curves of UL(θ) vs θ for k = 10,
w/kBT = −0.125 and square lattices of different
sizes. From their intersections one obtained θc. In
the inset, the data are plotted over a wider range
of densities.

transition.
The procedure of Fig. 4 was repeated for

−0.80 ≤ w/kBT < 0, showing that the values of
θc increase linearly with the magnitude of the lat-
eral couplings (see solid squares in Fig. 5). The
critical line (dotted line in the figure) was obtained
from the linear fit of the numerical data. As it is
possible to observe, the range of coverage at which
the transition occurs diminishes as w/kBT is in-
creased (in absolute value). This finding indicates
that the presence of attractive lateral interactions
between the rods does not favor the formation of
nematic order in the adlayer. The phenomenon
can be understood from the behavior of the sec-
ond virial coefficient, which will initially decrease
on introducing attractive w. This decrease implies
that the isotherms shift to lower values of chemical
potential, and consequently, the critical point shifts
to higher densities.

We did not assume any particular universality
class for the transitions analyzed here in order to
calculate their critical densities, since the analy-
sis relied on the order parameter cumulant’s prop-
erties. However, the fixed value of the cumu-
lants, U ∗ = 0.617(9), is consistent with the ex-
tremely precise transfer matrix calculation of U ∗ =

Figure 5: Temperature-coverage phase diagram
corresponding to attractive k-mers with k =
10. The inset in the upper-left (lower-right) cor-
ner shows a typical configuration in the nematic
(isotropic) phase.

0.6106901(5) [21] for the 2D Ising model. This find-
ing may be taken as an indication that the phase
transition belongs to the 2D Ising universality class.

With respect to the behavior of the system for
w/kBT < −0.80, the adsorbed layer “jumps” from
a low-coverage phase to a high-coverage phase.
This effect, which has been discussed in Fig. 2,
is represented in Fig. 5 by the dashed coexistence
line. The low-coverage phase is an isotropic state,
similar to that observed for w/kBT > −0.80 and
low density (see inset in the lower-right corner of
Fig. 5). On the other hand, the high-coverage
phase is also an isotropic state, but characterized by
the presence of local orientational order (domains
of parallel k-mers). A typical configuration in this
regime is shown in Fig. 6).

Finally, it is worth pointing out that: (1) the be-
havior of the order parameter in Fig. 3 clearly indi-
cates that the transition from the low-density dis-
ordered phase to the intermediate-density ordered
phase is an isotropic to nematic phase transition
(when all the words have the usual meaning). In
this case, the transition under study belongs to the
2D Ising universality class. It can also be thought
of as an unmixing or liquid-gas transition [15]. For
this reason we have called gas and liquid to the

010005-5



Papers in Physics, vol. 1, art. 010005 (2009) / P. Longone et al.

Figure 6: Typical configuration of the adlayer in
the high-coverage phase and w/kBT < −0.80.

phases reported in Fig. 5; and (2) even though it
has not been rigorously proved yet, a second phase
transition for non-interacting rods at high densi-
ties has been theoretically predicted [10] and nu-
merically confirmed [13]. This result has not been
confirmed for the case of attractive rods. An ex-
haustive study on this subject will be the object of
future work.

IV. Conclusions

We have addressed the critical properties of attrac-
tive rigid rods on square lattices with two allowed
orientations, and shown the dependence of the crit-
ical density on the magnitude of the lateral interac-
tions w/kBT . The results were obtained by using
MC simulations and FSS theory.

Several conclusions can be drawn from the
present work. On the one hand, we found that
even though the presence of attractive lateral in-
teractions between the rods does not favor the for-
mation of nematic order in the adlayer, the orien-
tational order survives in a range that goes from
w/kBT = 0 up to wc/kBT ≈ −0.80 (wc/kBT rep-
resents the critical value at which occurs a typical
transition of condensation in the adlayer). In this
region of w/kBT , the critical density increases lin-
early with the magnitude of the lateral couplings.
On the other hand, the evaluation of the fixed point
value of the cumulants U ∗ = 0.617(9) indicates

that, as in the case of non-interacting rods, the ob-
served phase transition belongs to the universality
class of the two-dimensional Ising model.

With respect to the behavior of the system for
w/kBT < −0.80, the continuous transition be-
comes first order. Thus, the adsorbed layer jumps
from a low-coverage phase, similar to that observed
for w/kBT > −0.80 and low density, to an isotropic
phase at high coverage, characterized by the pres-
ence of local orientational order (domains of paral-
lel k-mers)

Future efforts will be directed to (1) extend the
study to repulsive lateral interaction between the
k-mers; (2) obtain the whole phase diagram in
the space (temperature-coverage-rod’s size); (3) de-
velop an exhaustive study on critical exponents and
universality and (4) characterize the second phase
transition from a nematic order to a non-nematic
state occurring at high density.

Acknowledgements - This work was supported in
part by CONICET (Argentina) under project num-
ber PIP 112-200801-01332; Universidad Nacional
de San Luis (Argentina) under project 322000 and
the National Agency of Scientific and Technological
Promotion (Argentina) under project 33328 PICT
2005.

[1] L Onsager, The effects of shape on the inter-
action of colloidal particles, Ann. N. Y. Acad.
Sci. 51, 627 (1949).

[2] P J Flory, Thermodynamics of high polymer
solutions, J. Chem. Phys. 10, 51 (1942); P J
Flory, Principles of Polymers Chemistry, Cor-
nell University Press, Ithaca, NY (1953).

[3] M L Huggins, Some properties of solutions of
long-chain compounds, J. Phys. Chem. 46, 151
(1942); M L Huggins, Thermodynamic prop-
erties of solutions of long-chain compounds,
Ann. N. Y. Acad. Sci. 43, 1 (1942); M L Hug-
gins, Theory of solutions of high polymers, J.
Am. Chem. Soc. 64, 1712 (1942).

[4] J P Straley, Liquid crystals in two dimensions,
Phys. Rev. A 4, 675 (1971).

010005-6



Papers in Physics, vol. 1, art. 010005 (2009) / P. Longone et al.

[5] J Vieillard-Baron, Phase transitions of the
classical hard-ellipse system, J. Chem. Phys.
56, 4729 (1972).

[6] D Frenkel, R Eppenga, Evidence for algebraic
orientational order in a two-dimensional hard-
core nematic, Phys. Rev. A 31, 1776 (1985).

[7] K J Strandburg, Two-dimensional melting,
Rev. Mod. Phys. 60, 161 (1988).

[8] A J Phares, F J Wunderlich, Thermodynam-
ics and molecular freedom of dimers on plane
triangular lattices, J. Math. Phys. 27, 1099
(1986).

[9] A J Phares, F J Wunderlich, J D Curley, D
W Grumbine Jr, Structural ordering of inter-
acting dimers on a square lattice, J. Phys. A:
Math. Gen. 26, 6847 (1993).

[10] A Ghosh, D Dhar, On the orientational order-
ing of long rods on a lattice, Eur. Phys. Lett.
78, 20003 (2007).

[11] D A Matoz-Fernandez, D H Linares, A J
Ramirez-Pastor, Determination of the criti-
cal exponents for the isotropic-nematic phase
transition in a system of long rods on two-
dimensional lattices: Universality of the tran-
sition, Europhys. Lett. 82, 50007 (2008).

[12] D A Matoz-Fernandez, D H Linares, A J
Ramirez-Pastor, Critical behavior of long lin-
ear k-mers on honeycomb lattices, Physica A
387, 6513 (2008).

[13] D H Linares, F Romá, A J Ramirez-Pastor,
Entropy-driven phase transition in a system of
long rods on a square lattice, J. Stat. Mech.
P03013 (2008).

[14] D A Matoz-Fernandez, D H Linares, A
J Ramirez-Pastor, Critical behavior of long
straight rigid rods on two-dimensional lat-
tices: Theory and Monte Carlo simulations,
J. Chem. Phys. 128, 214902 (2008).

[15] T Fischer, R L C Vink, Restricted orientation
“liquid crystal” in two dimensions: Isotropic-
nematic transition or liquid-gas one(?), Euro-
phys. Lett. 85, 56003 (2009).

[16] J M Tavares, B Holder, M M Telo da Gama,
Structure and phase diagram of self-assembled
rigid rods: Equilibrium polydispersity and ne-
matic ordering in two dimensions, Phys. Rev.
E 79, 021505 (2009).

[17] H H Wensink, Columnar versus smectic order
in systems of charged colloidal rods, J. Chem.
Phys. 126, 194901 (2007).

[18] K Binder, Applications of the Monte Carlo
Method in Statistical Physics. Topics in cur-
rent Physics, Springer, Berlin (1984).

[19] V Privman, Finite Size Scaling and Numer-
ical Simulation of Statistical Systems, World
Scientific, Singapore (1990).

[20] T L Hill, An Introduction to Statistical
Thermodynamics, Addison Wesley Publishing
Company, Reading, MA (1960).

[21] G Kamieniarz, H W J Blöte, Universal ratio
of magnetization moments in two-dimensional
Ising models, J. Phys. A: Math. Gen. 26, 201
(1993).

010005-7