Phase diagram of a system of hard cubes on the cubic lattice

TL;DR

This study maps the phase diagram of 2x2x2 hard cubes on a cubic lattice, revealing four distinct phases and characterizing the nature of phase transitions through simulations and Landau theory.

Contribution

It provides the first detailed phase diagram of hard cubes on a lattice, identifying four phases and analyzing transition orders with a Landau theoretical framework.

Findings

Four phases identified: disordered, layered, sublattice ordered, columnar.

Disordered-layered transition is continuous; others are discontinuous.

Results near the disordered-layered transition align with the O(3) universality class.

Abstract

We study the phase diagram of a system of $2\times2\times2$ hard cubes on a three dimensional cubic lattice. Using Monte Carlo simulations, we show that the system exhibits four different phases as the density of cubes is increased: disordered, layered, sublattice ordered, and columnar ordered. In the layered phase, the system spontaneously breaks up into parallel slabs of size $2\times L \times L$ where only a very small fraction cubes do not lie wholly within a slab. Within each slab, the cubes are disordered; translation symmetry is thus broken along exactly one principal axis. In the solid-like sublattice ordered phase, the hard cubes preferentially occupy one of eight sublattices of the cubic lattice, breaking translational symmetry along all three principal directions. In the columnar phase, the system spontaneously breaks up into weakly interacting parallel columns of size $2\times 2\times L$ where only a very small fraction cubes do not lie wholly within a column. Within each column, the system is disordered, and thus translational symmetry is broken only along two principal directions. Using finite size scaling, we show that the disordered-layered phase transition is continuous, while the layered-sublattice and sublattice-columnar transitions are discontinuous. We construct a Landau theory written in terms of the layering and columnar order parameters, which is able to describe the different phases that are observed in the simulations and the order of the transitions. Additionally, our results near the disordered-layered transition are consistent with the $O(3)$ universality class perturbed by cubic anisotropy as predicted by the Landau theory.

Figures

Equations

$$\rho =\sum _{i=0}^7{\rho }_i.$$

$$\tilde{\eta }(k_x,k_y,k_z)=\frac{8}{L^3}\sum _{x,y,z}\eta (x,y,z){\mathrm{e}}^{i(k_{x}x+k_{y}y+k_{z}z)}.$$

$$\mathbf{L}=(L_x,L_y,L_z),$$

$$q_1$$

$$q_2$$

$$q_3$$

$$\mathcal{F}=a_L|\mathbf{L}{|}^2+b_L|\mathbf{L}{|}^4+2{\lambda }_L(L_x^2{L}_y^2+L_z^2{L}_y^2+L_x^2{L}_z^2),$$

$$2a_L{L}_x+4b_L{L}_x^3+4L_x(b_L+{\lambda }_L)(L_y^2+L_z^2)$$

$$2a_L{L}_y+4b_L{L}_y^3+4L_y(b_L+{\lambda }_L)(L_x^2+L_z^2)$$

$$2a_L{L}_z+4b_L{L}_z^3+4L_z(b_L+{\lambda }_L)(L_y^2+L_x^2)$$

$$2b_L{l}^2+a_L$$

$$(6b_L+4{\lambda }_L)s^2+a_L$$

$$(4b_L+2{\lambda }_L)c^2+a_L$$

$$\mathcal{H}({\mathbf{L}}_0{)}_{\alpha \beta }={\frac{{\partial }^2\mathcal{F}}{\partial L_{\alpha }\partial L_{\beta }}|}_{\mathbf{L}={\mathbf{L}}_0},$$

$$\mathcal{H}(0,0,0)=\left[\begin{array}{ccc}2a_L& 0& 0\\ 0& 2a_L& 0\\ 0& 0& 2a_L\end{array}\right],$$

$$\mathcal{H}(l,0,0)=\left[\begin{array}{ccc}-4a_L& 0& 0\\ & & \\ 0& -\frac{2a_L{\lambda }_L}{b_L}& 0\\ & & \\ 0& 0& -\frac{2a_L{\lambda }_L}{b_L}\end{array}\right],$$

$$\mathcal{H}(s,s,s)=\left[\begin{array}{ccc}-\frac{4a_L{b}_L}{3b_L+2{\lambda }_L}& -\frac{4a_L(b_L+{\lambda }_L)}{3b_L+2{\lambda }_L}& -\frac{4a_L(b_L+{\lambda }_L)}{3b_L+2{\lambda }_L}\\ & & \\ -\frac{4a_L(b_L+{\lambda }_L)}{3b_L+2{\lambda }_L}& -\frac{4a_L{b}_L}{3b_L+2{\lambda }_L}& -\frac{4a_L(b_L+{\lambda }_L)}{3b_L+2{\lambda }_L}\\ & & \\ -\frac{4a_L(b_L+{\lambda }_L)}{3b_L+2{\lambda }_L}& -\frac{4a_L(b_L+{\lambda }_L)}{3b_L+2{\lambda }_L}& -\frac{4a_L{b}_L}{3b_L+2{\lambda }_L}\end{array}\right],$$

$$\mathcal{H}(c,c,0)=\left[\begin{array}{ccc}-\frac{4a_L{b}_L}{2b_L+{\lambda }_L}& -\frac{4a_L(b_L+{\lambda }_L)}{2b_L+{\lambda }_L}& 0\\ & & \\ -\frac{4a_L(b_L+{\lambda }_L)}{2b_L+{\lambda }_L}& -\frac{4a_L{b}_L}{2b_L+{\lambda }_L}& 0\\ & & \\ 0& 0& -\frac{2a_L{\lambda }_L}{2b_L+{\lambda }_L}\end{array}\right],$$

$$\mathcal{F}$$

$$\mathcal{H}({\mathbf{L}}_0,{\mathbf{C}}_0{)}_{\alpha \beta }={\frac{{\partial }^2\mathcal{F}}{\partial {\varphi }_{\alpha }\partial {\varphi }_{\beta }}|}_{\mathbf{L}={\mathbf{L}}_0,\mathbf{C}={\mathbf{C}}_0},$$

$${\chi }_i$$

$$U_i$$

$$ϵ=\mu -{\mu }_c,$$

$$q_1(ϵ,L)$$

$${\chi }_1(ϵ,L)$$

$$U_1(ϵ,L)$$

$$Q_z(n){\mathrm{e}}^{i{\theta }_z(n)}={\varphi }_{er}^z(n)-{\varphi }_{or}^z(n)+i\left[{\varphi }_{ec}^z(n)-{\varphi }_{oc}^z(n)\right],$$

$${\mathcal{L}}^{\alpha }(z,z^{\prime })$$

$$F^{\alpha }(z,z^{\prime })$$

$${\mathcal{L}}_0^{||}(z)=[{\Omega }_p(L{)}^{L/2}{]}^{L/2},$$

$${\mathcal{L}}_1^{||}(z)=\frac{L^3}{4}[{\Omega }_p(L){]}^{\frac{L^2}{4}}\left({\left[\frac{{\Omega }_o(L-2)}{{\Omega }_p(L)}\right]}^2+{\left[\frac{{\Omega }_o(L-2)}{{\Omega }_p(L)}\right]}^4\right),$$

