The isotropic-nematic transition for hard rods on a three--dimensional (3D) cubic lattice

TL;DR

This study uses Monte Carlo simulations to analyze the isotropic-nematic transition of hard rods on a 3D lattice, revealing a weakly first-order transition with strong fluctuations and comparing results to theoretical predictions.

Contribution

First detailed Monte Carlo analysis of the isotropic-nematic transition for hard rods on a 3D lattice, highlighting the transition's weakly first-order nature and discrepancies with theory.

Findings

Transition occurs for rod length L >= 6.

Weakly first-order transition with no clear density jump.

Strong fluctuations observed near the transition.

Abstract

Using grand canonical Monte Carlo (GCMC) simulations, we investigate the isotropic-nematic phase transition for hard rods of size Lx1x1 on a 3D cubic lattice. We observe such a transition for L >= 6. For L = 6, the nematic state has a negative order parameter, reflecting the co--occurrence of two dominating orientations. For L >= 7, the nematic state has a positive order parameter, corresponding to the dominance of one orientation. We investigate rod lengths up to L = 25 and find evidence for a very weakly first-order isotropic-nematic transition, while we cannot completely rule out a second order transition. It was not possible to detect a density jump at the transition, despite using large systems containing several 10^5 particles. The probability density distributions P(Q) from the GCMC simulations near the transition are very broad, pointing to strong fluctuations. Our results complement earlier results on the demixing (pseudonematic) transition for an equivalent system in 2D, which is presumably of Ising--type and is present for L >= 7. We compare our results to lattice fundamental measure theory (FMT) and find that FMT strongly overestimates nematic order and consequently predicts a strong first order transition. The rod packing fraction of the nematic coexisting states, however, agree reasonably well between FMT and GCMC.

Figures

Equations

$$Q_i$$

$$S_i$$

$$\Xi (z_i)=\sum _{N_1=0}^{\infty }\sum _{N_2=0}^{\infty }\sum _{N_3=0}^{\infty }\frac{z_1^{N_1}{z}_2^{N_2}{z}_3^{N_3}}{N_1!N_2!N_3!}\mathcal{N}(\{N_j\}),$$

$$\Xi (z)=\sum _{N=0}^{\infty }\frac{z^N}{N!}{\mathcal{N}}^{\prime }(N),$$

$$⟨A⟩=\frac{1}{\Xi }\sum _{N=0}^{\infty }\frac{z^N}{N!}\sum _{\mathcal{S}}A(\{N_j\}).$$

$$f$$

$$\beta f^{\mathrm{id}}$$

$$\beta f^{\mathrm{ex}}$$

$${\Phi }^{0D}(\eta )$$

$$\displaystyle{\color[rgb]{0,0,0}\tilde{Q}_{i}}$$

$$\displaystyle{\color[rgb]{0,0,0}\tilde{S}_{i}}$$

$$g_M=2-\frac{⟨m_i^4⟩}{⟨m_i^2{⟩}^2},$$

Full text

I Introduction

Investigating phase transitions on lattice models has often helped formulate general concepts in phase transition theory, if not initiate them. Additionally, lattice models often take up the role of a simplified version of a continuum model of interest, in the hope they retain the basic physics, e.g. the type of phase transition. Popular examples are using the lattice gas to explore the liquid–vapor transition or lattice polymer models in order to reduce the large configurational space accompanying their continuum progenitors, enabling specific questions to be addressed.

Systems of hard, anisotropic particles have found steady interest over the past decades since they exhibit many types of phase transitions seen in molecular or colloidal materials, while being purely entropic in nature: These include transitions from an isotropic liquid phase to various smectic and crystalline phases via nematic and biaxial states Mar14 . Not only can hard–body models in the continuum be treated efficiently in the context of simulations All93 , but also advanced density functional methods can describe them accurately (e.g. fundamental measure theory) Mar14 ; Witt_thesis15 . It is hence interesting to see whether simple lattice analogs of these models show equivalent phase transitions, and likewise are accessible to density functional methods. In two dimensions (2D), hard rods with side–lengths $L\times 1$ on a square lattice define the most basic lattice model; in 3D, the corresponding lattice model entails hard rods with side–lengths $L\times 1\times 1$ on a cubic lattice.

