Phase Transitions in a System of Hard Rectangles on the Square Lattice

# Phase Transitions in a System of Hard Rectangles on the Square Lattice

## Abstract

The phase diagram of a system of monodispersed hard rectangles of size on a square lattice is numerically determined for and aspect ratio . We show the existence of a disordered phase, a nematic phase with orientational order, a columnar phase with orientational and partial translational order, and a solid-like phase with sublattice order, but no orientational order. The asymptotic behavior of the phase boundaries for large are determined using a combination of entropic arguments and a Bethe approximation. This allows us to generalize the phase diagram to larger and , showing that for , the system undergoes three entropy driven phase transitions with increasing density. The nature of the different phase transitions are established and the critical exponents for the continuous transitions are determined using finite size scaling.

###### pacs:
64.60.De, 64.60.Cn, 05.50.+q

## I Introduction

The study of entropy-driven phase transitions in a system of long hard rods has a long history dating back to Onsager’s demonstration Onsager (1949) that the three dimensional system undergoes a transition from an isotropic phase to an orientationally ordered nematic phase as the density of the rods is increased Onsager (1949); Flory (1956); de Gennes and Prost (1995); Vroege and Lekkerkerker (1992); Dhar et al. (2011). Further increase in density may result in a smectic phase that partially breaks translational symmetry, and a solid phase Frenkel et al. (1988); Bolhuis and Frenkel (1997). In two dimensional continuum space, the continuous rotational symmetry remains unbroken but the system undergoes a Kosterlitz-Thouless type transition from a low-density phase with exponential decay of orientational correlations to a high-density phase having quasi long range order Straley (1971); Frenkel and Eppenga (1985); Wojciechowski and Frenkel (2004); Khandkar and Barma (2005); Donev et al. (2006); Zhao et al. (2007); Vink (2009). Experimental realizations include tobacco mosaic virus Fraden et al. (1989), liquid crystals de Gennes and Prost (1995), carbon nanotube gels Islam et al. (2004), and brownian squares Zhao et al. (2011). The phenomenology is, however, much less clear when the orientations are discrete, and the positions are either on a lattice or in the continuum, when even the existence of the nematic phase has been convincingly seen in simulations only recently Ghosh and Dhar (2007); Disertori and Giuliani (2013).

Consider hard rectangles of size on a two-dimensional square lattice where each rectangle occupies () lattice sites along the short (long) axis. The limiting cases when either the aspect ratio or are better studied. When and , there are, remarkably, two entropy-driven transitions: from a low-density isotropic phase to an intermediate density nematic phase, and from the nematic phase to a high-density disordered phase Ghosh and Dhar (2007); Kundu et al. (2013). While the first transition is in the Ising universality class Matoz-Fernandez et al. (2008a); Fischer and Vink (2009), the second transition could be non-Ising Kundu et al. (2013), and it is not yet understood whether the high density phase is a re-entrant low density phase or a new phase Kundu et al. (2013); Kundu and Rajesh (2013). When (hard squares), the system undergoes a transition into a high density columnar phase. The transition is continuous for  Baxter (1982); Fernandes et al. (2007); Zhitomirsky and Tsunetsugu (2007); Feng et al. (2011); Ramola (2012), and first order for  Fernandes et al. (2007). When , keeping fixed, the lattice model is equivalent to the model of oriented rectangles in two dimensional continuum, also known as the Zwanzig model Zwanzig (1963). For oriented lines in the continuum (), a nematic phase exists at high density Fischer and Vink (2009). The only theoretical results that exist are when and (dimers), for which no nematic phase exists Heilmann and Lieb (1970); Kunz (1970); Gruber and Kunz (1971); Heilmann and Lieb (1972), , when the existence of the nematic phase may be proved rigorously Disertori and Giuliani (2013), and an exact solution for arbitrary on a tree like lattice Dhar et al. (2011); Kundu and Rajesh (2013).

Less is known for other values of and . Simulations of rectangles of size did not detect any phase transition with increasing density Barnes et al. (2009a), while those of parallelepipeds on cubic lattice show layered and columnar phases, but no nematic phase Casey and Harrowell (1995). In general, numerical studies of large rectangles are constrained by the fact that it is difficult to equilibrate the system at high densities using Monte Carlo algorithms with local moves, as the system gets jammed and requires correlated moves of several particles to access different configurations.

