Different phases of a system of hard rods on three dimensional cubic lattice

TL;DR

This study investigates the phase behavior of hard rods on a cubic lattice, revealing multiple density-driven phase transitions and the emergence of layered phases, with implications for understanding ordering in constrained systems.

Contribution

The paper introduces an efficient cluster algorithm to simulate high-density phases of hard rods and characterizes complex phase transitions, including layered and nematic phases, on a cubic lattice.

Findings

For k ≤ 4, the system remains disordered at all densities.

For k=5,6, a single transition to layered-disordered phase occurs.

For k ≥ 7, multiple transitions including isotropic to nematic phases are observed.

Abstract

We study the different phases of a system of monodispersed hard rods of length $k$ on a cubic lattice using an efficient cluster algorithm which can simulate densities close to the fully-packed limit. For $k\leq 4$, the system is disordered at all densities. For $k=5,6$, we find a single density-driven transition from a disordered phase to high density layered-disordered phase in which the density of rods of one orientation is strongly suppressed, breaking the system into weakly coupled layers. Within a layer, the system is disordered. For $k \geq 7$, three density driven transitions are observed numerically: isotropic to nematic to layered-nematic to layered-disordered. In the layered-nematic phase, the system breaks up into layers, with nematic order in each in each layer, but very weak correlation between the ordering direction between different layers. We argue that the layered-nematic phase is a finite-size effect, and in the thermodynamic limit, the nematic phase will have higher entropy per site.

Equations

$$\mathbf{Q}=|\mathbf{Q}|e^{i\theta }={\rho }_x+{\rho }_y{e}^{\frac{2\pi i}{3}}+{\rho }_z{e}^{\frac{4\pi i}{3}}.$$

$$Q_N$$

$$P_2$$

$${\Omega }_{\{\sigma \}}(\mu ,\mu ,{\mu }^{\prime })={\Omega }_{\{\sigma \}}(\mu ,\mu ,-\infty )\exp \left[A{z}^{\prime }+B{z}^{\prime 2}+...\right].$$

$${\alpha }^r{\beta }^s+{\alpha }^s{\beta }^r\le {\alpha }^{r+s}+{\beta }^{r+s}.$$

$$U_N=1-\frac{⟨|\mathbf{Q}{|}^4⟩}{2⟨|\mathbf{Q}{|}^2{⟩}^2}.$$

Full text

entropy driven, lattice systems, hard rods, nematic

I Introduction

The study of entropy driven phase transitions in systems of long hard rods has a long history dating back to Onsager’s demonstration of a phase transition from an isotropic phase to an orientationally ordered nematic phase with increasing density Onsager (1949). At even higher densities, a system of spherocylinders will show a smectic/columnar phases with partial translational order along with orientational order Polson and Frenkel (1997); Frenkel et al. (1988) and solid-like phases Bolhuis and Frenkel (1997). In two dimensions, there is no true ordered phase, but a system of needles undergoes a Kosterlitz-Thouless type transition into a high density phase with power law correlations Khandkar and Barma (2005); Frenkel and Eppenga (1985); Vink (2009). Physical systems where such transitions are observed include aqueous solutions tobacco mosaic viruses Fraden et al. (1989), liquid crystals De Gennes and Prost (1995), oxygen monolayers adsorbed on molybdenum surfaces Dünweg et al. (1991), carbon nanotube nematic gels Islam et al. (2004) and chlorine atoms adsorbed on silver Mitchell et al. (2001).