In 2D, several simulation studies have investigated the basic model introduced above. One can view it as a binary system of particles (rods oriented in $x$– and $y$–direction, respectively), thus a transition to a bulk state with preferred orientation corresponds to a demixing transition presumably of Ising type. This demixing is found for $L\geq 7$ in simulations Ghosh07 . The critical packing fraction for the onset of demixing scales approximately as $4.8/L$ for large $L$ Mat08 ; Kun16 . At very high packing fractions $\eta\approx 1$, theoretical arguments predict a reentrant transition from the demixed to a disordered state bearing some characteristics of a cubatic phase on a lattice Ghosh07 . This transition has been studied in more detail using simulations in Refs. Kun13 ; Kun14 . For rods with extensions $mL\times m$ (where $mL,m$ are integer and $L$ may be noninteger), the phase diagram has been investigated in Refs. Kun15 ; Kun16 where (for $m>1$) it is shown that a columnar phase appears between the demixed and high-density disordered phase.

Theoretical studies of the bulk phase of the basic models in 2D and 3D were sparked by DiMarzio DiMar61 , having approached the problem from the context of polymer theory. DiMarzio calculated the number of possible packings of rods—thus evaluated the entropy—through approximating the probability of inserting a new rod into a system already containing other rods. The approximation consists of treating all rods orthogonal to the rod that is to be inserted in a typical mean–field fashion: one such orthogonal rods (of length $L$) is considered a ‘cloud’ of $L$ independently–distributed obstacles (of size $1\times 1$ in 2D and $1\times 1\times 1$ in 3D) for the rod to be inserted. DiMarzio, however, did not evaluate the bulk state explicitly by maximizing the entropy. This appears to have first been done by Alben Alb71 for the 3D model, who found a strong first order nematic transition for rod lengths $L\geq 4$. The equivalence of the DiMarzio entropy to exact solutions on Bethe–like lattices has been investigated in Ref. Dhar11 (see also the much earlier Ref. Cot69 regarding such a connection). Among other results, Alben’s solution is recovered, and the 2D solution is given: demixing of the two rod species occurs for $L\geq 4$ above a critical packing fraction $\eta_{\rm c}(L)=2/(L-1)$. This is in contrast to the aforementioned simulations Ghosh07 , which found demixing for $L\geq 7$.

A completely different theoretical route for treating hard particles on lattices has been initiated by Lafuente and Cuesta Laf02 ; Laf04 using Fundamental Measure Theory (FMT). They provide a general recipe for how to construct a density functional for arbitrary hard particle mixtures, applicable to inhomogeneous situations. The construction of the functional begins with the requirement that the functional is exact for extreme confinement in so–called zero–dimensional cavities. A zero–dimensional cavity is a restricted domain on the lattice which can only hold one particle at a time and whose free energy is known. Such a cavity may consist of more than one point where the hard particle is positioned. For the basic rod models in 2D and 3D, the FMT functional and its bulk properties have been worked out in Ref. Oet16 . Surprisingly, it turns out that the bulk entropy is exactly the same as DiMarzio’s (which was not mentioned in Oet16 ), as well as that for rods on Bethe-like lattices with coordination numbers 4 and 6 (representing 2D and 3D, respectively) Dhar11 ; therefore, FMT likewise predicts a strong first order transition from an isotropic to a nematic state. The relative coexistence gap $2(\rho_{\rm nem}-\rho_{\rm iso})/(\rho_{\rm nem}+\rho_{\rm iso})$ ($\rho_{\rm nem}$ and $\rho_{\rm iso}$ are the coexistence number densities per lattice site in the nematic and isotropic state, respectively) is narrow for shorter $L$ (e.g., 2.5% for $L=4$), but quickly widens, becoming about 40% in the limit of very long rods. This is in line with the general expectations from Landau theory, which predicts (by means of general symmetry arguments) a first order liquid–nematic transition in 3D for hard bodies vanRoij16 .

It is interesting, therefore, to see whether simulations of the basic rod model in 3D can confirm this type of transition. This is the purpose of the present work. To our surprise, we have not found studies in the literature addressing the nematic transition for this basic 3D rod model. In Ref. Case95 , rods with size $15\times 3\times 3$ were studied and a layering transition was found. As shown below, our grand canonical simulations for rod lengths $L$ between 5 and 25 indicate that the isotropic–nematic phase transition sets in for $L=6$ (resembling the layering transition of Ref. Case95 ), moreover changes its character for $L\geq 7$ (explained below), and, rather peculiarly, is possibly continuous or only of very weak first–order type up to $L=25$, contrary to the theory above. The remainder of the paper is structured as follows: In Sec. II, we introduce the model and its order parameters. In Sec. III, the application of grand canonical Monte Carlo simulations for the model is described, and the FMT bulk free energy is recapitulated. The simulation results are described in Sec. IV and compared to FMT. A summary and a discussion of how the present results relate to those for continuum models is given in Sec. V. In the Appendix, the qualitative analogy of the 3D hard rod model with the three–state Potts model is exploited, and an investigation with an appropriate order parameter is done for $L=8$, including finite–size analysis.