$$\displaystyle\mathcal{L}^{ro}_{1}(z)=[\Omega_{p}(L)^{L/2}]^{\frac{L}{2}-2}\bigg{[}\left(\frac{L}{2}-2\right)\frac{L^{2}}{2}\times$$

$$\left\{[{\Omega }_o(L-2){\Omega }_p(L{)}^{\frac{L}{2}-1}{]}^2+[{\Omega }_o(L-2{)}^2{\Omega }_p(L{)}^{\frac{L}{2}-2}{]}^2\right\}$$

$$\displaystyle+2L^{2}\Omega_{o}(L-2)^{3}\Omega_{p}(L)^{L-3}\bigg{]}.$$

$$\Delta F^{||,ro}(z,z^{\prime })$$

$$\Delta F^{||,ro}(z,z^{\prime })=L^2{z}^{\prime }{\left({\left[\frac{{\Omega }_o(L-2)}{{\Omega }_p(L)}\right]}^2-\left[\frac{{\Omega }_o(L-2)}{{\Omega }_p(L)}\right]\right)}^2.$$

Full text

I Introduction

Models with only excluded volume interactions have been studied for a long time as the simplest statistical models of thermodynamic phase transitions. In these models, the phases and phase transitions are completely determined by the shape and density of the particles, and temperature plays no role. Well-known examples include the isotropic-nematic transition in long needles Onsager (1949); de Gennes and Prost (1995), the freezing transition in hard spheres Alder and Wainwright (1957); Wood and Jacobson (1957); Isobe and Krauth (2015), and phase transitions in lattice models with nearest neighbour exclusion Domb (1958); Burley (1960, 1961). Experimental systems exhibiting such entropy-driven phase transitions include those between nematic, smectic and cholesteric phases in liquid crystals de Gennes and Prost (1995), nanotube gels Islam et al. (2004), and suspensions of tobacco mosaic virus Fraden et al. (1989). Other examples exhibiting such transitions include adsorbed gas molecules on metallic surfaces Taylor et al. (1985); Patrykiejew et al. (2000); Dünweg et al. (1991), as well as colloidal suspensions such as polymethyl methacrylate (PMMA) suspended in poly-12-hydroxystearic acid Pusey and van Megen (1986). Despite a long history of study, a general understanding of the dependence of the nature of the emergent phases on the shapes of the particles, as well as the order of appearance of the phases with increasing density, is lacking.

The system of hard spheres in three dimensions was one of the first numerically studied systems Alder and Wainwright (1957); Wood and Jacobson (1957) to show such an entropy-driven phase transition. It undergoes a first-order transition from a fluid phase to a solid phase with face centred cubic packing Gast and Russel (1998); Zhu et al. (1997). More complicated shapes such as cubes, rhombohedra Pu et al. (2006) or in general three dimensional regular polyhedra or corner-rounded polyhedra Gantapara et al. (2013); Marechal et al. (2012); Batten et al. (2010) have been studied as more realistic models for experimental self-assembling systems Sun et al. (2000); Torquato and Jiao (2009); Hanrath et al. (2009), applications to drug delivery where shape of the carrier may decide its effectiveness Champion et al. (2007), biological material like immunoglobin Li et al. (2008), molecular logic gates Soe et al. (2011a, b); Godlewski et al. (2013), etc. Being able to predict the macroscopic material behaviour from knowing its constituent building blocks would help to engineer the synthesis of materials with prescribed properties Agarwal and Escobedo (2011); Damasceno et al. (2012).

Of the non-spherical shapes, the simplest is a cube, which has the additional feature that cubes can be packed to fill all space. Theoretical studies in the continuum have focused on two cases: unoriented cubes whose faces are free to orient in any direction, and parallel hard cubes whose axes are parallel to the coordinate axes. The system of unoriented cubes was shown, using Monte Carlo and event driven molecular dynamics simulations, to undergo a first order freezing transition from a fluid to a solid phase at a critical packing fraction $\eta\approx 0.51$ Smallenburg et al. (2012). Other simulations, however, found a cubatic phase that is sandwiched between the fluid and solid phases for packing fractions in the range $0.52<\eta<0.57$ Agarwal and Escobedo (2011). It has been claimed in Ref. Smallenburg et al. (2012) that the cubatic phase is a finite-size artifact. In the case of parallel hard cubes, early work focused on finding the equation of state using high-density expansions Hoover (1964), and low-density virial expansion up to the seventh virial coefficient Zwanzig (1956); Hoover and Rocco (1962). Monte Carlo simulations show that the system of parallel hard cubes undergoes a continuous freezing transition from a disordered fluid phase to a solid phase at density $\rho\approx 0.48$ Jagla (1998); Groh and Mulder (2001). The data near the critical point are consistent with the three-dimensional Heisenberg universality class Groh and Mulder (2001). These results are consistent with theoretical predictions using density functional theory Belli et al. (2012). Within this theory, the columnar phase is found to be not a stable phase at high densities Groh and Mulder (2001); Belli et al. (2012). Thus, it would appear that parallel hard cubes in the continuum show only one phase transition and the high density phase is crystalline.

Hard-core lattice gas models also provide interesting examples of entropy-driven phase transitions, and like in the continuum, have rich phase diagrams. Rigorous results are known for dimer gas Heilmann and Lieb (1970), hard triangles at full packing Verberkmoes and Nienhuis (1999), hard hexagons Baxter (1980), long rods Disertori and Giuliani (2013) and hard plates in three dimensions Disertori et al. (2018). For other shapes, Monte Carlo are more reliable than predictions based on approximate theories. Examples include rods Ghosh and Dhar (2007); Matoz-Fernandez et al. (2008); Kundu et al. (2013), pentamers Eisenberg and Baram (2000), tetronimos Barnes et al. (2009), squares Bellemans and Nigam (1966, 1967); Ramola and Dhar (2012); Ree and Chesnut (1966); Nath et al. (2016); Mandal et al. (2017), etc. In three dimensions, the results are much fewer and the detailed phase diagram is known only for long rods Gschwind et al. (2017); Vigneshwar et al. (2017). Monte Carlo simulations with local moves are often inefficient in equilibriating the system when the excluded volume is large or packing fraction is high. This puts a restriction on the kind of systems one can study. Recently, we have introduced an efficient Monte Carlo algorithm with cluster moves that has helped in overcoming these difficulties Kundu et al. (2012, 2013); Ramola et al. (2015). We have used this to determine the unexpectedly complex phase structure, and nature of the phase transitions in systems like rods in two Kundu et al. (2013) and three dimensions Vigneshwar et al. (2017), hard rectangles Kundu and Rajesh (2014, 2015a, 2015b); Nath et al. (2016, 2015), discretized discs Nath and Rajesh (2014, 2016), Y-shaped molecules Mandal et al. (2018), as well as mixtures Ramola et al. (2015); Kundu et al. (2016) of hard objects.