The corresponding problem on lattices where the orientation of rods are restricted to the lattice directions (also known as Zwanzig model) has also been studied in parallel. Consider a system of monodispersed rods of length $k$ that occupy $k$ consecutive lattice sites along any one of the lattice directions. Two rods cannot overlap. Early work based on virial expansion Zwanzig (1963), high density expansions Flory (1956) and the Guggenheim approximation DiMarzio (1961), predicted a transition from the low-density disordered phase to a nematic-ordered phase, as the density is increased. However, not much is known about the nature of the high-density phase in this model. At densities near full packing, nematic order is not expected to survive, as there are exponentially many disordered configurations De Gennes and Prost (1995); Ghosh and Dhar (2007). In two dimensions, numerical simulations have shown that the high-density disordered phase has no orientational long-range order, but it is not clear whether some subtle kind of order exists or not. In three dimensions, what is the nature of different phases in the lattice model of hard rods? This is the primary question that we address in this paper. We find that the high-density phase does not have nematic ordering, where rods preferentially align themselves parallel to each other. Instead, for $k>4$, they develop a smectic-like order, where two of the three orientations have large, nearly equal values, and the density of rods of third orientation is very small and the different layers are nearly independent. This way, entropy is maximised, while satisfying the packing constraints at high density. We do obtain a nematic phase, but only at intermediate densities, and only for $k\geq 7$. The determination of the rich structure of different possible orderings in this simple model is the main new finding of this paper. We expect that the behavior in dimensions $d>3$ would be qualitatively similar.

We first summarize known results for the model in two and three dimensions. In two dimensions, for $k=2$ (dimers), it may be shown rigorously that the system is disordered at all densities Heilmann and Lieb (1972), while at full packing, the system is power-law correlated  Fisher and Stephenson (1963); Huse et al. (2003). For $k\leq 6$, Monte Carlo simulations show that the system is disordered at all densities Ghosh and Dhar (2007). For $k\geq 7$, the system undergoes two transitions: first from a low-density disordered phase to an intermediate density nematic phase, and second from the nematic phase to a high density disordered phase Ghosh and Dhar (2007); Kundu et al. (2013); Kundu and Rajesh (2013). While the first transition belongs to the Ising and 3-state Potts universality classes on square Matoz-Fernandez et al. (2008a, b); Fischer and Vink (2009) and triangular lattices Matoz-Fernandez et al. (2008b, c) respectively, the universality class of the second transition has been difficult to resolve Kundu et al. (2013); Kundu and Rajesh (2013). For large enough $k$, the existence of the nematic phase may be rigorously proved Disertori and Giuliani (2013). The model may also be solved exactly on a tree-like lattice, corresponding to a Bethe approximation, and shows a nematic phase for $k\geq 4$, but does not exhibit a second transition Dhar et al. (2011). Density functional theory for rods give a similar result Oettel et al. (2016). The problem of rods may also be generalized to hard rectangles of size $m\times d$ or hard squares of size $m\times m$. For large aspect ratio, the system of rectangles shows four phases with increasing densities: disordered to nematic to columnar to solid-like phases. The detailed phase diagrams for different choices of $m$ and $d$, and the asymptotic behavior of the phase boundaries may be found in Ref. Kundu and Rajesh (2014, 2015a, 2015b); Nath et al. (2015). The $m\times m$ hard square models undergo a single density-driven transition from a disordered phase to a high-density columnar phase Bellemans and Nigam (1967); Ree and Chesnut (1966); Fernandes et al. (2007); Zhitomirsky and Tsunetsugu (2007); Feng et al. (2011); Ramola and Dhar (2012); Ramola et al. (2015); Nath et al. (2016); Mandal et al. (2017); Mandal and Rajesh (2017).

In three dimensions, dimer models ($k=2$) at full packing on bipartite lattices are known to show a Coulomb phase, with algebraic decay of orientational correlations  Huse et al. (2003), while for non-bipartite lattices, the correlations decay faster, and in some exactly solved cases, correlations are strictly zero beyond a finite range  Dhar and Chandra (2008). Not much is known for larger values of $k$. It would be expected that, like in the continuum, there will be an isotropic-nematic transition as the density is increased above zero. However, the minimum value of $k$ is for such a transition to occur is not known. Also, the nature of ordering at still higher densities is not studied much.