II Model

We consider a simple cubic lattice in 3D, the unit cell length is set to 1. Hard rods are parallelepipeds with extensions $L\times 1\times 1$ and corners sitting on lattice points. The position of a rod is specified by the corner whose lattice coordinates are minimal each. They are allowed to touch (i.e., share corners or faces), but forbidden to overlap. The cubic lattice restricts possible orientations to three, and we refer to rods oriented in $x$– resp. $y$– resp. $z$–direction as species 1 resp. 2 resp. 3. Species densities $\rho_{i}$ are defined as number of rods of species $i$ ($N_{i}$) per lattice site, $\rho=\rho_{1}+\rho_{2}+\rho_{3}$ is the total density and $\eta=\rho L\leq 1$ is the total packing fraction. We introduce bulk nematic order parameters $Q_{i}$ and biaxiality order parameters $S_{i}$ by

[TABLE]

where $(ijk)$ is a cyclic permutation of $(123)$.

An allowed configuration of the system is given by the set of all rod positions ${\mathcal{S}}=\{\mathbf{s}_{j,k}\}$ where $j=1,2,3$ indicates the species and $k=1...N_{j}$ the rod number and in which the rods do not overlap. The total number of allowed configuations is denoted by ${\mathcal{N}}(\{N_{j}\})$. The grand partition function of the system is defined by

[TABLE]

where $z_{i}=\exp(\beta\mu_{i})$ is the activity of species $i$ ($\beta=1/(k_{\rm B}T)$ is the inverse temperature and $\mu_{i}$ is the chemical potential of species $i$). In a bulk system all $z_{i}$ are equal, $z_{i}=z$. Note that the 3 possible particle orientations can be treated as an internal property of particles such that the grand partition function in the bulk can be written alternatively as a single sum over the total number of particles:

[TABLE]

where ${\mathcal{N}^{\prime}}(N)$ is the total number of allowed configurations having $N$ particles with three possible orientations each. The grand canonical average of an observable $A(\{N_{j}\})$ is defined by

[TABLE]

III Methods

III.1 Grand canonical Monte Carlo (GCMC) simulations

We have simulated the model on cubic lattices with size $V=M^{3}$. We have worked with lattice sizes between $M=50$ and $M=170$, depending on rod length. In each step, a particle insertion attempt or deletion attempt was chosen with probability $1/2$. For an insertion move, the orientation and the position of the particle was picked randomly. For the deletion move, one of the particles in the system was chosen randomly. Detailed balance was obeyed: the insertion move was accepted with a probability ${\rm min}(1,3zV/(N+1))$ if it lead to an allowed configuration, and the deletion move was accepted with a probability ${\rm min}(1,N/(3zV))$. $N$ indicates the number of particles before the insertion/deletion move was attempted.

To compute densities and orientation variables, we stored various “time” series of the rod numbers $N_{i}(t)$, measured every $10^{7}$ moves (otherwise noted). In isotropic states, $N_{i}(t)$ fluctuated around $\langle N\rangle/3$, and equilibrium was reached promptly. On the other hand, one species having a pronounced majority was a typical sign of being in e.g. the nematic state; there, intervals of alternating majority species were observed in the time series. We assumed full equilibration of the system if of the order of 10 such intervals occurred. This criterion was of course difficult to fulfill deep in the nematic phase. The collected time series were used to calculate averages $\langle A\rangle$ and probability density distributions $P(A)$ (histograms).

The signature of a first–order transition in a GCMC simulation is usually a double peak in the probability density distribution $P(N)$, and the peaks should have equal area at coexistence (activity $z_{\rm coex}$). In the vicinity of coexistence between an isotropic and a nematic state with excess of one species, the probability density distribution of one of the nematic order parameters (say $P(Q_{1})$) should furthermore exhibit three peaks: one peak centered at $Q_{1}\approx 0$ and two peaks centered at $Q_{1}=q>0$ and $Q_{1}=-q/2$, where the second peak corresponds to the species 2 or 3 being the majority species. Consider the distribution $P(Q_{\rm max})$ with $Q_{\rm max}=\max\limits_{i}(Q_{i})$: it should exhibit a double–peak structure located at $Q_{\rm max}\approx 0$ and $Q_{\rm max}=q>0$. On the other hand, a transition to a nematic state with excess of two species is best observed via the distribution $P(Q_{\rm min})$ with $Q_{\rm min}=\min\limits_{i}(Q_{i})$. It should likewise exhibit a double–peak structure located at $Q_{\rm min}\approx 0$ and $Q_{\rm min}=q<0$.