In addition to being the lattice version of the hard rods problem, the study of lattice models of hard rectangles is useful in understanding the phase transitions in adsorbed monolayers on crystal surfaces. The and planes of fcc crystals have square and rectangular symmetry and may be treated with lattice statistics if the adsorbate-adsorbate interaction is negligible with respect to the periodic variation of the corrugation potential of the underlying substrate Patrykiejew et al. (2000). For example, the critical behavior of a monolayer of chlorine (Cl) on Ag is well reproduced by the hard square model () and the high density structure of the Cl-adlayer may be mapped to the high density phase of the hard square problem Taylor et al. (1985). Structures such as , and are ordered structures. Lattice gas model with repulsive interaction up to forth nearest neighbor has been used to study the phase behavior of selenium adsorbed on Ni Bak et al. (1985). The results in this paper show that different phases with orientational and positional order may be obtained by only hard core exclusion.

At a more qualitative level, discrete models of hard rods have been used to obtain realistic phase diagram for polydispersed systems Clarke et al. (2000); Martínez-Ratón and Cuesta (2003), to describe orientational wetting of rods van Roij et al. (2000), to model and understand self-assembly of nano particles on monolayers Barnes et al. (2009a) and thermodynamics of linear adsorbates Ramirez-Pastor et al. (1999); Romá et al. (2003).

Lattice models of hard rectangles also falls in the general class of hard core lattice gas models of differently shaped particles. The study of these models also have been of continued interest in Statistical Physics as minimal models for the melting transitions. The different shapes studied include squares Bellemans and Nigam (1967); Pearce and Seaton (1988); Ramola and Dhar (2012), hexagons on triangular Baxter (1980); Heilmann and Praestgaard (1973) and square lattices Dickman (2012), triangles Verberkmoes and Nienhuis (1999), and tetrominoes Barnes et al. (2009a).

In this paper, we adapt and implement an efficient Monte Carlo algorithm with cluster moves that was very effective in studying the hard rod () problem on lattices Kundu et al. (2013, 2012). The hard rectangle model and the algorithm are described in detail in Sec. II. We observe four distinct phases at different densities: isotropic, nematic, columnar, and sublattice phases. These phases, suitable order parameters to characterize them, and other thermodynamic quantities are defined in Sec. III. From extensive large scale simulations, we determine the rich phase diagram for , and . The phase diagram for is discussed in Sec. IV. We find that all transitions except the isotropic-columnar transition for are continuous. The critical exponents and universality classes of the continuous transitions are determined. Section V contains the details about the phase diagram and the nature of the phase transitions for . In Sec. VI, we use a Bethe approximation and estimates of entropies for the different phases to determine the phase boundaries for large , allowing us to generalize the phase diagram to arbitrary and . In particular, it allows us to take the continuum limit , thus obtaining predictions for a system of oriented rectangles in the continuum. Sec. VII contains a summary and a discussion of possible extensions of the problem.

## Ii The Model and Monte Carlo algorithm

We define the model on a square lattice of size with periodic boundary conditions. Consider a mono-dispersed system of rectangles of size such that the aspect ratio is . A rectangle can be either horizontal or vertical. A horizontal (vertical) rectangle occupies sites along the ()-axis and sites along the ()-axis. No two rectangles may overlap. An activity is associated with each rectangle, where is the chemical potential.

We study this system using constant grand canonical Monte Carlo simulations. The algorithm that we implement is an adaptation of the algorithm with cluster moves that was introduced in Refs. Kundu et al. (2013, 2012) to study the problem of hard rods (m=1). We describe the algorithm below. Given a valid configuration of rectangles, in a single move, a row or a column is chosen at random. If a row is chosen, then all horizontal rectangles whose heads (bottom, left corner) lie in that row are removed, leaving the rectangles with heads in other rows untouched. The emptied row now consists of two kinds of sites: forbidden sites that cannot be occupied with horizontal rectangles due to the presence of vertical rectangles in the same row or due to rectangles with heads in the neighboring rows, and sites that may be occupied by horizontal rectangles in a valid configuration. An example illustrating the forbidden sites is shown in Fig. 1(a). It is clear that the sites that may be occupied are divided into intervals of contiguous empty sites. The problem of occupation of the emptied row with a new configuration now reduces to the problem of occupying the empty intervals. However, the empty intervals may be occupied independent of each other, as the occupation of one is not affected by the configuration of rectangles in the remaining ones. Thus, the re-occupation of the emptied row reduces to a problem of occupying a one dimensional interval with rods. This problem is easily solvable and the equilibrium probabilities of each new configuration may be easily calculated. We refer to Refs. Kundu et al. (2013, 2012) for the calculation of these probabilities. If a column is chosen instead of a row, then a similar operation is performed for the vertical rectangles whose heads lie in that column.