In this paper, we study the problem of a monodispersed system of rods of length $k$ using grand canonical Monte Carlo simulations that implements an algorithm with cluster moves. For $k\leq 4$, we find that the system remains disordered and is in the isotropic phase at all densities. When $k=5,6$, we observe a single transition from a low density disordered phase in which the the fractional number of different orientations is nearly equal, to a layered-disordered phase in which the fractional number of one orientation becomes very small, and system develops a layer-like structure, where each layer is a plane with most of the rods being of two orientations lying within the plane, and very weak correlations between different layers. When $k\geq 7$, at intermediate densities, we numerically observe two other phases: a nematic phase, and a new phase that we call the layered-nematic phase. In the layered-nematic phase, each plane has two dimensional nematic order, but there is no overall bulk nematic order. We argue that the existence of this layered-nematic phase is a finite size effect, and that in the thermodynamic limit, parallel orientation between different layers will be regained.

The rest of the paper is organized as follows. In Sec. II, we define the model precisely and describe the grand canonical Monte Carlo scheme that is used to simulate the system. Section III describes the different phases – isotropic, nematic, layered-nematic and layered-disordered – that we observe in our simulations. In Sec. IV, we use perturbation theory to argue that the layered-nematic phase observed in simulations is an artifact of finite system sizes, and the observed behavior should cross over to nematic order for length-scales greater that some crossover scale $L^{*}(\rho)$, where $\rho$ is the density of covered sites. Section V consists results of detailed simulations for systems with $k=2,3,\ldots,7$. The minimum length of rods that is needed for each of the phases to exist is determined. The critical densities and chemical potentials, and other critical parameters are determined for $k=5,6,7$. We end with a summary and discussion of results in Sec. VI.

II Model Description and Monte-Carlo Algorithm

Consider a cubic lattice of size $L\times L\times L$ with periodic boundary conditions. The lattice sites may be occupied by rods that occupy $k$ consecutive lattices sites in any one of the three mutually orthogonal directions. The rods interact only through excluded volume interactions, i.e., a lattice site may be occupied by at most one rod. We associate a weight $e^{\mu}$ with each rod, where $\mu$ is the chemical potential rescaled by temperature. We will call a rod oriented in the $x$-, $y$- and $z$-directions as $x$-mer, $y$-mer, and $z$-mer respectively. The site of a rod with the smallest $x$-, $y$-, and $z$-coordinates will be called its head.

The aim of the paper is to determine the different phases in the hard rod model as the density, or the rod-length $k$ is varied. This is done primarily through grand canonical Monte Carlo simulations. Conventional algorithms with local evaporation and deposition moves fail to equilibrate the system (within available computer time) at large densities because the system gets stuck in long-lived metastable states. Instead, we implement a Monte Carlo algorithm with cluster moves Kundu et al. (2012, 2013) that has recently proved useful in equilibrating systems of hard particles with large excluded volume interactions at densities close to one Kundu et al. (2012, 2013) or at full packing Ramola et al. (2015).

The Monte Carlo algorithm that we use is the following: remove all the $x$-mers, leaving all $y$-mers and $z$-mers undisturbed. The empty intervals in each row in the $x$-direction, separated from each other by $y$-mers or $z$-mers, is now re-occupied by $x$-mers with the correct equilibrium probabilities. The calculation of these probabilities reduces to a one dimensional problem which may be solved exactly (see Refs. Kundu et al. (2013); Ramola et al. (2015); Nath and Rajesh (2014) for details). The evaporation and deposition move satisfies detailed balance as the transition rates depend only on the equilibrium probabilities of the new configuration. Following evaporation and deposition of $x$-mers, we repeat the set of steps with $y$-mers, and then with $z$-mers.