In this paper, we study a system of $2\times 2\times 2$ hard cubes on the cubic lattice using this cluster algorithm Kundu et al. (2012, 2013); Ramola et al. (2015). The positions of cubes are now discrete. The differences between the phase structure of this discrete problem and the corresponding continuum one are not well understood. Earlier Monte Carlo studies of the discrete problem Panagiotopoulos (2005) found no phase transitions for densities upto full packing (in Ref. Panagiotopoulos (2005), the problem of cubes correspond to $\sigma=2$). On the other hand, for cubes with sides of length two, the approximate density functional theory predicts that there should be a transition from a disordered phase to a layered phase at low densities and from a layered phase to a columnar phase at higher densities. When the length of a side is six, the theory predicts a transition from a disordered phase to a solid, and then to two types of columnar phases Lafuente and Cuesta (2003). Simulations of a mixture of cubes of sizes two, and four or six, show a demixing transition Dijkstra and Frenkel (1994) in contradiction to predictions from density functional theory Lafuente and Cuesta (2002). However, the prediction for a pure system of cubes have not, to our knowledge, been tested in large scale simulations.

Here, we find that this system of cubes goes through four distinct phases as the density of cubes is increased: disordered, layered, sublattice ordered, and columnar ordered. In the layered phase, the system spontaneously breaks up into parallel slabs of size $2\times L\times L$ which are preferentially occupied by cubes. Within each slab, the cubes are disordered; translation symmetry is thus broken along exactly one principal axis. In the solid-like sublattice ordered phase, the hard cubes preferentially occupy one of eight sublattices of the cubic lattice, breaking translational symmetry along all three principal directions. In the columnar phase, the system spontaneously breaks up into weakly interacting parallel columns of size $2\times 2\times L$ which are preferentially occupied by cubes. Within each column, the system is disordered, and the columns break translational symmetry along along two principal directions. By studying systems of different sizes, we argue that the disordered-layered phase transition is continuous, while the layered-sublattice and sublattice-columnar transitions are discontinuous. We construct a Landau theory written in terms of the layering order parameter $\mathbf{L}$ and columnar order parameter $\mathbf{C}$ which is able to describe the different phases that are observed in the simulations and the order of the transitions. Additionally, our results near the disordered-layered transition are consistent with the Landau theory prediction of scaling behaviour in the $O(3)$ universality class perturbed by cubic anisotropy.

The remainder of the paper is organized as follows. Section II defines the model precisely and describes the grand canonical Monte Carlo scheme that is used to simulate the system. Section III describes the different phases – disordered, layered, sublattice ordered and columnar ordered– that we observe in our simulations. In Sec. IV, we understand our simulation results in terms of a Landau theory approach. Sections V - VII characterize the different phase transitions that occur in this system. Section VIII discusses the long lived metastable states that we observe at densities close to full packing. This section also presents a perturbation expansion that allows us to argue that the high density phase will be columnar. Finally, Sec. IX contains a discussion of our results.

II Model & Algorithm

Consider a $L\times L\times L$ cubic lattice with periodic boundary conditions and even $L$. The lattice may be occupied by cubes of size $2\times 2\times 2$ (i.e having side-length of $2$ lattice spacings) whose positions are in registry with the lattice sites. We associate a weight $z=e^{\mu}$ with each cube, where $z$ is the activity and $\mu$ is the chemical potential. The cubes interact through only excluded volume interaction, i.e. no two cubes can overlap in volume. For a cube, we identify the vertex with minimum $x$-, $y$-, and $z$-coordinates as its head. The configuration of the system can thus be fully specified by the spatial coordinates of the heads of all the cubes in the system.

We study the model using grand canonical Monte Carlo simulations implementing an algorithm which includes cluster moves Kundu et al. (2012, 2013); Ramola et al. (2015); Nath and Rajesh (2014) that help in equilibrating systems of hard particles with large excluded volume at densities close to full packing Kundu et al. (2012, 2013) or at full packing Ramola et al. (2015). Below, we briefly summarise the algorithm.

Choose at random one of the $3L^{2}$ rows, where each row consists of $L$ consecutive sites in any one direction. Evaporate all the cubes whose “heads” lie on this row. The row now consists of empty intervals separated from each other by sites that cannot be occupied by the head of a cube due to the hard constraints arising from cubes in neighbouring rows. The empty intervals are reoccupied by new configurations of cubes with the correct equilibrium probabilities. The calculation of these probabilities reduces to a one dimensional problem which may be solved exactly by standard transfer matrix methods (see Refs. Kundu et al. (2013); Ramola et al. (2015); Nath and Rajesh (2014); Vigneshwar et al. (2017) for details). This evaporation and deposition move satisfies detailed balance as the transition rates depend only on the equilibrium probabilities of the new configuration. We use a parallelised version of the algorithm described above, which exploits the fact that the rows separated from each other by a distance two can be updated independently and concurrently. We check for equilibration by taking different initial configurations of the system that correspond to different phases and confirming that the results are independent of the initial configuration. We find that the algorithm is able to equilibrate systems with density upto $\approx 0.95$ for $L\gtrsim 100$, though slightly larger densities may be attained for smaller systems.

III Different Phases

We first define and describe the phases that we observe in our simulations. To do so, it is convenient to divide the lattice into 8 sublattices, depending on whether each $x$-, $y$-, and $z$- coordinates are odd or even. A site with coordinates $(x,y,z)$ belongs to sublattice whose binary representation is ($x$ mod $2$) ($y$ mod $2$) ($z$ mod $2$), as shown in Fig. 1. We define $\rho_{i}$ to be the fraction of lattice sites occupied by the cubes whose heads lie on the sublattice $i$. The total density of the system $\rho$ is then

[TABLE]

Further, let $\eta(x,y,z)$ be equal to $1$ if $(x,y,z)$ is occupied by a head of the cube, and be equal to [math] otherwise. Consider the Fourier transform

[TABLE]

We define the order parameter $\mathbf{L}$ as

[TABLE]

where $L_{x}=\tilde{\eta}(\pi,0,0)$, $L_{y}=\tilde{\eta}(0,\pi,0)$ and $L_{z}=\tilde{\eta}(0,0,\pi)$. A non-zero value in $L_{x}$ will imply that there is translational order of period two in $x$ direction. Similar interpretations hold for $L_{y}$ and $L_{z}$. The $\mathbf{L}$ vector is thus a measure of the layering tendency of the system in each cartesian direction, and we shall refer to it as the layering vector. In a layered phase, only one cartesian component of $\mathbf{L}$ is expected to be nonzero in the thermodynamic limit. In contrast, a columnar-ordered phase is characterized by a layering vector with two nonzero cartesian components. Finally, a solid-like sublattice-ordered phase is characterized by a layering vector with all three components nonzero.

Note that $\mathbf{L}$ serves as a “faithful” order parameter for each of these three phases: In the layered case, it correctly distinguishes between the six symmetry-related states of the system (corresponding to the two possible layered states for layering along each of the three cartesian directions) by taking on the six symmetry related values $(\pm|L|,0,0)$, $(0,\pm|L|,0)$, and $(0,0,\pm|L|)$. In the columnar-ordered case, it correctly distinguishes between the twelve symmetry-related columnar states by taking on the twelve symmetry-related values $(\pm|L|,\pm|L|,0)$, $(0,\pm|L|,\pm|L|)$, and $(\pm|L|,0,\pm|L|)$. Finally, the eight symmetry-related sublattice ordered states are correctly described by the eight symmetry-related values $(\pm|L|,\pm|L|,\pm|L|)$.