In addition to the above evaporation-deposition move, we find that the autocorrelation time is reduced considerably by introducing a flip move. In this move, a site is picked at random. If it is occupied by the head of a horizontal rectangle, then we check whether sites are occupied by the heads of horizontal rectangles. If that is the case, we call this set of aligned rectangles a rotatable plaquette of horizontal rectangles. In the flip move such a rotatable plaquette of size , containing horizontal rectangles, is replaced by a similar plaquette of vertical rectangles. An example of the flip move is shown in Fig. 1(b). If is occupied by the head of a vertical rectangle and a rotatable plaquette of vertical rectangles is present, then it is replaced by a plaquette of aligned horizontal rectangles. A Monte Carlo move corresponds to evaporation–deposition moves and flip moves. It is easy to check that the algorithm is ergodic and obeys detailed balance.

We implement a parallelized version of the above algorithm. In the evaporation-deposition move, we simultaneously update all rows that are separated by . Once all rows are updated in this manner, the columns are updated. We also parallelize the flip move. The lattice is divided into blocks of size . The flipping of each of these blocks is independent of the other and may therefore be flipped simultaneously. We flip a rotatable plaquette with probability . The parallelization and efficiency of the algorithm allows us to simulate large systems (up to ) at high densities (up to ).

We check for equilibration by starting the simulations with two different initial configurations and making sure that the final equilibrium state is independent of the initial condition. One configuration is a fully nematic state, where all rectangles are either horizontal or vertical and the other is a random configuration where rectangles of both vertical and horizontal orientations are deposited at random.

## Iii The Different Phases

Snapshots of the different phases that we observe in simulations are shown in Fig. 2. First is the low density isotropic (I) phase in which the rectangles have neither orientational nor translational order [see Fig. 2(c)]. Second is the nematic (N) phase in which the rectangles have orientational order but no translational order [see Fig. 2(d)] . In this phase, the mean number of horizontal rectangles is different from that of vertical rectangles. The third phase is the columnar (C) phase, having orientational order and partial translational order [see Fig. 2(e)]. In this phase, if majority of rectangles are horizontal (vertical), then their heads, or bottom left corners, preferably lie in rows (columns) that are separated by . Thus, it breaks the translational symmetry in the direction perpendicular to the orientation but not parallel to the orientation. Clearly, there are symmetric C phases. In this phase the rectangles can slide much more along one lattice direction. The fourth phase is the crystalline sublattice (S) phase with no orientational order [see Fig. 2(f)]. We divide the square lattice into sublattices by assigning to a site a label . The sublattice labeling for the case is shown in Fig. 2(a). In the S phase, the heads of the rectangles preferably occupy one sublattice, breaking translational symmetry in both the directions.

From the symmetry of the system we would expect up to phases. The orientational symmetry could be present or broken while the translational symmetry could be unbroken, broken along only one or both - and - directions. If the orientational symmetry is broken, then the translational symmetry could be broken either parallel or perpendicular to the preferred orientations. Out of the possibilities, we do not observe (i) a phase with no orientational order but partial translational order, (ii) a phase with orientational order and complete translational order (iii) a smectic like phase in which orientational order is present and translational symmetry parallel to the orientation is broken.

To distinguish among the four different phases we define the following order parameters:

 Q1 = m2k⟨Nh−Nv⟩N, (1a) Q2 = m2k⟨|∑m2−1j=0nje2πijm2|⟩N, (1b) Q3 = m2k⟨|∑m−1j=0rje2πijm|−|∑m−1j=0cje2πijm|⟩N, (1c) Q4 = m2k⟨n0−n1−n2+n3⟩N, (1d)

where and are the total number of horizontal and vertical rectangles respectively, is the number of rectangles whose heads are in sublattice , are the number of rectangles whose heads are in row , and are the number of rectangles whose heads are in column . All four order parameters are zero in the I phase. is non-zero in the N and C phases, is non-zero in the C and S phases, is non-zero only in the C phase, and is non-zero only in the S phase. in Eq. (1d) has been defined for . Its generalization to is straightforward.

We now define the thermodynamic quantities that are useful to characterize the transitions between the different phases. ’s second moment , compressibility and the Binder cumulant are defined as

 χi = ⟨Q2i⟩L2, (2a) κ = [⟨ρ2⟩−⟨ρ⟩2]L2, (2b) Ui = 1−⟨Q4i⟩3⟨Q2i⟩2. (2c)

The transitions are accompanied by the singular behavior of the above thermodynamic quantities at the corresponding critical densities. Let , where is the critical chemical potential. The singular behavior is characterized by the critical exponents , , , defined by , , , , and , where is the correlation length, , and represents any of the order parameters. Only two exponents are independent, others being related to them through scaling relations.