To reduce equilibration and autocorrelation times at high densities, we also implement a flip move. If there is a $k\times k$ square, that is fully covered by $k$ parallel $k$-mers, then we can flip the orientation of $k$-mers, within this square, without affecting any other rods, as shown in the schematic diagram in Fig. 1. Clearly, the flip move does not violate the hard-core constraint and satisfies detailed balance. We define one Monte Carlo time step as updating every row in the $x$-, $y$- and $z$- directions (total of $3L^{2}$ rows), and $L^{3}$(in case of small system sizes) or $L^{3}/k^{2}$(in case of large system sizes) flip moves.

The flip move is crucial for equilibrating the system at densities close to full packing. Figure 2 shows the time evolution of density $\rho$ for a system with $k=7$, starting from nematic initial conditions in which most of the rods lie in the $x$-direction, using the evaporation-deposition algorithm with and without the flip move. The value of $\mu$ is such that the equilibrium configuration does not have nematic order (see Sec. V for details). When the flip move is present, $\rho$ reaches its equilibrium value in about $3\times 10^{5}$ Monte Carlo steps. On the other hand, when the flip move is absent, the system does not reach equilibrium even after $10^{7}$ Monte Carlo steps.

The algorithm is easily parallelized as all the rows can be updated simultaneously. The flip move may also be parallelized by choosing a plane and then choosing one of the $k^{2}$ sublattices randomly. All $k\times k$ squares with left bottom corner lying in this sublattice may be updated simultaneously. All the data presented in the paper is obtained through a parallelized implementation of the algorithm.

III Different phases

In this section, we describe and define the different phases that we observe in our Monte Carlo simulations. Let $\rho_{x}$, $\rho_{y}$ and $\rho_{z}$ be the density of sites occupied by $x$-mers, $y$-mers and $z$-mers respectively. Consider the two-dimensional vector

[TABLE]

We define the nematic order parameters as

[TABLE]

where $\langle\cdots\rangle$ denotes average over the equilibrium probabilities and $\rho=\rho_{x}+\rho_{y}+\rho_{z}$ is the total fraction of sites occupied by $k$-mers.

Isotropic phase: In the isotropic phase, the system is disordered with $\rho_{x}\approx\rho_{y}\approx\rho_{z}$. The probability distribution of $\mathbf{Q}$ is centered about the origin and the order parameters take the value $Q_{N}\approx 0$ and $P_{2}\approx 0$.

Nematic phase: In the nematic phase, a majority of the rods are of one orientation, while the rods of the other two orientations have smaller, roughly equal densities. If $x$- is the preferred direction, then $\rho_{x}\gg\rho_{y}\approx\rho_{z}$, as can be seen in the temporal evolution of three densities shown in Fig. 3(a). A snapshot of a randomly chosen $xy$ plane, as shown in Fig. 3(b), clearly shows that most rods are $x$-mers. In the nematic phase, $Q_{N}\approx\rho$ and $P_{2}\approx 1$.

Layered-Nematic phase: In the layered-nematic phase, there is a spontaneous symmetry breaking, and one of the $xy$, $yz$ or $zx$ planes is selected, and the density of rods that are oriented perpendicular to this plane is suppressed [see Fig. 4(a)], making the system layered. If the chosen plane is the $xy$ plane, then $\rho_{x}\approx\rho_{y}\gg\rho_{z}$. In the layered-nematic phase, within a $xy$-plane, the rods have two-dimensional nematic order. This may seen in the snapshots, shown in Fig. 4(b)–(d), of three randomly chosen $xy$ planes. Each of the planes has two-dimensional nematic order, but could be majority $x$-mers or$y$-mers. To quantify further, we show the time evolution of the local nematic order parameter $n_{x}(z)-n_{y}(z)$, where $n_{x}(z)$ and $n_{y}(z)$ are the densities of sites occupied by $x$-mers and $y$-mers in layer $z$, for four planes in Fig. 4(e). Each of the planes has nematic order which could change sign during the time evolution. There are roughly equal number of planes with majority $x$-mers and majority $y$-mers, as may be seen from the double-peaked probability distribution function $P(n_{x}-n_{y})$, shown in Fig. 4(f), which is obtained by averaging over the different $xy$ planes and over time. If the system has layered-nematic phase, $Q_{N}\approx\rho/2$ and $P_{2}\approx-1$.