A negative nematic order parameter in continuum models for rods would correspond to rods preferentially orienting perpendicular to a director. This has not been reported in the literature in the case of uniaxial hard rods.

III.2 Lattice FMT

A general FMT functional for hard rod mixtures on lattices with arbitrary dimensions has been derived by Lafuente and Cuesta Laf02 ; Laf04 . For the specific case of rods with length $L$ and width 1, the basic definitions and examples for the functionals and their equilibrium properties in two and three dimensions are provided in Ref. Oet16 . We only need the free energy density in the homogeneous case (bulk) for the present work.

[TABLE]

Here, $\beta=1/(k_{\rm B}T)$ is the inverse temperature which is set to 1 and

[TABLE]

is the excess free energy of a zero–dimensional cavity (which can hold no or only one particle) depending on its average occupation $\eta\in[0,1]$. The equilibrium bulk state is found by minimizing $f$ at constant total density $\rho$ with respect to the order parameters $Q\equiv Q_{i}$ and $S\equiv S_{i}$ (for one specific $i$). We have found stable states only for $S=0$ and $Q\geq 0$. Isotropic–nematic coexistence is determined by equating the chemical potential $\mu=\mu_{1}=\mu_{2}=\mu_{3}$ (with $\mu_{i}=\partial f/\partial\rho_{i}$) and pressure $p=\mu\rho-f$ between the isotropic and nematic states.

IV Results

IV.1 GCMC results

For all rod lengths up to $L=25$, we did not see any sign of a double peak in $P(N)$, indicating that a hypothetical density difference between isotropic and nematic states would be well below 1%.

For $L=5$, the distributions $P(Q_{\rm min})$ and $P(Q_{\rm max})$ show a single peak near zero for packing fractions up to 0.84 ($z=25$), indicating a stable isotropic phase. For higher packing fractions, our simple algorithm did not equilibrate the system on a reasonable time scale. Therefore, we cannot exclude the possibility for a nematic phase of some sort at higher packing fraction.

For $L=6$, upon varying $z$ between 2.2 and 2.9, $P(Q_{\rm min})$ transforms from one fairly sharp peak near zero, via a very broad peak (resulting from an overlap of two broad peaks), and finally to a moderately sharp peak at $Q_{\rm min}<0$. This indicates a transition from an isotropic state to a nematic state with two excess species (negative order parameter). Results for $P(Q_{\rm min})$ for several $z$ (illustrating this transition) are shown in Fig. 1(a), and a series of snapshots in Fig. 1(c), which display cuts in the plane defined by the two excess species for $z=3.8$ (in the nematic phase). The snapshots show that the system separates into weakly coupled layers that are essentially only populated by the excess species (i.e. 2D systems). Within these layers, there is no dominance of one of the excess species, in accordance with the absence of demixing for rods with $L=6$ in 2D Ghosh07 . We note that $Q_{\rm min}$ (or $Q_{\rm max}$)—in contrast to $N$—is particulary sensitive to incomplete thermalization of the system (i.e. insufficient exploration of all possible equilibrium configurations), which happens with the standard algorithm used once this layering occurs.

For $L\geq 7$, upon varying $z$, $P(Q_{\rm max})$ transforms from having one peak near zero, via a broad two–peak structure, and finally to a single peak with its maximum at $Q_{\rm max}>0$. This indicates a transition from an isotropic state to a nematic state with one excess species (positive order parameter). $P(Q_{\rm max})$ and $P(\eta)$ are shown in Fig. 2 for rod lengths 8 and 25, among several values of $z$, respectively. The vicinity of coexistence is characterized again by two peaks in $P(Q_{\rm max})$ smeared out broadly, yet the corresponding peaks in $P(\eta)$ show no sign of splitting or of significant broadening. Therefore, the transition must entail a very small density gap, pointing to a weak first order transition, at most.