To characterize these phases, it is also useful to define two other measures of spontanously broken symmetry: the columnar vector $\mathbf{C}$ whose components are given by $C_{x}=\tilde{\eta}(0,\pi,\pi)$, $C_{y}=\tilde{\eta}(\pi,0,\pi)$ and $C_{z}=\tilde{\eta}(\pi,\pi,0)$, and the sublattice scalar $\phi=\tilde{\eta}(\pi,\pi,\pi)$. In contrast to $\mathbf{L}$, neither $\mathbf{C}$ nor $\phi$ fully distinguish between the broken symmetry states in which they are nonzero in the thermodynamic limit. This is clear since $\mathbf{C}$ is expected to be nonzero in the columnar ordered phase and the sublattice ordered phase, but does not correctly distinguish between the twelve symmetry-related columnar states or the eight symmetry-related sublattice ordered states of the system. Similarly, $\phi$ is expected to be nonzero in the sublattice ordered phase, but does not correctly distinguish between the eight symmetry-related sublattice ordered states.

The underlying reason for this distinction between $\mathbf{L}$ on the one hand, and $\mathbf{C}$ and $\phi$ on the other, is clarified considerably if we pass from the globally defined quantities $\mathbf{L}$, $\mathbf{C}$ and $\phi$, to the corresponding local fields $\mathbf{L}(\vec{r})$, $\mathbf{C}(\vec{r})$ and $\phi(\vec{r})$. These local fields should be thought of as being the coarse-grained variables (coarse-grained over a linear scale of a few lattice spacings) whose sum over the entire volume gives the corresponding global variables. Thinking in terms of these local fields, we see that $C_{x}(\vec{r})\sim L_{y}(\vec{r})L_{z}(\vec{r})$ (and similarly for the other components). This is related to the fact that the composite variable $L_{y}L_{z}$ acts as a field that couples linearly to $C_{x}$ in a Landau-type description of spontaneous symmetry breaking. Likewise, $\phi(\vec{r})\sim L_{x}(\vec{r})L_{y}(\vec{r})L_{z}(\vec{r})$. Thus, $\mathbf{C}(\vec{r})$ and $\phi(\vec{r})$ take on mean values set by composite variables formed from the components of the local layering vector, which emerges as the fundamental quantity for describing the broken symmetries of the system. It is therefore not surprising that the corresponding global variables ${\mathbf{C}}$ and $\phi$ do not fully distinguish between different symmetry-related states with spontaneous columnar or sublattice order. This also suggests that a Landau theory for all three broken symmetry phases should involve $\mathbf{L}$ as the key variable, although we shall see below that the symmetry-allowed couplings between $\mathbf{L}$ and $\mathbf{C}$ and $\phi$ can also play a crucial role in determining the structure of the phase diagram.

In our simulations, we measure the magnitudes of the global variables $\mathbf{L}$, $\mathbf{C}$ and $\phi$:

[TABLE]

In addition, we monitor the joint probability distribution (histogram) of $L_{x}$, $L_{y}$ and $L_{z}$ in order to visualize the nature of the symmetry breaking present in various ordered states.

The variation of $q_{i}$ with density $\rho$ is shown in Fig. 2. For low densities $q_{i}\rightarrow 0$ in the thermodynamic limit for $i=1,2,3$ and the system is in a disordered phase. As the density is increased, $q_{1}$ becomes non-zero in the thermodynamic limit when the density crosses $\rho\approx 0.718$, signalling the onset of spontaneous layering, while $q_{2}$ and $q_{3}$ continue to be zero in the thermodynamic limit. Upon further increasing the density, both $q_{2}$ and $q_{3}$ become nonzero in the thermodyamic limit when the density increases beyond $\rho\approx 0.79$. This corresponds to the onset of a cystalline phase with spontaneous sublattice ordering. Finally, when the density goes beyond $\rho\approx 0.957$, $q_{3}$ becomes zero, while $q_{2}$ and $q_{1}$ remain nonzero, corresponding to columnar order.

Below, we describe the behaviour of the system in each of these phases in some more detail.

Disordered Phase: At low densities, the cubes are in a disordered phase in which the cubes are far apart and there is no ordering. All the mean sublattice densities are equal, i.e., $\rho_{i}\approx\rho/8$, for $i=0,\dots,7$. In the disordered phase, all components of $\mathbf{L}$ tend to zero in the limit of large system sizes [see Eq. (3) for definition] and the system retains all the symmetries of the underlying cubic lattice.

Layered Phase: In the layered phase, translational symmetry is broken in only one direction. The cubes preferentially occupy either odd or even planes normal to this direction. This may be seen by examining snapshots of randomly chosen pairs of even and odd planes in the three directions as shown in Fig. 3, where the eight different colours represent cubes whose heads on a particular sublattice. Grey colour represents sites that are occupied by cubes whose heads are on neighbouring planes. In Fig. 3(a)-(d), showing the snapshots of randomly chosen even and odd $yz$ and $xz$ planes, there are approximately equal number of coloured cubes and grey cubes, showing both odd and even $yz$- and $xz$-planes are equally occupied. On the other hand, it can be seen that Fig. 3(e), showing snapshot of a randomly chosen even $xy$ plane, has much larger number of coloured squares than grey squares, while Fig. 3(f), showing snapshot of a randomly chosen odd $xy$ plane, is mostly grey, showing that in this configuration, the heads of cubes preferentially occupy even $xy$-planes.

The breaking of translational invariance is also quantitatively reflected from the time evolution of the eight sublattice densities and $L_{x},L_{y},L_{z}$, as shown in Fig. 4(a) and (b) respectively. From Fig. 4(a), we see that four sublattices are preferentially occupied. From Fig. 4(b), we also see that one of the components of $\mathbf{L}$ is larger than the other two, i.e., $|L_{z}|\gg|L_{x}|\approx|L_{y}|$, confirming that the system is layered in the $z$-direction. More precisely, one component of $\mathbf{L}$ remains nonzero in the thermodynamic limit, and the other two are zero.

Sublattice Phase: In the sublattice phase, translational symmetry is broken in all three principal directions of the cubic lattice. In this phase, the cubes preferentially occupy one of the eight sublattices. This may be seen by examining the snapshots of randomly chosen pairs of even and odd planes in the three directions as shown in Fig. 5. It may be seen that in each of the directions, one of the planes has a larger number of cubes, compared to the grey squares. We see that in this case the cubes preferentially occupies simultaneously even $yz$, odd $xz$ and even $xy$-planes, which implies most of the cubes occupy sublattice 2.

The breaking of translational symmetry is reflected in the time evolution of the eight sublattice densities and $|L_{x}|,|L_{y}|,|L_{z}|$ as shown in Fig. 6(a) and (b) respectively. In Fig. 6(a), sublattice $2$ is preferentially occupied over the seven. From Fig. 6(b), we also see that $|L_{x}|,|L_{y}|,|L_{z}|$ are non-zero and equal, i.e., $|L_{x}|\approx|L_{y}|\approx|L_{z}|\gg 0$, confirming that the system has sublattice order.