Layered-Disordered phase: In the layered-disordered phase, like in the layered-nematic phase, majority of the rods lie in one of the $xy$, $yz$ or $zx$ planes [see Fig. 5(a)]. Let the chosen plane be the $xy$ plane, i.e., $\rho_{x}\approx\rho_{y}\gg\rho_{z}$. In the layered-disordered phase, unlike the layered-nematic phase, the rods within a $xy$ plane do not have nematic order, i.e., $n_{x}(z)\approx n_{y}(z)$ for each layer $z$. This may seen in the snapshots, shown in Fig. 5(b)–(d), of three randomly chosen $xy$ planes, where in each of the planes, there are roughly equal number of $x$-mers and $y$-mers present. The nematic order in each plane fluctuates about zero, as may be seen from the time evolution of the nematic order of four planes as shown in Fig. 5(e)] as well as probability distribution [see Fig. 5(f)] of the local nematic order parameter $n_{x}(z)-n_{y}(z)$. In the layered-disordered phase, $Q_{N}\approx\rho/2$ and $P_{2}\approx-1$.

The order parameters $Q_{N}$ and $P_{2}$ fail to distinguish between the layered-nematic and layered-disordered phases and take the values $Q_{N}\approx\rho/2$ and $P_{2}\approx-1$ for both phases. Though we observe both these phases in our simulations, we argue in the next section that the layered-nematic phase has lower entropy per site that the nematic phase, and is thus a metastable phase.

IV The instability of the layered-nematic phase

In this section, we discuss the instability of the layered-nematic phase. We argue that the layered-nematic phase seen in our simulations is the result of finite size of our samples, and different layers would be expected to develop alignment and hence the usual nematic order if we could study samples of much larger sizes.

We start by considering a system in which the chemical potential of rods in different directions are different: in the $x$- and $y$-directions, it is $\mu$, but in the $z$-direction it is $\mu^{\prime}$. We start with the case when $z^{\prime}=e^{\mu^{\prime}}=0$. In this case, different $z$-layers decouple, and the problem reduces to the problem of $k$-mers on a two dimensional square lattice. We assume that $\mu$ is such that in each layer there is nematic ordering, but in different layers, it may be in different directions. We consider the spontaneous-symmetry broken state $\{\sigma\}$, where the ordering direction in the layer $z=i$ is $\sigma_{i}$, taking values $\pm 1$, depending on the mean orientation being in the $x$- or $y$-directions. There are $2^{L}$ such states, for the $L\times L\times L$ lattice, say with fixed boundary conditions that enforce the specified layered order.

When $z^{\prime}=0$, the states with different $\{\sigma\}$ are degenerate. Now, we develop a perturbation theory for the partition function $\Omega_{\{\sigma\}}(\mu,\mu,\mu^{\prime})$ in powers of $z^{\prime}=\exp(\mu^{\prime})$ Ramola and Dhar (2012); Nath et al. (2015).

[TABLE]

Explicit expressions for the coefficients $A$ and $B$ can be written down in terms of expectation values of appropriate operators in the state $\{\sigma\}$. We find that $A=\epsilon^{k}$, independent of $\{\sigma\}$, where $\epsilon$ is the density of holes in the problem. Hence any difference between different layer-orderings only shows up in the $B$.