The critical exponents , , and are obtained by finite size scaling of the different quantities near the critical point:

 U ≃ fu(ϵL1/ν), (3a) Q ≃ L−β/νfq(ϵL1/ν), (3b) χ ≃ Lγ/νfχ(ϵL1/ν), (3c) κ ≃ Lα/νfκ(ϵL1/ν), (3d)

where , , , and are scaling functions.

## Iv Phase Diagram and Critical Behavior for m=2

In this section, we discuss the phase diagram for the case and aspect ratio . The critical exponents characterizing the different continuous transitions are determined numerically.

### iv.1 Phase Diagram

The phase diagram obtained from simulations for and integer are shown in Fig. 3. The low density phase is an I phase for all . The case is different from other . In this case, the problem reduces to a hard square problem and orientational order is not possible as there is no distinction between horizontal and vertical rectangles. The hard square system undergoes only one transition with increasing density, it being a continuous transition from the I phase to a C phase Baxter (1982); Bellemans and Nigam (1967); Pearce and Seaton (1988); Ramola and Dhar (2012). This transition belongs to the Ashkin Teller universality class (see Refs. Fernandes et al. (2007); Zhitomirsky and Tsunetsugu (2007); Feng et al. (2011); Ramola (2012) for recent numerical studies). For , we find that the system undergoes one continuous transition directly from the I phase to a crystalline S phase. On the other hand, the system with may exist in I, C, or S phases. With increasing density, the system undergoes two phase transitions: first from the I to a C phase which could be continuous or first order, and second, from the C to a S phase which is continuous. For , we observe three continuous transitions with increasing density: first from the I to the N phase, second into the C phase and third into the S phase. By confirming the existence of the N and C phases for , we expect the phase behavior for to be similar to that for .

The system undergoes more than one transition only for . We now present some supporting evidence for this claim. In Fig. 4 (a) we show the probability distribution of the nematic order parameter , when , for different values of and fixed , close to the I-C and the C-S transitions. For lower values of , the distribution is peaked around zero corresponding to the I phase. With increasing , the distribution becomes flat and two symmetric maxima appear at ( also becomes nonzero simultaneously), corresponding to a C phase. On increasing further, the two maxima continuously merge into a single peak at =0, corresponding to the S phase ( also becomes zero and becomes nonzero). Fig. 4 (b) shows the distribution of for three different system sizes at a fixed value of for which has two symmetric maxima at . The two peaks become sharper and narrower with increasing . We find the similar behavior for also. From the above, we conclude that the C phase exists for albeit for a very narrow range of . For and we do not observe the existence of a columnar phase and find that the probability distributions of and are peaked around zero for all . Hence, we conclude that .

The N phase exists only for . This is also true for  Ghosh and Dhar (2007). To see this, notice that the I-C transition for is first order (see Fig. 3). If a nematic phase exists for , then the first transition would have been continuous and in the Ising universality class Matoz-Fernandez et al. (2008a).

### iv.2 Critical behavior for the isotropic-sublattice (I-S) phase transition

The system of rectangles with undergoes a direct I-S transition for . At this transition, the translational symmetry gets broken along both - and - direction but the rotational symmetry remains preserved. We study this transition using the order parameter [see Eq. (1b)]. is non-zero in the S phase and zero in the I phase. In this case and remains zero for all values of . The data collapse of , , and for different values of near the I-S transition are shown in Fig. 5 for and in Fig. 6 for . From the crossing of the Binder cumulant data for different , we estimate the critical chemical potential () for and () for . The order parameter increases continuously with from zero as is crossed, making the transition continuous. Since the phase has a four fold symmetry due to the four possible sublattices, we expect the transition to be in the Ashkin-Teller universality class. Indeed, we find a good collapse with and . Numerically, we find for and for . being larger than , we do not observe any divergence in . This transition could have also been studied using the order parameter .

### iv.3 Critical behavior of the isotropic-columnar phase (I-C) transition

The I-C transition is seen for . When for , for the I-C and the C-S transitions are close to each other, making unsuitable for studying the critical behavior. We, therefore, study the I-C transition for ( rectangles) and ( rectangles).

The critical behavior is best studied using the order parameter [see Eq. (1c)]. is non-zero only in the C phase. First, we present the critical behavior for . The simulation data for different system sizes are shown in Fig. 7. From the crossing of the Binder cumulant curves, we obtain (). The transition is found to be continuous. There are four possible columnar states: majority of heads are either in even or odd rows (when horizontal orientation is preferred), or in even or odd columns (when vertical orientation is preferred). Due to this four fold symmetry, we expect the I-C transition to be in the Ashkin-Teller universality class. The data for different collapse with , and (see Fig. 7), confirming the same. Unlike the I-S transition for , and lies between the Ising and q=4 Potts points. At the I-C transition partial breaking of translational symmetry and complete breaking of rotational symmetry occur simultaneously.