We additionally define the variances ${\rm Var}(Q_{\rm max}(z))=\langle Q_{\rm max}^{2}\rangle-\langle Q_{\rm max}\rangle^{2}$ (likewise for ${\rm Var}(Q_{\rm min})$). The maximum of ${\rm Var}(Q_{\rm max}(z))$ or ${\rm Var}(Q_{\rm min}(z))$ signals the point of strongest orientational fluctuations, and is hence used to locate the coexistence activity.

For $L=8$ as an example, we have investigated finite–size effects on the behavior of $\langle Q_{\rm max}(z)\rangle$ and ${\rm Var}(Q_{\rm max}(z))$ more closely to examine their influence on the isotropic–nematic transition point. Fig. 3 shows these functions when varying the lattice extension from $M=50$ to 168. For a first order transition, $\langle Q_{\rm max}(z)\rangle$ should exhibit a jump at $z_{\rm coex}$, whereas $z_{\rm coex}$ would be a bifurcation (critical) point for a second order transition. In any case, a finite system size smears the jump or bifurcation. With increasing $M$, $\langle Q_{\rm max}(z)\rangle$ becomes steeper in the coexistence region (see Fig. 3(a)). The curves intersect at a common point for lattice extensions between 50 and 64, indicative of a second order transition. However, the intersection shifts to a lower value of $z$ for $M=100$ and moreso for $M=168$, pointing toward a weak first order transition, consistent with the smeared two–peak structure in $P(Q_{\rm max})$. The corresponding behavior of ${\rm Var}(Q_{\rm max}(z))$ is shown in Fig. 3(b): its peak sharpens and slightly shifts to lower $z$ for increasing $M$. The peak position for the largest lattice extension $M=168$ is what we use to define $z_{\rm coex}$. The peak height in ${\rm Var}(Q_{\rm max}(z))$ shrinks with increasing $M$. This is not quite consistent with a first–order transition for which the peak height should stay constant. One sees that the average $\langle Q_{\rm max}\rangle$ at coexistence (Fig. 3(a)) decreases for increasing lattice size. The order parameter is below 0.1 at coexistence for $M=168$. We additionally performed a finite–size analysis using two–dimensional order parameters known from three–state Potts models and a scaling analysis for the volume susceptibility in Appendix A for the system with $L=8$. It points towards a weak first–order transition, but we cannot rule out a continuous transition with certainty.

Figures 3(c) and (d) show $\langle Q_{\rm max}(z)\rangle$ and ${\rm Var}(Q_{\rm max}(z))$ in the transition region for rod length $L=25$ and large lattice extension ($M=170$). Here, the peak in the variance ${\rm Var}(Q_{\rm max}(z))$ is readily pronounced for this lattice size, and the average value $\langle Q_{\rm max}(z)\rangle$ is around 0.2 at coexistence, i.e. it is larger than for $L=8$, but the increase is very moderate. There is no evidence that a sizeable jump in nematic order takes place at coexistence.

*Other findings: * Our simulations showed a particular finite-size effect that appeared for small systems (e.g. $M=32$ and $L=8$): the system was composed of stacked layers populated with preferably one species, whereby the type of species varied randomly across the layers. We estimated that the entropic gain per particle for this layering configuration decreases inversely proportional to the area $M^{2}$, thus this effect should vanish with larger box sizes (as confirmed by the simulations). We have also implemented biased sampling methods in addition to the standard GCMC algorithm. Specifically, we used successive umbrella sampling Vin05 , which samples equilibrium configurations strictly within a particular interval of a given obervable (e.g. $N$ or $Q$) via biasing—shifting successively over a range of such intervals—to obtain histograms $P(N)$ or $P(Q)$ with higher resolution at the tails of the distribution (the method captures the statistics of rare configurations better than standard GCMC). As $P(N)$ or $P(Q)$ did not show a clear signature for a first–order phase transition, however, these investigations did not lead to better results.

IV.2 Comparison to FMT results

Lattice FMT predicts a strong first order transition for $L\geq 4$ to a nematic state with one excess species (positive order parameter) Oet16 . We did not find stable states with a negative order parameter. The phase diagram resulting from FMT and GCMC is shown in Fig. 4. As described before, we were not able to detect a density gap between isotropic and nematic states in our GCMC simulations. Surprisingly, GCMC simulations and the FMT show similar packing fractions for the coexistence state. We have investigated this more closely by examining $\langle\eta(z)\rangle$ and $\langle Q_{\rm max}\rangle(\langle\eta\rangle)$.