Explicit expression for $B$ involves unperturbed probability weight of configurations that have holes that allow two $z$-rods to be put in. Consider two z-mers that have in plane coordinates $(x,y)$ and $(x+\Delta_{1},y+\Delta_{2})$. The probability that we can place these two rods breaks up into products of weights that in the $i$th layer, the sites $(x,y,i)$ and $(x+\Delta_{1},y+\Delta_{2},i)$ are both not occupied by in-plane rods. If this probability is $\alpha(\vec{r})$ or $\beta(\vec{r})$ for $x$ and $y$-orderings in the plane, with $\vec{r}=(\Delta_{1},\Delta_{2})$, the term $B$ is a sum of terms of the form $\alpha^{r}\beta^{s}$, where $r$ is the number of planes with $x$-ordering that intersect both rods, and $s$ is the corresponding number of planes with $y$-ordering. By symmetry in the $x$ and $y$ directions, there will also be a term $\alpha^{s}\beta^{r}$ for the same $\{\sigma\}$ corresponding to separation $\vec{r}^{\prime}=(\Delta_{2},\Delta_{1})$ between the rods. But we notice that, for all $\alpha,\beta\geq 0$, and integers $r,s\geq 0$,

[TABLE]

This implies that the second correction term, when all in-plane nematic orientations are parallel is greater than the term when they are not. Thus, the concentration of $z$-rods induces an effective aligning interaction between nearby layers. Note that this interaction term is proportional to the volume of the system, and would dominate over the degeneracy $2^{L}$ term coming from the number of different states $\{\sigma\}$. This is an order-by-disorder mechanism, where the degeneracy between different equal-weight states $\{\sigma\}$ is lifted, once the perturbation $z^{\prime}$ is introduced.

However, the excess free energy in the ordered nematic state per unit volume is only of ${\mathcal{O}}([\rho^{\prime}/k]^{2})$, where $\rho^{\prime}/k$ is the number density of $z$-mers. In our simulations, for larger values of $\mu$, $\rho^{\prime}$ becomes very small, for instance in $L=112$ and $k=7$, from Fig. 6, we see that $\rho^{\prime}(\mu=5.42)=0.023$ and $\rho^{\prime}(\mu=5.43)=0.0009$, just beyond onset of the layered nematic phase, representing an order of magnitude decrease in $\rho^{\prime}$ as $\mu$ is increased, so that one expects to see configurations with non-parallel nematic order between layers with significant weight if the disordering term $L\ln 2$ is of same order as the the ordering term $L^{3}\rho^{\prime 2}$. Hence, we expect that for $L>L^{*}\sim k/\rho^{\prime}$, the ordering term will win, and nematic ordered state will dominate. However, in our simulations, for $\mu=5.42,k=7$, $\rho^{\prime}\approx 0.023$, and so, $L^{*}$ is of order $300$. For lower values of $\mu$, $\rho^{\prime}$ is larger, and we do see the nematic order. In the other case of $\mu=5.43$, $L^{*}(\rho^{\prime}=0.0009)\sim 7000$. This is much higher than the system size $L=112$, implying that the layered-nematic phase is favored for $\mu=5.43$.

V Phase diagram and critical behavior

We numerically determined the order parameters $Q_{N}$ [see Eq. (2)] and $P_{2}$ [see Eq. (3)] as function of $\rho$ for different $k$ as shown in Fig. 7. We first determine $k_{min}$, the minimum value of $k$ required for each of the phases to appear.

V.1 $k_{min}$

From Fig. 7, it is evident that for $k\leq 4$, both $Q_{N}$ and $P_{2}$ are zero for all values of $\rho$. There are no phase transitions and the system is in the disordered isotropic phase for all densities.

For $k=5$ and $k=6$, $Q_{N}$ increases from [math] to $0.5$ at high densities, while $P_{2}$ simultaneously decreases from [math] to $-1$. These values are indicative of the layered phase, and show that the system undergoes a single transition from an isotropic phase to a layered phase. Thus, for observing a layered phase $k\geq k^{layered}_{min}=5$. We note that there is no nematic phase when $k=5,6$. The critical values for the isotropic-layered transition are: $\mu_{c}(5)\approx 3.82$ and $\rho_{c}(5)\approx 0.874$ and $\mu_{c}(6)\approx 1.0$ and $\rho_{c}(6)\approx 0.68$.