The I-C transition for is also continuous [see Fig. 4 (a)], and is therefore expected to be in the Ashkin-Teller universality class. However, there is no reason to expect that will be the same as that for .

For , the I-C transition is surprisingly first order. Fig. 8(a) shows the time profile of density near the I-C transition. alternates between two well defined densities, one corresponding to the I phase and the other to the C phase. This is also seen in the probability distribution for density [see Fig. 8(b)]. Near the I-C transition it shows two peaks corresponding to the I and the C phases. Thus, at , the density has a discontinuity, which is shown by the shaded region in the phase diagram (see Fig. 3). The probability distribution of the order parameter shows similar behavior [see Fig. 8(c)]. Near the distribution shows three peaks: one at corresponding to the I phase and the other two at , corresponding to the C phase. At the three peaks become of equal height and the order parameter jumps from zero to a non-zero value. These peaks sharpen with increasing system size [see Fig. 8(d)]. These are typical signatures of a first order transition. Hence, we conclude that the I-C transition may be continuous or first order depending on .

### iv.4 Critical behavior of the isotropic-nematic phase (I-N) transition

We find that the nematic phase exists only for . We study the I-N phase transition for using the order parameter [see Eq. (1a)]. is non-zero in the N and C phases and zero in the I phase. We confirm that the ordered phase is an N phase by checking that , which is non-zero only in the C phase, is zero. In the nematic phase the rectangles may choose either horizontal or vertical orientation. Thus, we expect the transition to be in the Ising universality class. When , this has been verified using extensive Monte Carlo simulations Matoz-Fernandez et al. (2008a). Here, we confirm the same for . The data for , and for different collapse onto one curve when scaled with the two dimensional Ising exponents , , and (see Fig. 9). We find (). We note that the value of at the point where the curves for different cross is slightly smaller than the Ising value . This suggests that larger system sizes are necessary for better collapse of the data.

### iv.5 Critical behavior of the nematic-columnar phase (N-C) transition

The N-C transition is also studied for , using the order parameter . is zero in the nematic phase but nonzero in the columnar phase. At the I-N transition orientational symmetry gets broken. If the nematic phase consists of mostly horizontal (vertical) rectangles, then there is no preference over even and odd rows (columns). In the columnar phase the system chooses either even or odd rows (columns), once the orientational symmetry is broken. Due to the two broken symmetry phases we expect this transition to be in the Ising universality class. We indeed find good data collapse when , and for different system sizes are scaled with Ising exponents (see Fig. 10). The critical chemical potential or critical density is obtained from the crossing point of the binder cumulant for different . We find () for this transition. We expect the critical behavior to be same for .

### iv.6 Critical behavior of the columnar-sublattice phase (C-S) transition

The C-S transition exists for . This transitions is studied by choosing . We characterize the C-S transition using the order parameter which is non zero only in the S phase. In the C phase the system chooses one particular orientation and either even or odd rows or columns, depending on the orientation. This corresponds to two sublattices being chosen among four of them. In the C-S transition the translational symmetry gets broken completely by choosing a particular sublattice, but along with that the orientational symmetry gets restored. This transition is found to be continuous. The data of , and for different near the C-S transition collapse well when scaled with the exponents belonging to the Ashkin-Teller universality class. The estimated critical exponents are , and (see Fig. 11). Binder cumulants for different system sizes cross at (). We expect similar behavior for and but possibly with different . The C-S transition occurs at very high density. With increasing , the relaxation time becomes increasingly large, making it difficult to obtain reliable data for the C-S transition when .

## V Phase Diagram and Critical Behavior for m=3

### v.1 Phase diagram

The phase diagram that we obtain for , is shown in Fig. 12. When , the corresponding hard square system has a single, first order transition from the I phase into the C phase Fernandes et al. (2007). The shaded region between two points denotes a region of phase coexistence. For , the system undergoes two first order transitions with increasing density: first an I-C transition and second a C-S transition. This is unlike the case , where for and we find only one transition. For , we find three transitions as in the case. The first transition from I to N phase is continuous while the second from N to C phase appears to be first order. Although we cannot obtain reliable data for the third transition into the S phase, we expect it to be first order. We note that the minimum value of beyond which the nematic phase exists is for both and , and matches with that for  Ghosh and Dhar (2007).

### v.2 The isotropic-columnar phase (I-C) transition