The behavior of $\langle\eta\rangle(z)$ is captured very well by FMT, except for some mild disagreement in the coexistence region, see Fig. 5. A system with very short rods ($L=2$) resembles rather a hard lattice gas (no phase transition), and so the FMT results lie on top of the GCMC data. With longer rods $L=8$ and 25 the high– and low–density limits render good agreement with GCMC, but near coexistence the FMT data show a van der Waals loop characteristic for a first order transition (visible for $L=8$). This is absent in the simulations.

A side–by–side comparison of $\langle Q_{\rm max}\rangle(\langle\eta\rangle)$ unveils stronger disagreement, see Fig. 6. For $L=8$, FMT predicts much more pronounced nematic order for a given $z$. The closeness of the values for $\eta_{\rm coex}$ (GCMC) and $\eta_{\rm coex}$ of the nematic state (FMT) thus might only be a serendipitous accident. For $L=25$, the FMT and GCMC results agree somewhat more. It is actually not clear whether FMT and GCMC would agree in the limit $L\to\infty$ in the vicinity of the nematic transition. FMT renders the correct second and third virial coefficient, and deviates from the fourth on. Although the nematic transition shifts to lower packing fraction with increasing $L$, higher virial coefficients affect the particular location of the transition, as Zwanzig analyzed in the case of a hard rod model with restricted orientational and continuous translational degrees of freedom Zwa63 . This is not the case for hard rod models with continuous orientations, where it can be expected that a second–virial approximation is sufficient for $L\to\infty$ Zwa63 .

V Summary and discussion

Using GCMC simulations, we have found interesting characteristics of the nematic phase transition for hard rods of size $L\times 1\times 1$ on a 3D cubic lattice. We have observed a nematic transition for $L\geq 6$. For exactly $L=6$, the nematic state has a negative order parameter, reflecting the occurrence of two dominating orientations. For $L\geq 7$, the nematic state has a positive order parameter, corresponding to the dominance of one orientations. We have investigated rod lengths up to $L=25$ and have found evidence that the isotropic–nematic transition is of very weak first order for all rod lengths. It was not possible to decide on the type of the transition unambiguously, despite using large systems containing several $10^{5}$ particles. The nematic order parameter $Q$ at the transition is very small, around 0.1, and the probability distribution $P(Q)$ from the GCMC simulations are very broad, pointing to strong fluctuations. Our results complement earlier results on the demixing (pseudonematic) transition for an equivalent system in 2D Ghosh07 , which is presumably of Ising–type and occurs for $L\geq 7$. We have compared our results to lattice FMT (which in the bulk is equivalent to the DiMarzio entropy DiMar61 ) and find that FMT overestimates nematic order strongly. This, however, points to a route for how to improve the FMT. In its current form, the FMT functional is based on the exact description of 0D cavities that can hold up to one particle. It appears sensible to extend this to cavities that can hold two or more particles. Such cavities would capture the correlations between particles lying perpendicular to each other more accurately, and their inclusion would also improve upon the mean–field assumption of DiMarzio.

It is interesting to evaluate the use of lattice models of rods for a possible simplified description of continuum models. We conclude our study by summarizing some findings for continuum models from the literature and relating these to lattice models.

*Continuum hard rods in 2D: * In 2D continuum models with anisotropic particles, simulations have addressed the case of hard needles Fre85 and hard ellipses Fre90 ; Xu13 . Since the particle orientation is indeed a continuous variable—in contrast to lattice models—fluctuations in the average orientation may destroy the long–range order of a nematic state (leading to a quasinematic state); however, anisotropic hard particles (for which position and orientation degrees of freedom are coupled) may still exhibit true, nematic long–range order, see the discussion in Ref. Fre85 . Thus a continuous isotropic–quasinematic transition of Kosterlitz–Thouless type is possible, as well as an isotropic–nematic transition of first order. Interestingly, both types had been found in the hard ellipse system Fre90 (first order for aspect ratio 4 and continuous for aspect ratio 6), though, Ref. Luo14 contradicts this and finds only continuous transitions. This topic may still be an open one. Nevertheless, since the 2D lattice rod model shows an Ising–type, demixing transition, we may conclude that in 2D lattice and continuum models show qualitatively very different behavior with respect to the type of phase transition.