Columnar Phase: The system is in a columnar phase at large densities. In the columnar phase, the system breaks translational symmetry along two directions and the heads of the cubes preferentially occupy two sublattices. This may be seen by examining the snapshots of the planes in the three directions, as shown in Fig. 7. From Figs. 7(d) and (f), corresponding to snapshots of odd $xz$-plane and odd $xy$-plane respectively, it may be seen that these planes contain very few heads of cubes. Thus, most cubes have heads with even $y$-coordinate and even $z$-coordinate. If now the $x$-coordinate has no definite parity, then the phase will be columnar, else it will be a sublattice phase. From the snapshots of even and odd $yz$-planes, shown in Figs. 7(a) and (b), it can be seen that both planes have roughly equal number of heads of cubes, showing that the $x$-coordinate has no definite parity. This feature may also be observed from the snapshots shown in Figs. 7(c) and (e) of even $xz$- and even $xy$- planes, where two colours are seen in each snapshot corresponding to even $x$ and odd $x$.

From Fig. 8(a), we see that the sublattice [math] and $4$ are preferentially occupied over the six sublattices corresponding to the heads of most of the cubes having odd $y$ and $z$-coordinates. From Fig. 8(b), we see that two order parameters are large compared to one, i.e., $|L_{y}|\approx|L_{z}|\gg|L_{x}|$, this implies that translation symmetry is broken in both the $y$- and $z$-directions. The columnar phase can be visualised as a set of tubes extending along the $x$-direction, in which the cubes can slide along.

IV Landau theory for $2\times 2\times 2$ cubes

In this section, we formulate a Landau-type theoretical description of the phases found in Sec. III. As noted earlier, spontaneous layering, sublattice ordering and columnar ordering are faithfully described by the layering vector $\mathbf{L}$, defined in Eq. (3), with the columnar vector $\mathbf{C}$ and the sublattice scalar $\phi$ more naturally thought of as composite objects constructed from the local layering order parameter field. It is therefore natural to try to construct the Landau theory in terms of the order parameter vector $\mathbf{L}$. Here, we demonstrate that while such a Landau theory correctly captures the low density disordered phase, the layered phase, and the sublattice phase, as well as phase transitions between them, it does not allow for the possibility of a columnar phase. To account for the columnar phase, we include the symmetry-allowed couplings to the columnar vector $\mathbf{C}$ and write down a coupled theory for $\mathbf{L}$ and $\mathbf{C}$. This augmented Landau theory correctly predicts the existence of a columnar phase, as well as the nature of the transition to the columnar phase.

We start by constructing the functional only in terms of $\mathbf{L}$. The symmetries of the Landau functional $\mathcal{F}(\{L_{\alpha}\})$ are that it is invariant under $\{L_{\alpha}\leftrightarrow-L_{\alpha}\}$ for $\alpha=(x,y,z)$ and cyclical permutations of the indices $(x,y,z)$. With these constraints, the most general functional is

[TABLE]

where we have truncated the expansion upto fourth order. To make sure that $\mathcal{F}$ goes to $+\infty$ when $|\mathbf{L}|\rightarrow\infty$, we require that $b_{L}>0$ and $\lambda_{L}>-3b_{L}/2$. The Landau theory in Eq. (7) is that of $O(3)$ model with a cubic anisotropy.

The equilibrium phase is obtained from the global minimum of $\mathcal{F}$, and is obtained from the solutions of $\nabla_{L}\mathcal{F}=0$. In component form, these equations are

[TABLE]

The solutions to Eqs. (8)–(10) may be found in closed form. We find that the solutions are of the form $(0,0,0)$, $(l,0,0)$, $(s,s,s)$ and $(c,c,0)$ or its cyclic permutations. Substituting into Eqs. (8)–(10), the equations satisfied by $l,s,c$ are

[TABLE]

The stability of the phases is determined by examining the Hessian $\mathcal{H}(\mathbf{L}_{0})$ defined as

[TABLE]

where $\alpha$ and $\beta$ run over the indices $(x,y,z)$. For $\mathbf{L}_{0}$ to be local minimum or locally stable, we require that the three eigenvalues of the Hessian, calculated at $\mathbf{L}_{0}$, are all positive.

For the disordered phase $(0,0,0)$, the Hessian is

[TABLE]

whose three eigenvalues are all equal to $2a_{L}$. For the eigenvalues to be positive, we require that $a_{L}>0$.

For the layered solution $(l,0,0)$, the Hessian is

[TABLE]

whose eigenvalues are the diagonal entries in Eq. (16). For the eigenvalues to be positive, we require that $a_{L}<0$ and $\lambda_{L}>0$.

For the sublattice solution $(s,s,s)$, the Hessian is

[TABLE]

whose eigenvalues are $-4a_{L}$, $4a_{L}\lambda_{L}/(3b_{L}+2\lambda_{L})$, and $4a_{L}\lambda_{L}/(3b_{L}+2\lambda_{L})$. For the eigenvalues to be positive, we require that $a_{L}<0$ and $\lambda_{L}<0$.

For columnar solution $(c,c,0)$, the Hessian is

[TABLE]

whose eigenvalues are $-4a_{L}$, $4a_{L}\lambda_{L}/(2b_{L}+\lambda_{L})$, and $-2a_{L}\lambda_{L}/(2b_{L}+\lambda_{L})$. The ratio of the second and third eigenvalues is $-2$. This implies that the three eigenvalues cannot be made simultaneously positive. Thus, the columnar solution is not a stable solution.

From the above analysis, we find that there exists a unique stable solution for each choice of $a_{L}$ and $\lambda_{L}$. For $a_{L}>0$ and any $\lambda_{L}$, the only stable phase is the disordered phase where $\mathbf{L}=0$. For $a_{L}<0$ and $\lambda_{L}>0$, we find that stable solution is a layered phase, where $\mathbf{L}$ is a one-component vector of the form $(l,0,0)$. For the case where $a_{L}<0$ and $\lambda_{L}<0$, the stable solution is a sublattice phase, where $\mathbf{L}$ is vector of the form $(s,s,s)$.

These observations are summarized in the phase diagram shown in Fig. 9. The disordered-layered transition and disordered-sublattice transitions are both continuous and, within Landau theory, belong to the universality class of the $O(3)$ model with cubic anisotropy. On the other hand, the sublattice-layered transition is discontinuous, as the orientation of the $\mathbf{L}$ vector changes abruptly from along one of the axes to $(1,1,1)$ or an equivalent direction.

The simplest Landau-type theory described in Eq. (7) predicts disordered, layered and sublattice phases, but disallows a columnar phase. To construct a minimal theory that predicts all the phases that are seen in the simulations, we extend the functional in Eq. (7) to explicitly depend on the columnar vector $\mathbf{C}(\vec{r})$ defined earlier. Here, $\mathbf{C}(\vec{r})$ should be thought of as an independent coarse-grained vector field since, $\langle C_{x}(\vec{r})\rangle=\left\langle L_{y}(\vec{r})L_{z}(\vec{r})\right\rangle$, cannot be fully determined in terms of $\langle L_{y}(\vec{r})\rangle$ and $\langle L_{z}(\vec{r})\rangle$. The Landau functional $\mathcal{F}(\{L_{\alpha},C_{\alpha}\})$ should be invariant under $L_{\alpha}\leftrightarrow-L_{\alpha}$, and $C_{\alpha}\leftrightarrow-C_{\alpha}$ for $\alpha=(x,y,z)$ and under cyclical permutations of the indices $(x,y,z)$. The augmented functional, truncated upto fourth order is:

[TABLE]

where $\mu$ couples the $\mathbf{L}$ and $\mathbf{C}$ vectors. We restrict ourselves to $\mu>0$ since this correctly describes the situation in a columnar ordered configuration of our system of cubes (with our definitions of these vectors, it is easy to see that $C_{x}$ has the same sign as the product $L_{y}L_{z}$ in a columnar ordered state with columnar vector pointing in the $x$ direction, and similarly for columnar vectors pointing in the other cartesian directions). As we demonstrate below, this augmented Landau theory now accounts for the presence of a stable columnar phase in addition to the other stable phases already obtained by thinking entirely in terms of $\mathbf{L}$

The extended functional in Eq. (19) has seven independent parameters and six variables. Deriving analytic equations of the phase boundary is not possible as that would require us to simultaneously solve six coupled equations. Rather, we focus on showing that there are parameter regimes for which the columnar phase, as well as the other phases exist and are stable. This is achieved by assigning numerical values to the parameters values and solving the coupled equations for equilibrium numerically. The stability is checked using the Hessian $\mathcal{H}(\mathbf{L},\mathbf{C})$:

[TABLE]

where $\phi$ runs over the components of the vectors $\mathbf{L}$ and $\mathbf{C}$, is now a $6\times 6$ matrix.

For analysing the functional in Eq. (19) to determine its global minima, we consider the four different cases discussed below. For each of these cases, we set $\mu=2$ and fix the parameters $b_{L}$ and $b_{c}$ to be large and positive ($b_{L}=b_{c}=8$).

Case 1. $a_{L}>0$, $a_{c}>0$: In this case, in the absence of the coupling ($\mu=0$), both $\mathbf{L}=0$, and $\mathbf{C}=0$. For small positive $\mu$, we expect the system to be still in the disordered phase. We confirm this by setting $a_{L}=a_{c}=1.2$, and treating $\lambda_{L}$ and $\lambda_{c}$ as free parameters. For this case, whatever be the values and sign of $\lambda_{L}$ and $\lambda_{c}$, we find that the disordered phase is the only stable phase.

Case 2. $a_{L}>0$, $a_{c}<0$: In this case, we expect that $\mathbf{L}=0$, and $\mathbf{C}\neq 0$. Such solutions are unphysical, and we expect that the mapping from the microscopic variables of the model to the parameters of the Landau theory is such that, this regime is never reached.

Case 3. $a_{L}<0$, $a_{c}>0$: In this case, in the absence of coupling ($\mu=0$), the $\mathbf{L}$ shows both layered and sublattice phases depending on the sign of $\lambda_{L}$. When $\mu\neq 0$, for these phases to be valid, $\mathbf{C}$ should be zero in the layered phase and have three non-zero components in the sublattice phase. We confirm that this is indeed the case by determining numerically the phase diagram for $a_{L}=-1.2$, $a_{c}=1.2$. We find that the system is layered when $\lambda_{L}>0$, and has sublattice order when $\lambda_{L}<0$, irrespective of the sign of $\lambda_{c}$. The schematic $\lambda_{L}$-$\lambda_{c}$ phase diagram for this case, obtained by minimising the free energy at different sample phase points, is summarised in Fig. 10.

Case 4. $a_{L}<0$, $a_{c}<0$: In this case, we expect that when $\lambda_{c}>0$, then $\mathbf{C}$ could exist in a layered phase, making it possible for the system to be in a columnar phase. We determine the phase diagram for $a_{L}=a_{c}=-1.2$. We find that for $\lambda_{c}>0$ and $\lambda_{L}>0$, there is a regime where the system is in a columnar phase. For large $\lambda_{c}$ and $\lambda_{L}$, there are spurious unphysical solutions. For the other cases, the system is in a sublattice phase. The schematic $\lambda_{L}$-$\lambda_{c}$ phase diagram for this case, obtained by minimising the free energy at different sample phase points, is summarised in Fig. 11.

V Disordered-Layered Transition

In this section, we describe the phase transition from disordered phase to layered phase. To do so, we first define susceptibility $\chi_{i}$, and Binder cumulant $U_{i}$ associated with the order parameter $q_{i}$ [see Eq. (4)–(6) for definition] as

[TABLE]

where $c_{1}=c_{2}=9/15$ and $c_{3}=1/3$. The values of $c_{i}$ are chosen so that the Binder cumulant is zero in the disordered phase. We also define the deviation from the critical point as

[TABLE]

where $\mu_{c}$ is the critical value of the chemical potential.

We find that the disordered-layered transition is continuous. A suitable order parameter to study this transition is $q_{1}$ which is zero in the disordered phase and non-zero in layered phase. The critical behaviour may be obtained by studying the non-analytic behaviour of the different physical quantities, which, near the transition, is captured by the finite size scaling behaviour:

[TABLE]

where $f_{q}$, $f_{\chi}$ and $f_{U}$ are scaling functions, and $\nu$, $\beta$, $\gamma$ and $\alpha$ are the usual critical exponents.

From the Landau theory presented in Sec. IV, we expect that the transition belongs to the universality class of three dimensional $O(3)$ model with cubic anisotropy. In $O(N)$ models with cubic anisotropy, the phase transition is in the symmetric $O(N)$ universality class if $N<N_{c}$ and is in the cubic anisotropic universality class if $N>N_{c}$. Early work, using perturbative renormalisation group theory Ketley and Wallace (1973); Wallace (1973); Aharony (1973); Nelson et al. (1974); Brézin et al. (1974), high temperature series expansion Ferer et al. (1981) and non-perturbative RG calculations Tissier et al. (2002) suggests that $3<N_{c}<4$. Further work using RG calculations upto three loops Maier and Sokolov (1988); Shpot (1989), four loops Mayer et al. (1989); Varnashev (2000), five loops Pakhnin and Sokolov (2000); Kleinert et al. (1997); Kleinert and Schulte-Frohlinde (1995); Kleinert and Thoms (1995) find $N_{c}\lesssim 3$, while six loop RG calculations suggest $N_{c}\approx 2.89$. Monte Carlo simulations are consistent with $N_{c}=3$ Caselle and Hasenbusch (1998). However, in three dimensions, the exponents for the model with cubic anisotropic critical are very close to the exponents for the Heisenberg model Manuel Carmona et al. (2000). Therefore, we use the exponents for the three-dimensional Heisenberg model to analyse the data. In the following, we check that the data near the critical point are consistent with these exponents.

The critical point may be determined by the crossing point of the data for Binder cumulant for different system sizes. From this criterion, we find that the critical parameters are $\mu\approx 2.063$ corresponding to $\rho\approx 0.718$ [see Fig. 12(a)]. The data for Binder cumulant, susceptibility, and order parameter for different system sizes collapse onto one curve when scaled as in Eqs. (24)-(26) with the critical exponents for Heisenberg model in three dimensions, as shown in Fig. 12. We use the numerical estimates for the critical exponents $\nu=0.704$, $\beta=0.362$ and $\gamma=1.389$ Holm and Janke (1993). We conclude that the disordered-layered transition belongs to the universality class of the $O(3)$ model.