When $k=7$, it may be seen from Fig. 7 that $Q_{N}$ increases from zero to $\approx\rho$ and then decreases to $Q_{N}\approx\rho/2$. Simultaneously, $P_{2}$ increases to $1$ and then drops sharply to $-1$. These values are indicative of nematic and layered phases. We conclude that a nematic phase exists for $k\geq k^{nematic}_{min}=7$. Numerically we find that the layered phase may be further divided into layered-nematic and layered-disordered phases, which is presumably an artifact of the small sizes of our system, as discussed in Sec. IV.

V.2 $k=5,6$

Rods of length $k=5$ are the smallest to show the layered-disordered phase at high densities. We first show that this phase is stable and that the Monte Carlo algorithm equilibrates the system at these densities. To show the stability, we compare the order parameters $Q_{N}$ and $P_{2}$ for two different system size in Fig. 8. The data has only a very weak dependence on the system size, showing that the finite size effects are not important and that the layered phase is stable in the thermodynamic limit. The critical values for the transition is $\mu_{c}(5)\approx 3.82$ and $\rho_{c}(5)\approx 0.874$.

To check that at the high values of $\mu$ and densities $\rho$, our simulations do not suffer from slow down due to jamming problems, we observed the evolution with two different initial conditions: one corresponding to a nematic phase and the other corresponding to an isotropic phase and check that the final state is independent of the initial conditions. The time evolution of $|\mathbf{Q}|$ is shown in Fig. 9 for both of these initial conditions. Clearly, the system loses memory of the initial conditions quite rapidly, and the order parameter reaches a value close to $0.5$, indicative of the layered phase.

We now study of the isotropic to layered-disordered transition. There are three symmetric ordered states. By analogy to the three state Potts model, we expect that this transition should be first order. The numerical data is consistent with a first order transition. First, we show in Fig. 10(a) the probability distribution of the order parameter $Q_{N}$ near the transition point. The distribution has two peaks for values of $\mu$ close to the transition point, one near $Q_{N}\approx 0$, corresponding to the isotropic phase and the other close to $Q_{N}\approx 0.25$, corresponding to the layered phase. Double peaked distribution are a signature of first order transitions and co-existence. This can be further confirmed by looking at two-dimensional density plots of $P(\mathbf{Q})$, as shown in Fig 10(b)–(e), where as $\mu$ is increased, the simultaneous presence of peaks at the origin and and at $\pi/3$, $\pi$ and $5\pi/3$ can be seen in Fig 10(c) and (d).

Further evidence of the first order nature may be obtained by studying the Binder cumulant

[TABLE]

$U_{N}$ is zero in the isotropic phase and $1/2$ in the completely ordered phase. The variation of $U_{N}$ with $\mu$ is shown in Fig. 11. The Binder cumulant becomes negative near the transition point, with its minimum decreasing with system size. Binder cumulant becoming negative is a strong signature of the transition being first order.

The results for $k=6$ are very similar to that for $k=5$. The system undergoes a single transition from isotropic to layered phase with critical parameters are approximately $\mu_{c}(6)\approx 1.0$ and $\rho_{c}(6)\approx 0.68$. We note that the critical values are smaller than that for $k=5$.

V.3 $k=7$

When $k=7$, the systems undergoes a transition from an isotropic phase to a nematic phase at low densities and from nematic phase to a layered phase at high densities (see Fig. 7). Here, we analyze the nature of the transitions as well as the nature of the layered phase. We first discuss the isotropic-nematic transition. There are three symmetric nematic phases corresponding to the three different orientations. By analogy with the three state Potts model, we expect that the transition will be first order in nature. The dependence of the order parameter $Q_{N}$ and the Binder cumulant $U_{N}$ on $\mu$ for different system sizes are shown in Fig. 12. $Q_{N}$ does not show any sign of a discontinuity, nor does the Binder cumulant become negative, both being signatures of a first order transition. Likewise, the probability distribution for $Q_{N}$, shown in Fig. 13 does not show a bimodal distribution for all values of $\mu$ near the critical point. From the crossing of the Binder cumulants, we conclude that the critical parameters are $\mu_{c}\approx-0.23$ corresponding to $\rho_{c}\approx 0.556$. The three state Potts model in 3-dimensions has a very weak first order transition that is difficult to detect in numerical simulations and we expect that the same difficulty holds for the problem of rods. Our simulations do not five a clear evidence of the nature of this phase transition.