*Restricted-oriented rods and boards in 3D (Zwanzig model): * Zwanzig Zwa63 initiated investigating the packing of hard rods with 3D mutually–orthogonal orientations (a restricted–orientations model) in his seminal work. For long rods, he found that the virial coefficients scaled very differently with aspect ratio compared to those of hard rod models with unrestricted orientations. Therefore, a second–virial approximation is not sufficient to locate the isotropic–nematic transition. This finding should similarly apply to the 3D lattice model investigated by us. Martinez–Raton Mart04 calculated the full phase diagram using an approximate FMT. There, a nematic state with negative order parameter is found in conjunction with layering. Thus the appearance of such nematic states in Mart04 and in the present work seems to be a direct consequence of restricting particle orientations.

*Boards and cuboids in the 3D continuum: * These models release the restriction on orientation, but retain the particle anisometry and shape of the aforementioned Zwanzig models. Refs. Esco05 ; Esco08 simulated bulk properties of hard tetragonal parallelpipeds, havign found a parquet phase for the case of cuboids with two short, symmetrical axes (i.e. rods). This is absent in the Zwanzig models. Ref. Perou13 investigated more deeply the appearance of biaxial phases in these systems.

*Continuum hard rods (spherocylinders) in 3D: * Spherocylinders have proven to be the most sought-after model with regards to hard anisotropic particles in the 3D continuum. On the simulation side, Bolhuis and Frenkel Bol96diss ; Bol97 used MC simulations with $\sim 500$ particles in the canonical ensemble. The nematic phase was found to be stable for aspect ratios $L/D\geq 3.8$. The isotropic–nematic transition has been assumed to be first order but for small aspect ratios $L/D\leq 5.0$ it was not possible to detect a density jump between coexisting states. The order parameter $Q^{\rm nem}(\rho)$ showed weak hysteresis, but a jump in $Q^{\rm nem}$ between coexisting states (a clear signature of a first order transition) was not determined. The transition density is estimated from the kink in the $Q^{\rm nem}$–$\rho$–curve at the (arbitrary) value $Q^{\rm nem}=0.4$. For $L/D\geq 15$ a density jump of $>10\%$ had been found, the absolute value of the density jump was maximal at around $L/D=20$. The nematic order parameter of the nematic phase at coexistence was found to be 0.784 at $L/D=15$ and further increases when the aspect ratio lengthens. So, for $L/D\geq 15$, the transition is very strongly first order. Vink et al. Vin05 employed grand canonical MC techniques. For $L/D=15$ they confirmed the sizable density jump, and at coexistence the probability distribution $P(\rho)$ showed two clearly separated peaks.

On the theory side, the most accurate density functionals for hard anisotropic particles have been derived from Fundamental Measure Theory (FMT) Han09 ; Han10 ; Witt14a ; Witt14b ; Witt15 ; Witt16 . Wittmann et al. used mixed measures to derive a functional that is exact in the low–density limit; for higher densities, typical FMT approximations were employed Witt15 . The corresponding phase diagram for hard spherocylinders showed almost quantitative agreement with Ref. Bol97 . However, the FMT results showed an almost constant density jump between the coexisting isotropic and nematic states for aspect ratios starting from $L/D\gtrsim 3.5$. Thus the transition is strongly first order for all aspect ratios.

It hence appears that for moderate aspect ratios strong orientational fluctuations cause the isotropic–nematic transition to be weakly first order, both on the lattice and in the continuum. The lattice model exhibits in addition the peculiarity of a nematic state with negative order parameter (for $L=6$). The lattice and the continuum stand in contrast to each other for large aspect ratios since we have not found evidence for a strong first order transition. Moreover, it does not seem that a smooth crossover in the topology of the phase diagram is possible from the lattice to the continuum: this was shown by a study using a second–virial density functional with variable orientational degrees of freedom Shu04 .

In the present work, we have not investigated the possibility of high–density phases in this system, owing to limitations of our simple grand–canonical algorithm. These high–density phases could include a completely disordered phase of cubatic–type, which could be similar to the high–density phase in 2D Ghosh07 ; Kun13 . This problem should be treated with optimized algorithms as has been done in 2D Kun13 .

**Note added: After submission of this work, a study on the same system was published on arXiv Raj17 . There, the authors used optimized algorithms for higher densities and could show that the layered nematic phase with $q<0$ not only appears for $L=6$ but also for $L=5$ and $L=7$ at high packing fractions of around 0.9. The corresponding findings for the isotropic–nematic($q>0$) transition are similar to ours. The systems studied were smaller than ours, therefore this transition appeared to be critical in the size regime used. **