The I-C transition exists when . We study the this transition for , using the order parameter . Now there are six possible choices for the C phase: heads are predominantly in one of the rows or (mod ) with all the columns equally occupied (if horizontal orientation is preferred) or in one of the columns or (mod ) and all the rows are equally occupied (if vertical orientation is preferred). Making an analogy with the six state Potts model, we expect the I-C transition to be first order. The probability distribution of the density and the order parameter for near the I-C transition is shown in Fig. 13. The distributions are clearly double peaked at and near the transition point, one corresponding to the I phase and the other to the C phase. We find that these peaks become sharper with increasing system size. This is suggestive of a first order phase transition with a discontinuity in both density and order parameter as crosses . The discontinuity in the density is denoted by the shaded regions of Fig. 12. The chemical potential at which the I-C transition occurs is given by, . Similar behavior is seen near the I-C transition for rectangles of size with and . We observe that the discontinuity in the density increases with .

### v.3 Critical behavior of the isotropic-nematic phase (I-N) transition

As for , for we find the existence of the nematic phase only for . We study the I-N transition for with the order parameter . It is expected to be in the Ising universality class since there are two possible choices of the orientation: either horizontal or vertical. We are unable to obtain good data collapse for , and as the relaxation time increases with increasing and . Instead, we present some evidence for the transition being continuous and belonging to the Ising universality class. In Fig. 14 (a), the distribution of the order parameter near the I-N transition is shown. The two symmetric peaks of the distribution come closer with decreasing and merge to a single peak, this being a signature of a continuous transition. The Binder cumulant for different system sizes cross at () [see Fig. 14 (b)]. The value of at is very close to the value () for the Ising universality class.

### v.4 The nematic-columnar phase (N-C) transition

The N-C transition is studied for using the order parameter . Contrary to our expectation that the N-C transition should be in the q=3 Potts universality class, we observe a first order transition. The temporal dependence of the density near the N-C transition is shown in Fig. 15(a). Density jumps between two well separated values corresponding to the two different phases near the coexistence. Fig 15(b) shows the discontinuity in the order parameter near the transition. shows two peaks of approximately equal height near . However, we are limited in our ability to obtain reliable data for rectangles for larger system sizes, and the observed first order nature could be spurious.

### v.5 The columnar-sublattice phase (C-S) transition

The C-S transition is studied by choosing . We use the order parameter which is nonzero only in the S phase. The probability distribution of the density and the order parameter for rectangles near the C-S transitions is shown in Fig. 16. The distributions are again double peaked at and near the transition point, making the C-S transition first order. These peaks become sharper with increasing system size. The discontinuity in the density near the C-S transition is very small and can also be seen in the shaded portions of Fig. 12. We estimate . Similar behavior near the C-S transitions is also observed for , but the relaxation time increases with .

## Vi Estimation of the phase boundaries using analytical methods

In this section we obtain the asymptotic behavior of the phase diagram for large using theoretical arguments.

### vi.1 The isotropic–nematic phase boundary

The critical density for the I-N phase transition, for fixed and may be determined by making an analogy with the continuum problem. The limit , keeping fixed corresponds to the system of oriented lines in the continuum. For this problem  Ghosh and Dhar (2007); Matoz-Fernandez et al. (2008b). The constant is estimated to be  Ghosh and Dhar (2007); Matoz-Fernandez et al. (2008b). Thus, we expect , where is independent of . For , we observe only a weak dependence of on with the critical density being , and .

### vi.2 The nematic–columnar phase boundary

To obtain the asymptotic behavior of the N-C phase boundary, we use an ad hoc Bethe approximation scheme for rods due to DiMarzio DiMarzio (1961), adapted to other shapes Sokolova and Tumanyan (2000). To estimate the phase boundary of the nematic–columnar transition of rectangles on the square lattice with sites, we require the entropy as a function of the occupation densities of the types of rows/columns. The calculations become much simpler, if we consider a fully oriented phase with only horizontal rectangles. Now, the nematic phase corresponds to the the phase where there is equal occupancy of each of the types of rows, while the columnar phase breaks this symmetry and preferentially occupies one type of row. For this simplified model with only one orientation, we estimate the entropy within an ad hoc Bethe approximation as detailed below. We present the calculation for , classifying the rows as even and odd rows. Generalization to higher values of is straight forward.

Let there be () number of rectangles whose heads (left bottom site of the rectangle) occupy even (odd) rows. We first place the rectangles one by one on the even rows. Given that rectangles have already been placed, the number of ways in which the rectangle can be placed may be estimated as follows. The head of the rectangle has to be placed on an empty site of an even row. We denote this site by (see Fig. 17). The site can be chosen at random in () ways, being the number of sites in even rows and being the number of occupied sites in the even rows by the rectangles. We now require that the consecutive sites to the right of are also empty. The probability of this being true is , where is the conditional probability that (see Fig. 17) is empty given that is empty. In terms of and , is given by

 Px(B|A)=M2−2kjeM2−2kje+je. (4)