We now examine the transition from nematic to layered phase. For $k=7$ and $L=112$, we find a range of $\mu$ ($5.43<\mu<5.60$) for which the system finds itself in the layered-nematic phase. We check that this phase is stable by simulating systems with initial condition that is isotropic, nematic and layered-disordered. Also, we notice that the transition from nematic to layered nematic is accompanied by a is a sharp decrease in $\rho^{\prime}$. However, as we do not expect this to be a thermodynamic phase transition, we did not undertake a detailed study of the layered-nematic phase.

VI Summary and discussion

To summarize, we studied the problem of monodispersed hard rods on a three dimensional cubic lattice using grand canonical Monte Carlo simulations and theoretical methods to obtain the phases for rods of length $k$. We showed that for $k\leq 4$, the system is in a disordered isotropic phase at all densities $\rho$, and there are no phase transitions. For $k=5,6$, the system undergoes a single transition into a high density layered-disordered phase, where the system breaks up into two dimensional layers, but disordered within a layer. For $k=7$, we find that as density is increased, the system makes a transition into a nematic phase. Further increase of density results in a layered-disordered phase. We also observe a layered-nematic phase between the nematic and layered-disordered phases, which we have argued is a finite-size effect.

For values of $k>7$, we expect that the phase diagram remains qualitatively the same as that for $k=7$. We expect the critical density for the isotropic-nematic transition to decrease with increasing $k$, as is confirmed by Monte Carlo simulations of systems with $k=8,9,10$. As seen from Fig. 14, $\rho_{c}$ decreases from 0.556 for $k=7$ to 0.364 for $k=10$. For large $k$, we expect that $\rho_{c}\sim k^{-2}$, as can be seen by assuming that near the transition point, number of rods per $k\times k\times k$ cube should be of of order 1. The nematic-layered transition is essentially a two dimensional transition, as different layers are nearly independent. Thus, we expect that the critical density for this transition varies as $1-a/k^{2}$ for large $k$ as in two dimensions Ghosh and Dhar (2007).

These arguments are easily extended to higher dimensions. For large enough $k$, we will expect a isotropic-nematic transition at a critical density that scales as $k^{-(d-1)}$. A nematic phase may be thought of a union of parallel lines, with hard core constraint along a line, and the problem becomes essentially one dimensional, with weak correlations between different lines. In the high-density phase, we expect that the system will break to two-dimensional layers, with only weak interaction between different layers. This critical density will vary as $1-a/k^{2}$ for large $k$. Preliminary simulations in four dimensions are consistent with the above observations.

From the results of the paper, it is clear that for $k\geq 5$, the fully packed phase shows spontaneous symmetry breaking by selecting the layering plane. It is thus qualitatively different from the $k\leq 4$. Extending the problem of rods to cuboids would result in a much richer phase diagram, as expected from the corresponding case of hard rectangles in two dimensions. However, simulations of such systems is a challenging task.

Acknowledgements.

The simulations were carried out in single node cluster machines Annapurna(2 x Intel Xeon Processor X5570 2.93 GHz) and Nandadevi(2 x Intel Xeon E5-2667 3.3 GHz) using OpenMP C language. We thank Joyjit Kundu and Kedar Damle for helpful discussions. DD’s work is partially supported by the grant DST-SR-S2/JCB-24/2005 of the Government of India.