**Acknowledgments: ** We thank Richard Vink (Göttingen) for very valuable discussions. This work is supported within the DFG/FNR INTER project “Thin Film Growth” by the Deutsche Forschungsgemeinschaft (DFG), project OE 285/3-1, and by the Landesgraduiertenförderung Baden–Württemberg. Y. Ai thanks the German Academic Exchange Service (DAAD) for a RISE grant in summer 2015.

Appendix A Finite size analysis for $L=8$

Here we apply a general approach for constructing multidimensional order parameters, developed for $q$–state Potts models Bind93 . The discussion in Bind93 , done in the canonical ensemble, can be mapped directly to the purely entropic lattice model treated in the grand–canonical ensemble. This can be done by identifying the temperature with the chemical potential, and the spin states with the rod species. Regarding the symmetries of the order parameter, our lattice model is equivalent to a three–state Potts model where the symmetry–broken phase consists of three equivalent states (each of the Cartesian axes can be a preferred direction). According to Ref. Bind93 , the order parameter dimensionality is two, and orthogonal axes in order parameter space are formed by pairs $(\tilde{Q}_{i},\tilde{S}_{i})$ of unnormalized nematic and biaxial order parameters (Eqs. (1) and (2)):

[TABLE]

where $(ijk)$ is a cyclic permutation of $(123)$. In 2D histograms $P(\tilde{S_{i}},\tilde{Q_{i}})$, a single peak around the origin signals an isotropic state, whereas three peaks with the same distance from the origin and with a threefold symmetry signal a nematic phase. The simultaneous occurrence of all four peaks would correspond to nematic–isotropic coexistence. In Fig. 7, we show $P(\tilde{S_{i}},\tilde{Q_{i}})$ for three values of $z$ in a system with rod length $L=8$ and lattice extension $M=100$. While the isotropic and nematic state can be clearly identified, the histogram near the transition (middle graph in Fig. 7) shows a broad distribution. The four peaks (expected for coexistence) cannot be identified clearly.

In order to facilitate a distinction between first order and continuous transitions, Ref. Bind93 introduces the reduced cumulant

[TABLE]

where $m_{i}$ is the radial distance in the order parameter plane, $m_{i}=\sqrt{\tilde{Q}_{i}^{2}+\tilde{S}_{i}^{2}}$. Cumulants for different lattice extensions $M$ will intersect at $z_{\rm coex}$. In the case of a first order transition, $g_{M}$ will develop a minimum for increasing $M$ below $z_{\rm coex}$, and $z_{\rm coex}-z_{\rm min}$ (where $z_{\rm min}$ is the position of the minimum) scales as $1/M^{3}$. In Fig. 8 we show $g_{M}$ for the system with $L=8$. The common intersection at $z\approx 0.38$ confirms the transition activity identified via the position of the maximum in Var($Q_{\rm max}$) for the largest system ($M=168$). However, our largest system $M=168$ is just large enough to see an indication of a minimum in $g_{M}$; therefore, we cannot do a finite size analysis for the position of the minimum. Nevertheless, the existence of the minimum would point toward a weak first order transition. Investigations much more extensive than the present one would be needed to establish the order of the transition unambiguously.

A scaling analysis corroborates the transition to noncritical behavior for large system sizes. Consider the unnormalized volume susceptibility $\xi_{\rm max}={\rm Var}(\max\limits_{(ijk)}\{N_{i}-(N_{j}+N_{k})/2\})$ which corresponds to the normalized version expressed by ${\rm Var}(Q_{\rm max})$. For a second order transition, one would expect to find a scaling function $\tilde{\xi}(x)$ with the scaling variable $x=(1-z/z_{\rm coex})M^{1/\nu}$ and the relation $\tilde{\xi}(x)=M^{-3-\gamma/\nu}\xi(x)$ Barkema . Here $\gamma$ and $\nu$ are critical exponents in standard notation. The such rescaled susceptibilities for $L=8$ and for different $M$ are shown in Fig. 9 using 3D Ising critical exponents ($\gamma\approx 1.237,\nu\approx 0.630$). The critical exponents of the 3D antiferromagnetic Potts model on a cubic lattice are similar Kol95 . For lattice extensions $M=64...136$, the data collapse on a single scaling curve for $x<0$ (isotropic side). For $x>0$ (nematic side) the collapse is not perfect with deviations coming mainly from small systems. However, the data for $M=168$ clearly deviate from the approximate scaling function, indicating that the transition is not critical anymore. This matches very well with the conclusions drawn from the cumulant analysis.