To place the rectangle, we also require that the site (see Fig. 17) and the consecutive sites to the right of are also empty. The probability of this being true is , where is the conditional probability that is empty given is empty, and is the conditional probability that (see Fig. 17) is empty given that both and are empty. Sites C and D belongs to an odd row. Since both A and B are empty, C and D can be occupied only by rectangles with heads in the same odd row. But, there are no such rectangles. Therefore, , and . Thus, given that rectangles have been placed, the rectangle may be placed in

 νje+1=(M2−2kje)×[Px(B|A)]2k−1 (5)

ways. Hence, the total number of ways of placing rectangles with heads on even rows is

 Ωe = 1Ne!Ne−1∏je=0νje+1 (6) = 1Ne!Ne−1∏je=0(M2−2kje)2k(M2−2kje+je)2k−1.

Keeping the rectangles with heads on even rows, we now place rectangles one by one on the odd rows. Given that rectangles have been placed on the odd rows, the number of ways of placing the rectangle may be estimated as follows. The head of the rectangle must be placed on an empty site on an odd row. We denote this site by (see Fig. 17). may be chosen in ways, where we have ignored correlations between rectangles. Here is the number of occupied sites in the odd rows due to the rectangles on the even rows, and the is the number of sites occupied by rectangles in odd rows. We now require that the consecutive sites to the right of are also empty. The probability of this being true is , where is the conditional probability that (see Fig. 17) is empty given that is empty. is given by

 Px(F|E)=M/2−2kNe−2kjoM/2−2kNe−2kjo+Ne+jo, (7)

where we have again ignored all correlations.

For placing the rectangle, we also require that the site (see Fig. 17) and the consecutive sites to the right of are also empty. The probability of this being true is , where is the conditional probability that is empty given is empty, and is the conditional probability that (see Fig. 17) is empty given that both and are empty. Ignoring correlations, is given by

 Py(G|E)=M/2−2kNe−2kjoM/2−2kjo. (8)

If we calculate following the procedure developed by DiMarzio in Ref. DiMarzio (1961), then the resultant entropy is not symmetric with respect to and and depends on the order of placement. To overcome this shortcoming, we follow the Bethe approximation proposed in Ref. Sokolova and Tumanyan (2000) as follows:

 P(H|F∩G) = P(F∩G|H)P(H)Pxy(G|F)P(F), (9) ≈ Px(G|H)Py(F|H)Pxy(G|F). (10)

where in Eq. (9), we used and in Eq. (10), we replaced by , which is an approximation.

In Eq.(10), from symmetry, it is easy to see that and and can be read off from Eqs. (7) and (8). To obtain an expression for , the probability that is empty, given that the site is empty, we again ignore correlations. We then obtain

 Pxy(G|F)=M/2−2kNe−2kjoM/2−2kjo+jo. (11)

The number of ways of placing the rectangle is

 νjo+1 = (M2−2kNe−2kjo)×Px(F|E)]2k−1 (12) × Py(G|E)×[P(H|G∩F)]2k−1.

Substituting for each of the quantities on the right hand side, we obtain the total number of ways of placing the rectangles on the odd rows as

 Ωo = 1No!No−1∏jo=0νjo+1 (13) = 1No!No−1∏jo=0(M2−2kNe−2kjo)4k(M2−2kjo)2k ×(M2−2kjo+jo)2k−1(M2−2kNe−2kjo+Ne+jo)4k−2.

We would like to express the entropy in terms of the total density and the densities of occupied sites in even and odd rows, given by and respectively. Clearly,

 ρe = 4kNeM, (14) ρo = 4kNoM, (15) ρ = ρe+ρo. (16)

The entropy per site in the thermodynamic limit is given by

 s(ρe,ρo)=limM→∞1Mln(Ω0Ωe). (17)

Substituting for and from Eqs. (6) and (13), we obtain

 s(ρe,ρo) = −∑i=o,eρi4klnρi2k−(1−ρ)ln(1−ρ)+(1−ρ+ρ2k)ln(1−ρ+ρ2k) (18) +∑i=o,e1−ρi2ln(1−ρi)−12∑i=o,e(1−ρi+ρi2k)ln(1−ρi+ρi2k).

We express the entropy in terms of density and the order parameter , defined as

 ψN−C=ρe−ρoρ. (19)