VI Layered-Sublattice Transition

In this section, we study the the nature of the layered-sublattice phase transition. We analyze the transition using the order parameter $q_{3}$, defined in Eq. (6), which is zero in the layered phase and non-zero in the sublattice phase.

Figures 13(a) and (b), show the time evolution of density and $q_{3}$ after equilibration. For clarity, we have also superimposed a running average of density, where each point has been averaged over 40 consecutive data points. Both $\rho$ and $q_{3}$ exhibit two states, one in which density is higher and $q_{3}$ is non-zero and another where density is lower and $q_{3}$ is approximately zero. The system fluctuates in time between these two states. This is characteristic of a first order transition where both the sublattice and layered phases have the same free energy at the transition point.

The probability distribution for density, $P(\rho)$, and $q_{3}$, $P(q_{3})$ for three different values of $\mu$, close to the transition point, are shown in Figs. 13(c) and (d) respectively. Note that we have used the smoothened density to obtain the distribution. As the transition point is crossed, it can be seen that the distribution changes from having more area at the lower density to having more area at the higher density. From the value of $\mu$ for which $P(\rho)$ has roughly same height, we conclude that the critical chemical potential is $\mu\approx 2.67$. Similar features may be seen for $P(q_{3})$. Finally, the effect of system size on the distributions at the critical activity are studied in Figs. 13(e) and (f). It can be seen that the peaks become higher and sharper with increasing system size. The jump in density at the transition is $\approx 0.001$. These are again characteristic of a first order transition, and we conclude that the layered-sublattice transition is discontinuous.

VII Sublattice-Columnar Transition

In this section, we study the the nature of the layered-sublattice phase transition. We analyze the transition using the order parameter $q_{3}$, as defined in Eq. (6). $q_{3}$ is zero in the columnar phase and non-zero in the sublattice phase.

Figures 14(a) and (b), show the time evolution of density and $q_{3}$ after equilibration. For clarity, we have also superimposed a running average of density, where each point has been averaged over 10 consecutive data points. Both $\rho$ and $q_{3}$ exhibit two states, one in which density is higher and $q_{3}$ is non-zero and another where density is lower and $q_{3}$ is approximately zero. The system fluctuates in time between these two states. This is characteristic of a first order transition or co-existence, where both the sublattice and columnar phases have the same free energy at the transition point. The jump in density across the transition is $\approx 0.0025$.

The probability distribution for density, $P(\rho)$, and $q_{3}$, $P(q_{3})$ for three different values of $\mu$, close to the transition point, are shown in Figs. 14(c) and (d) respectively. Note that we have used the smoothened density to obtain the distribution. As the transition point is crossed, it can be seen that the distribution changes from having more area at the lower density to having more area at the higher density. From the value of $\mu$ for which $P(\rho)$ has roughly same height, we conclude that the critical chemical potential is $\mu\approx 5.395$ for $L=60$. Similar features may be seen for $P(q_{3})$. Finally, the effect of system size on the distributions at the critical activity are studied in Figs. 14(e) and (f). We find that the critical point has a strong finite size dependence. For instance, if the critical density $\rho_{c}(L)$ is defined as the midpoint between the two peaks in the distribution $P(\rho)$, then we find $\rho_{c}(50)\approx 0.9522$, $\rho_{c}(60)\approx 0.9553$, and $\rho_{c}(70)\approx 0.9572$. Since the difference in $\rho_{c}(L)$ are of the order of the jump in $\rho$ (approximately $0.0025$), we plot the the probability distribution of $\Delta\rho$ where $\Delta\rho=\rho-\rho_{c}(L)$. From Fig. 14(e), we find that, with increasing system size, the peaks become higher and sharper. The same features are seen for $P(q_{3})$ [see Fig. 14(f)]. These are characteristics of a first order transition, and we conclude that the sublattice-columnar transition is weakly first order.

VIII Stability of Columnar Phase

In our Monte Carlo simulations, we are unable to equilibrate the system efficiently in the high density phase for densities larger than $\rho\approx 0.96$ for system sizes larger than $L\geq 100$. For these densities, we find that the system often gets stuck in very long-lived metastable states, which consist of layers of size $2\times L\times L$, each layer having a two dimensional columnar order. However, columnar order in consecutive layers may have different orientations. We illustrate with a typical example that was obtained in simulations. Consider a system that is layered in the $z$-direction. We define the columnar order parameter for layer $n$, $\mathbf{Q}_{z}(n)=Q_{z}(n)\mathrm{e}^{i\theta_{z}(n)}$, as

[TABLE]

where $\phi^{z}_{er}(n)$ and $\phi^{z}_{or}(n)$ are the packing fraction of the cubes with heads lying on even and odd rows of the $n$-th plane respectively, while $\phi^{z}_{ec}(n)$ and $\phi^{z}_{oc}(n)$ are the corresponding packing fractions on even and odd columns. For a layer with perfect columnar order, $\theta_{z}$ takes one of four values [math], $\pi/2$, $\pi$, $3\pi/2$. Similar definitions hold for $Q_{x}(n)$ and $Q_{y}(n)$. In Fig. 15, we show the variation of $Q_{z}(n)$ and $\theta_{z}(n)$ of a configuration that is layered in the $z$-direction (even planes are occupied), obtained after equilibrating for $10^{7}$ Monte Carlo steps. It can be see that while the magnitude remains constant across the even ayers, $\theta_{z}(n)=0$ or $\theta_{z}(n)=\pi/2$, showing that the columnar order in different planes have different orientations. We find that the system remains stuck in this meta stable phase for upto and beyond $10^{7}$ Monte Carlo steps.

Though such metastable states exist, we will show below that the true equilibrium state is one with the same columnar order in all the layers. The presence of cubes that common to adjacent layers tend to create an aligning interaction. This may be demonstrated though a perturbative calculation. This calculation is similar to the high density expansion developed for hard squares and rectangles Bellemans and Nigam (1967); Ramola and Dhar (2012); Nath et al. (2015); Nath and Rajesh (2014); Nath et al. (2016).

We start by assuming that the system is layered in the $z$-direction (even planes), and that $z$ is large enough so that there is perfect columnar order in each layer. We set up a perturbation expansion based on number of cubes with heads in odd planes. To do so, we introduce two activities, $z$ for cubes with heads on even planes and $z^{\prime}$ for cubes with heads on odd planes. When $z^{\prime}=0$, the layers are independent, and the problem reduces to that of $2\times 2$ hard square lattice gas model. We will determine, to first order in perturbation theory, the difference in free energy between a state in which all planes have even-row order and states in which one plane is misaligned with either odd-row order or column-ordered. We will denote these states by $S_{\alpha}$, whose partition function $\mathcal{L}^{\alpha}(z,z^{\prime})$ and the free-energy $F^{\alpha}(z,z^{\prime})=-\ln\mathcal{L}^{\alpha}(z,z^{\prime})$ can be formally written as

[TABLE]

For $S_{||}$, for which all planes have even-row order, the partition function, when there are no defect cubes (cubes with heads on odd planes), is

[TABLE]

where $\Omega_{p}(L)$ is the partition function of a periodic column of size $2\times 2\times L$. Consider a single defect cube that is placed in any of the $L/2$ odd planes. Wiithin a plane, it can choose any one of $L^{2}$ sites. The partition function for $S_{||}$ in the presence of one defect cube is

[TABLE]

where $\Omega_{o}(L-2)$ is the partition function of an open column of size $2\times 2\times(L-2)$. In Eq. (31), the first term represents the correction coming when the head of the defect cube is placed on an even row and second term corresponds to when the head of the defect cube is placed on an odd row.

Consider now a state $S_{ro}$ in which one of the planes (say $z=0$) is odd-row ordered. Since the partition function with no defect is identical to that for $S_{||}$, the difference in partition function appears only in the first-order correction term:

[TABLE]

The difference in free-energies, $\Delta F^{||,ro}(z,z^{\prime})$, may be written as

[TABLE]

Simplifying Eq. (33), we obtain

[TABLE]

Since the right hand side is a perfect square, $\Delta F^{||,ro}(z,z^{\prime})>0$ for any $L$. Thus, the state with one misaligned row-ordered state has higher free energy, and we conclude that introduction of defect cubes results in an effective aligning interaction that tends to make all the planes have columnar order with the same orientation.

The large $z$ behaviour of $\Delta F^{||,ro}(z,z^{\prime})$ may be determined by noting that for large $L$, $\Omega_{p}(L)=a_{p}\lambda^{L}$ and $\Omega_{o}(L-2)=a_{o}\lambda^{L-2}$ where $\lambda$ is the largest root of the equation $x^{2}-x-z=0$. From Ref. Nath et al. (2016), we obtain $a_{p}=1$, $a_{o}=\lambda/(2\lambda-1)$ and $\lambda=(1+\sqrt{1+4z})/2$. Setting $z^{\prime}=z$ and evaluating in the limit of $z\gg 1$, we obtain

[TABLE]

It is straightforward to generalize this calculation to the state $S_{ce}$ in which one of the planes has column-order. We omit the calculation, but we obtain a similar increase in free energy when a layer is misaligned.

We now ask whether metastable states, as seen in Fig. 15, are due to finite size effects or due to the algorithm being unable to equilibrate the system at high densities within available computer time. Though a misaligned plane results in a rise in free energy as Eq. (35), there is a gain in entropy $\ln 2$ per column when there are no defect cubes. Thus, we can identify a crossover length $L^{*}$ at which the free energy gained by alignment of a plane is balanced by the entropy lost due to alignment of such a plane. Equating the two free energies, $\Delta F^{||,ro}(z,z^{\prime})\sim\ln 2$, we obtain

[TABLE]

For the metastable state shown in Fig. 15, $z=403.43$ for the state , we obtain $L^{*}\approx 33$. For system sizes smaller than this length, the misaligned phases are favoured, but are a finite size effect. However, since the system lengths that we have simulated are much larger than $L^{*}$, we conclude that the presence of such metastable states are due to an inability of the Monte Carlo algorithm to equilibrate states with misaligned planes for large $z$, due to large entropic barriers.

IX Summary and conclusions

In this paper we studied the phases and the phase transitions in a system of $2\times 2\times 2$ hard cubes on a three dimensional cubic lattice. We show the existence of four different phases. In order of increasing density, these are a disordered phase, a layered phase in which the system breaks up into $L/2$ interacting slabs of size $2\times L\times L$ each having fluid-like order, a solid-like sublattice phase, where the cubes preferentially occupy one sublattice, and a columnar phase in which the system breaks up into $L^{2}/4$ columns of dimension $2\times 2\times L$ with a fluid-like order within a column. The disordered-layered transition is shown to be a continuous transition that is consistent with the universality class of the three dimensional $O(3)$ model with cubic anisotropy. The other two transitions – layered-sublattice and sublattice-columnar – are shown to be discontinuous. The phase diagram is summarized in Fig. 16.

We formulated a Landau theory, consistent with the symmetries of the system, that is able to describe all the phases seen in the Monte Carlo simulations. Within the minimal functional, as described in Eq. (7), based only on the layering vector $\mathbf{L}$, we find that the columnar phase is unstable for the full range of the parameters, while it predicts the existence and stability of disordered, layered and sublattice phases. It also predicts the disordered-layered and disordered-sublattice transitions to be continuous and belong to the universality class of the three-dimensional $O(3)$ model with cubic anisotropy. The layered-sublattice transition is predicted to be first order. To obtain a stable columnar phase, we extended the free energy functional to explicitly depend on the columnar vector $\mathbf{C}$ [see Eq. (19)]. Within this extended functional, it is possible to show the existence and stability of all the different phases observed in simulations.

The results in this paper are not consistent with the theoretical predictions of density functional theory. Density functional theory predicts that $2\times 2\times 2$ hard cubes cubes undergo transitions from a disordered phase to layered phase to a columnar phase at high densities Lafuente and Cuesta (2003). However, it does not predict the sublattice phase, that is seen in our Monte Carlo simulations. Understanding why the theory fails, and how it should be modified to give the correct predictions is a promising area for future study. For $6\times 6\times 6$ cubes, the theory predicts a transition from a disordered to solid to two types of columnar phase Lafuente and Cuesta (2003). Testing these predictions in simulations would also be of interest.

The results in this paper are also in contradiction to earlier Monte Carlo simulations Panagiotopoulos (2005), wherein no phase transitions were found even though the simulations were performed close to full packing (in Ref. Panagiotopoulos (2005), the problem of cubes correspond to $\sigma=2$). This discrepancy could be due to small system sizes that were studied ($L=18,24$) in Ref. Panagiotopoulos (2005) compared to the systems sizes studied in the current paper ($L$ upto $200$).

The existence of a sublattice phase is quite surprising. It would appear that as you introduce vacancies at full packing, the sublattice phase gets destabilised in favour of the columnar phase. However, a larger number of vacancies somehow stabilises the sublattice phase. Also, it implies that the interactions between the different layers in the layered phase are not weak. If they were weak, then we would expect that once the system becomes layered, the problem becomes effectively a problem of hard squares in two dimensions. This lower dimensional system does not exhibit a sublattice phase.

In the continuum, the system of parallel hard cubes undergoes a continuous freezing transition from a disordered fluid phase to a solid phase as density is increased Jagla (1998); Groh and Mulder (2001), consistent with theoretical predictions using density functional theory Belli et al. (2012). The continuum limit may be reached by determining the phase diagram for $k\times k\times k$ cubes and extrapolating for large $k$. Preliminary simulations for $3\times 3\times 3$ cubes suggest that the sublattice phase does not exist, but the layered and columnar phases exist. Thus, for larger $k$, the layered-columnar transition would appear to be similar to the disordered-columnar transition in $k\times k$ hard squares. For this model, high density expansions suggest that the critical density tends to an asymptotic value that is less than one for large $k$ Nath et al. (2015). If this were true, the continuum problem should have two transitions. Re-examining the problem of parallel hard cubes in the continuum is a promising area for future study.

Acknowledgements.

The simulations were carried out in single node cluster machine Nandadevi (2 x Intel Xeon E5-2667 3.3 GHz) using OpenMP parallelisation. DD’s research was partially supported by the J. C. Fellowship, awarded by the Department of Science and Technology, India, under the grant DST-SR-S2/JCB-24/2005.