is zero in the nematic phase and non-zero in the columnar phase. For a fixed value of , the equilibrium values of and are determined by maximizing the entropy with respect to . In Fig. 18 we show the variation of entropy with for different densities. For small values of , the entropy is maximized by , , . Beyond a critical density , is maximized by , , . grows continuously with for , and thus the transition for is continuous.

The expansion of in powers of has only even powers of since is invariant when . Thus, the critical density is obtained from the condition . This gives

 ρN−Cc = −1+4k−√1−4k+8k22k−1 (20) = (2−√2)+Ak+O(k−2),

where . We note that as , tends to a independent value and that the transition exists for all .

We can similarly calculate the entropy for and then generalize the expression of entropy for arbitrary and . Now, there are densities , , , , corresponding to the types of rows. In terms of them, the entropy is given by

 s({ρi}) = −m∑i=1ρim2klnρimk−(1−ρ)ln(1−ρ)+(1−ρ+ρmk)ln(1−ρ+ρmk) (21) +m∑i=11−ρ+ρimln(1−ρ+ρi)−1mm∑i=1(1−ρ+ρi+ρ−ρimk)ln(1−ρ+ρi+ρ−ρimk).

Here, we define the order parameter to be

 ψN−C=ρ1−ρ2ρ, (22)

where we set . Now, is not invariant when . Thus, when expanded in powers of , has cubic terms, making the transition first order. This is illustrated in Fig. 19 which shows the variation of entropy with for different near the N-C transition. For low densities exhibit a single peak at , but with increasing a secondary maximum gets developed at . For the maximum at and becomes of equal height. Beyond the global maximum of jumps to , making the N-C transition to be first order.

Unlike the case, there is no way to obtain an analytic expression for . For m=3, the numerically determined for different is shown in Fig. (20). From the data, we obtain

 ρN−Cc=0.62875+0.107k,m=3. (23)

We note that this expression has the same form as for [see Eq. (20)].

For , we proceed as follows. The transition density is bounded from above by the spinodal density , the density at which the entropy at changes from a local maximum to local minimum. is obtained from the condition and we obtain

 ρN−Cs = −m+2km2−m√1−4k+4k2m2(1−m−km+km2)≥ρN−Cc, (24) = √m1+√m+B1(m)k+O(k−2), k→∞, (25) = 1−1√m+B2(k)m+O(m−3/2), m→∞, (26)

where and . The spinodal density is compared with in Fig. 20. From Eq. (26), it follows that and tends to one when . The limit corresponds to the continuum limit. In this limit . Thus, it is not clear whether the nematic-columnar phase transition will exist in the continuum.

### vi.3 The columnar–sublattice phase boundary

The dependence of the C-S phase boundary on and may be determined by estimating the entropy for the C and S phases close to full packing. We approximate the entropy of the C phase by the entropy of the fully aligned C phase. Since the heads of the rectangles are all in either even or odd rows/columns, by ignoring the unoccupied rows/columns, the calculation of entropy reduces to an one dimensional problem of rods. The mean number of holes in a row is , and the mean number of rods (of length ) in a row is . There are such rows. The number of ways of arranging the rods and holes on a row is

 Ωrow=[L(1−ρ)+ρLmk]![L(1−ρ)]![ρLmk]!, (27)

such that the total number of ways of arranging the rectangles is . Hence, , the entropy per site of the columnar phase is given by , which for densities close to is

 SC≈1−ρmln[ekm(1−ρ)]+O[(1−ρ)2]. (28)

We now estimate the entropy for the sublattice phase. At full packing, the head of each rectangle is on one of sublattices. Ignoring the sites belonging to other sublattices, it is easy to see that each configuration of rectangles can be mapped on to a configuration of rods of length on a lattice of size . When , by solving the full packed problem of rods on strips, it is known that the entropy per unit site is  Ghosh and Dhar (2007). For rectangles of size we obtain

 SS(ρ=1)≈lnkm2k2, k≫1, (29)

where is the entropy per site of the S phase. For densities close to 1 () we estimate the correction term to by removing fraction of rectangles at random from the fully packed state. Here, we ignore the entropy of the holes, assuming that the holes form bound states. This gives the entropy of the sublattice phase, close to the full packing to be approximately

 SS≈lnkm2k2−1m2k[(1−ρ)ln(1−ρ)+ρlnρ],  k≫1, (30)

Comparing Eqs. (28) and (30) up to the leading order we obtain the critical density for the C-S transition to be

 ρC−Sc≈1−A2mk2, k≫1, (31)

where is a constant.

Given the above asymptotic behavior of the phase boundaries, we expect that the phase diagram for will be similar to that obtained for for large . For small , we check that for rectangles of size , there are two transitions just as in rectangles. Thus, we are led to conjecture that the phase diagram for will be qualitatively similar to that for for all