The columnar-disorder phase boundary in a mixture of hard squares and dimers

TL;DR

This paper investigates the phase transition in a lattice model of hard squares and dimers, estimating the boundary between disordered and ordered phases using analytical methods that align well with simulations.

Contribution

It introduces an analytical approach to estimate the phase boundary in a mixture of hard squares and dimers, validated by Monte Carlo simulations.

Findings

Phase boundary estimates agree with simulations

Identified a Kosterlitz-Thouless transition

Characterized the nature of ordered and disordered phases

Abstract

A mixture of hard squares, dimers and vacancies on a square lattice is known to undergo a transition from a low-density disordered phase to high-density columnar ordered phase. Along the fully packed square-dimer line, the system undergoes an Kosterliz-Thouless type transition to a phase with power law correlations. We estimate the phase boundary separating the ordered and disordered phases by calculating the interfacial tension between two differently ordered phases within two different approximation schemes. The analytically obtained phase boundary is in good agreement with Monte Carlo simulations.

Equations

$$e^{-\sigma L_y}=\frac{{\sum }_I{Z}^{(I)}}{Z^{(0)}},$$

$$Z^{(I)}=Z_L^{(I)}{Z}_R^{(I)},$$

$$Z_L^{(I)}=$$

$$Z_R^{(I)}=$$

$$Z^{(0)}=\prod _{i=1}^{L_y/2}\Omega (L_x,0).$$

$$G(y,\Delta )=\sum _{\ell =0}^{\infty }\Omega (\ell ,\Delta )y^{\ell +\Delta /2},$$

$$G(y,0)$$

$$G(y,1)$$

$$G(y,0)$$

$$G(y,1)$$

$$f(y)=$$

$$-(z_h+z_v+z_0^2)y+1,$$

$$\Omega (\ell ,\Delta )=a(\Delta ){\lambda }^{\ell }[1+O(e^{-c\ell })],\ell \gg 1,$$

$$\lambda =\frac{1}{y_1}.$$

$$a(0)=$$

$$a(1)=$$

$$\Omega (\ell ,\Delta )=z_0\Omega (\ell ,\Delta -1)+z_h\Omega (\ell ,\Delta -2),\Delta \ge 2.$$

$$H(\ell ,x)=\sum _{\Delta =0}^{\infty }\Omega (\ell ,\Delta )x^{\Delta }.$$

$$H(\ell ,x)=\frac{\Omega (\ell ,1)x+\Omega (\ell ,0)(1-z_{0}x)}{1-z_{0}x-z_h{x}^2}.$$

$$x_{\pm }=\frac{-(z_0\pm \sqrt{z_0^2+4z_h})}{2z_h}.$$

$$\Omega (\ell ,\Delta )=\sum _{i=\pm }\frac{a(1)+a(0)z_h{x}_i}{(z_0+2z_h{x}_i)x_i^{\Delta }}{\lambda }^{\ell },\ell \gg 1.$$

$$a(\Delta )=\frac{p_{+}}{x_{+}^{\Delta }}+\frac{p_{-}}{x_{-}^{\Delta }},\Delta \ge 0,$$

$$p_{\pm }=\frac{a(1)+a(0)z_h{x}_{\pm }}{z_0+2z_h{x}_{\pm }}.$$

$$Z_L^{(I)}$$

$$Z_R^{(I)}$$

$$\max ({\eta }_i,{\eta }_{i+1})=\frac{|{\eta }_i-{\eta }_{i+1}|+{\eta }_i+{\eta }_{i+1}}{2},$$

$$Z^{(I)}$$

$$\sum _I{Z}^{(I)}$$

$$\displaystyle\bigg{(}p_{+}\frac{x_{+}\sqrt{\lambda}+1}{x_{+}\sqrt{\lambda}-1}+p_{-}\frac{x_{-}\sqrt{\lambda}+1}{x_{-}\sqrt{\lambda}-1}\bigg{)}\Bigg{]}^{L_{y}/2}.$$

$$Z^{(0)}=\bigg{[}\Omega(L_{x},0)\bigg{]}^{L_{y}/2}=\bigg{[}a(0)\lambda^{L_{x}}\bigg{]}^{L_{y}/2}.$$

$$e^{-\sigma L_{y}}=\Bigg{[}(\frac{z_{v}}{\lambda}+\frac{z_{s}}{\lambda^{2}})\bigg{(}p_{+}\frac{x_{+}\sqrt{\lambda}+1}{x_{+}\sqrt{\lambda}-1}+p_{-}\frac{x_{-}\sqrt{\lambda}+1}{x_{-}\sqrt{\lambda}-1}\bigg{)}\Bigg{]}^{\frac{L_{y}}{2}}.$$

$$\left[\frac{z_v}{\lambda }+\frac{z_s}{{\lambda }^2}\right]\left[p_{+}\frac{x_{+}\sqrt{\lambda }+1}{x_{+}\sqrt{\lambda }-1}+p_{-}\frac{x_{-}\sqrt{\lambda }+1}{x_{-}\sqrt{\lambda }-1}\right]=1.$$

$$z_s^{1/4}+z_d^{1/2}+z_0=1,$$

$${\rho }_s$$

Full text

I Introduction

It is well-known that the universal features of continuous phase transitions may be determined by studying simple model systems. Among such models, hard core lattice gases (HCLGs) are one of the simplest. In these lattice models, particles interact only through excluded volume interactions and phase transitions, if any, depend only on the shape of the constituent particles and do not depend on temperature. Hence, they have often been termed as geometrical or entropy-driven transitions. Despite is simplicity, HCLGs exhibit different broken-symmetry phases like solid-like sublattice order, columnar or smectic phase with partial translational order, and nematic phase with orientational order, depending on the particle shape. Many monodispersed systems of differently shaped particles have been studied in the literature. Examples include triangles Verberkmoes and Nienhuis (1999a), squares Bellemans and Nigam (1967, 1966); Ree and Chesnut (1966); Ramola and Dhar (2012); Nath et al. (2016); Mandal et al. (2017), Y-shaped molecules Ruth et al. (2015), dimers Kasteleyn (1961); Temperley and Fisher (1961); Huse et al. (2003); Nicholls et al. (2017), trimers Ghosh et al. (2007), tetrominoes Mao et al. (2002); Barnes et al. (2009), rods Ghosh and Dhar (2007); Kundu et al. (2013); Oettel et al. (2016); Vigneshwar et al. (2017); Gschwind et al. (2017), rectangles Kundu and Rajesh (2014, 2015); Nath et al. (2015); Gurin et al. (2017), discs Fernandes et al. (2007); Nath and Rajesh (2014), and hexagons Baxter (1980), the latter being the only exactly solvable case. Polydispersity can hardly be avoided in experiments. However, compared to the monodispersed systems, the phenomenology of HCLG models of mixtures of particles of different shapes is less understood.

Amongst mixtures, the best studied example is that of depletion interaction in mixtures of particles with small excluded volume and particles with larger excluded volume. When the excluded volumes are on-site and first nearest neighbour, then in two dimensions, it is known from different numerical studies that there is a critical line ending in a tricritical point separating a low-density disordered fluid-like phase from a high-density solid-like sublattice phase Oliveira and Stilck (2015); Poland (1984); Liu and Evans (2001); van Duijneveldt and Lekkerkerker (1993); Verberkmoes and Nienhuis (1999b); Frenkel and Louis (1992). The nature of the transition is similar to that of the transition observed in the system with only larger particles. Similar demixing transition occurs in binary mixture of large and small cubes Dijkstra and Frenkel (1994); Dijkstra et al. (1994); Lafuente and Cuesta (2003). Other examples of mixtures that have been studied on lattices include bidispersed rods Stilck and Rajesh (2015); Kundu et al. (2016) in which the phase diagram is richer than the monodispersed case, but without the appearance of any new phase. In a recent paper Ramola et al. (2015), a mixture of hard squares and dimers was studied and, very interestingly, it was demonstrated that the phase diagram consists of two critical lines with continuously varying exponents that meet at a point called as the Askin-Teller-Kosterlitz-Thouless point. It provided the first example of a HCLG whose critical properties vary with the composition. In this paper, we focus on aspects of this mixture of dimers and squares.

We first discuss the known results for monodispersed systems of dimers or squares. The dimer model Kasteleyn (1961); Temperley and Fisher (1961); Fisher (1961); Fisher and Stephenson (1963); Heilmann and Lieb (1972); Huse et al. (2003); Heilmann and Lieb (1970), in which each particle occupies $1\times 2$ or $2\times 1$ sites on the square lattice is the simplest model for anisotropic particles. It is exactly solvable in the fully packed limit. In this limit, the system is critical and the correlations between dimers decay with separation $r$ as as $r^{-2}$ Fisher and Stephenson (1963). At densities different from the fully packed limit, the system is disordered Heilmann and Lieb (1972, 1970).

The hard square model Bellemans and Nigam (1967, 1966); Ree and Chesnut (1966); Ramola and Dhar (2012); Mandal et al. (2017) in which each particle occupies $2\times 2$ sites on the square lattice is the prototypical model to show columnar order. As density is increased, the system undergoes a continuous phase transition from a low-density disordered phase to a high density columnar ordered phase with a four-fold symmetry. In the columnar phase, the squares preferentially occupy either even or odd rows or even or odd columns, thus breaking translational order in only one of the two directions. The continuous transition belongs to the Ashkin-Teller universality class with the correlation length exponent $\nu\approx 0.92$ Zhitomirsky and Tsunetsugu (2007); Feng et al. (2011); Ramola et al. (2015).

The mixture of hard squares and dimers was studied both numerically and analytically in Ref. Ramola et al. (2015). The system undergoes a transition from a square-rich columnar phase to a dimer-rich disordered phase across a critical line along which the critical exponents continuously vary depending on the composition of the mixture, consistent with the Ashkin-Teller universality class. On the fully packed line, it was shown that the system undergoes a Korterlitz-Thouless type transition from a columnar phase to a power law correlated phase as the dimer density is increased. Along the fully packed line, the configurations of dimers and squares may be mapped onto a height field. By writing an effective Hamiltonian for the two dimensional height field, it was possible to theoretically explain the numerically obtained results, along the fully packed line Ramola et al. (2015). However, the height mapping does not allow the phase boundary to be determined, as the relation between the rigidity and microscopic parameters is not easy to establish. In this paper, we determine the phase diagram within an approximation scheme, and compare with the numerically obtained phase boundary.

Estimates of phase boundaries in systems showing columnar order, obtained from standard approximation schemes like density functional theory, high density expansions, Flory-type approximations, etc., are quite poor [see Ref. Mandal et al. (2017) for a tabulation of results for the hard square model]. In recent work Nath et al. (2016); Mandal et al. (2017), we described a systematic way of determining the interfacial tension between two differently ordered columnar phases in terms of number of defects and overhangs in the interface. The estimates obtained from the interfacial tension are in good agreement with the numerical results for the hard square gas Nath et al. (2016); Mandal et al. (2017) as well as the hard rectangle gas Nath et al. (2016). In this paper, we use the same method to obtain the phase boundary for the mixture of squares and dimers. We estimate the interfacial tension between two different columnar ordered phases and by setting it to zero, we obtain limiting condition for the stability of columnar phase. First, we assume that the interface between the two ordered phases has no overhangs and that the ordered phases have perfect order. We improve the estimate for the phase boundary by allowing the interface to have overhangs of height one. The results are summarized in Fig. 6. For instance, along the fully packed line our estimates for the critical density for squares are within $8\%$ of the numerical result.

The remainder of the paper is organized as follows. In section II, we define the model precisely and give an outline of the calculation. The details of the calculations are presented in Sec. III and Sec. IV, and the results are compared with results from Monte Carlo simulations. We end with a discussion in Sec. V.

II Model and outline of calculation

Consider a square lattice of size $L_{x}\times L_{y}$. The lattice may be occupied by particles of three different shapes: squares, horizontal dimers and vertical dimers of size $(2\times 2)$, $(2\times 1)$ and $(1\times 2)$ respectively. The particles interact only through excluded volume interaction, i.e. no two particles may overlap. We associate with each square, horizontal dimer, vertical dimer and vacancy ($1\times 1$) activities $z_{s}$, $z_{h}$, $z_{v}$ and $z_{0}$ respectively. We will refer to the bottom left corner of a particle as its head.

Depending on the values of the activities, the system may exist in a disordered fluid like phase or in an ordered phase which has columnar order Ramola et al. (2015). In the columnar phase, the heads of squares and vertical (horizontal) dimers preferably occupy either even or odd rows (columns) with equal fraction on an average in even or odd columns (rows).

The aim of the paper is to determine the phase boundary between the columnar and disordered phases. This is done by estimating the interfacial tension $\sigma(z_{s},z_{h},z_{v},z_{0})$ between two ordered columnar phases, and by equating it to zero, we obtain estimates for the critical activities and densities. Let the phase in which majority of heads of squares and vertical dimers are in even (odd) rows be called even (odd) phase. To compute $\sigma(z_{s},z_{h},z_{v},z_{0})$, we impose an interface in the system by fixing squares at the left boundary to be even and those at the right boundary to be odd. A snapshot of a typical equilibrium configuration seen in a Monte Carlo simulation of a system with these boundary conditions is shown in Fig. 1. There is a left phase with even squares and vertical dimers separated from a right phase with odd squares and vertical dimers by an interface. The horizontal dimers could be even or odd in both phases. To define an unique position of the interface for any allowed configuration of particles, we adopt the convention that the boundary between the left or even phase and right or odd phase is placed as far left as possible. With this convention, there is a well-defined interface. The bulk phases have only few defects. By defects, we mean particles of the wrong type (odd in even phase or even in odd phase). When the defects are removed, a fully ordered columnar phase is recovered.

Let $Z^{(0)}$ and $Z^{(I)}$ be the partition function of the system without and with an interface $I$ respectively. The interfacial tension $\sigma$ is defined as

[TABLE]

Due to the nature of interactions between particles being hard core, we can write the partition function of the system in the presence of an interface as a product of the partition functions of the left and right phases, i.e.

[TABLE]

where $L$ and $R$ denote left and right.

$Z^{(I)}$ cannot be calculated for an arbitrary interface $I$. We therefore calculate it within two approximations. As a first approximation, we consider the simplest case where we ignore overhangs in the interface and defects in the bulk. In this simplified model, the interface is defined by the position of right boundary of the left phase, and denoted by $\eta_{i}$ (see Fig. 2). Since we assume perfect columnar order for the left and right phases, the partition functions for both left and right phases are a product of partition functions of tracks made up of two adjacent rows ($L_{y}/2$ of them):

[TABLE]

where $\Omega(\ell,\Delta)$ is the partition function of a track of two rows with a shape as shown in Fig. 3. In the region corresponding to $\Delta$, only horizontal dimers can be placed. The right hand side of Eq. (3) follows from the fact that, for the left phase, there must be either a vertical dimer or square touching the interface. The partition function of the system without an interface is

[TABLE]

In the second approximation, we allow the interface to have overhangs of height one. We still do not allow defects in the bulk. The calculation is on the same lines as that described above. This allows us to obtain an improved estimate of the critical parameters.

The calculation of the interfacial tension involves determining the partition function $\Omega(\ell,\Delta)$ of a track of two rows, and is done in the next section.

III Calculation of two-row partition function $\Omega(\ell,\Delta)$

In this section we calculate $\Omega(\ell,\Delta)$, the partition function of a track of two rows with shape as shown in Fig. 3. Consider the generating function defined as

[TABLE]

where the power of $\sqrt{y}$ is the number of sites in the system.

First, consider the case $\Delta=0,1$. $G(y,0)$ and $G(y,1)$ obey simple recursion relations which are shown diagrammatically in Fig. 4 and may be written as

[TABLE]

These are easily solved to yield

[TABLE]

where the function $f(y)$ in the denominator is

[TABLE]

a third order polynomial in $y$. Let $y_{1}$ be smallest root of $f(y)=0$. Then, it is clear that

[TABLE]

where

[TABLE]

The prefactor $a(\Delta)$ for $\Delta=0,1$ is determined by calculating the coefficient of $y^{\ell+\Delta/2}$. It is easily checked that

[TABLE]

Equations (12)–(15) determine $\Omega(\ell,\Delta)$ for $\Delta=0,1$.

We now calculate $\Omega(\ell,\Delta)$ for $\Delta\geq 2$. The recursion relation obeyed by $\Omega(\ell,\Delta)$ is shown diagrammatically in Fig. 5, and may be written as

[TABLE]

To solve Eq. (16), define the generating function

[TABLE]

Multiplying Eq. (16) by $x^{\Delta}$ and summing over $\Delta$, we obtain a linear equation for $H(\ell,x)$ that may be solved to yield

[TABLE]

The generating function $H(\ell,x)$ has two simple poles determined by the zeros of the quadratic equation $1-z_{0}x-z_{h}x^{2}=0$. These are

[TABLE]

Expanding the denominator of Eq. (18) about $x_{+}$ and $x_{-}$ we obtain the coefficient of $x^{\Delta}$ to be

[TABLE]

Using Eqs. (12) and (20), we can write

[TABLE]

where

[TABLE]

IV Interfacial tension and critical parameters

We now calculate interfacial tension $\sigma(z_{s},z_{h},z_{v},z_{0})$. The results for the interface without overhangs is in Sec. IV.1 and for the interface with overhangs in Sec. IV.2.

IV.1 Without overhangs

For large $L_{x}$ and $\eta_{i}$, using Eqs. (3), (4) and (12), the partition functions of the left and right domains may be written as

[TABLE]

Taking the product of Eqs. (22) and (23) and using the relation

[TABLE]

we obtain the partition function of the system with interface

[TABLE]

Note that $Z^{(I)}$ depends only on the difference $|\eta_{i}-\eta_{i+1}|$ and not on $\eta_{i}$. Summing over all configurations $I$ is equivalent to summing over all differences $(\eta_{i}-\eta_{i+1})$. Performing the sum and using Eq. (21), we obtain

[TABLE]

The partition function of the system without an interface can be easily calculated using Eq. (5)

[TABLE]

The interfacial tension, as defined in Eq. (1) may be determined from Eqs. (26) and (27) to be

[TABLE]

The phase boundary corresponds to the values of the parameters at which the interfacial tension vanishes. This immediately gives

[TABLE]

The phase boundary obtained from Eq. (29) is shown by the magenta line in Fig. 6 and Fig. 6. Here, the activities have been normalized using

[TABLE]

where $z_{d}=z_{h}=z_{v}$, such that a two dimensional phase diagram may be obtained. It shows transitions at $z_{s}^{1/4}=0.725$ along square-vacancy (SV) line and at $z_{s}^{1/4}=0.616$ along square-dimer (SD) line. These should be compared with results from Monte Carlo simulations Ramola et al. (2015): $z_{s}^{1/4}=0.759$ along SV line and at $z_{s}^{1/4}=0.692$ along SD line.

The density of sites occupied by squares $\rho_{s}$, horizontal dimers $\rho_{h}$ and vertical dimers $\rho_{v}$ may be calculated from the partition function $Z^{(0)}$ as:

[TABLE]

where the factor 4 and 2 accounts for the area of a square and a dimer respectively. Substituting for $Z^{(0)}$ from Eq. (27), the densities may be written in the thermodynamic limit $L_{x}\rightarrow\infty$, $L_{y}\rightarrow\infty$ as:

[TABLE]

where total dimer density $\rho_{d}=\rho_{h}+\rho_{v}$ and density of vacancy $\rho_{0}=1-\rho_{d}-\rho_{s}$.

In the density plane, it shows transitions at $\rho_{s}=0.928$ along SV line and at $\rho_{s}=0.714$ along SD line. The Monte Carlo results show transitions at $\rho_{s}=0.932$ along SV line and at $\rho_{s}=0.843$ along SD line. Our estimation of critical activities along SV and SD lines agree satisfactorily with the numerical results (see Sec. IV.3).

IV.2 With overhangs

In the calculation presented in Sec. IV.1, the interface was modelled as having no overhangs. We now allow the interface to have a certain subclass of overhangs, and obtain an improved estimate for the phase boundary. To be able to do the calculation in the presence of overhangs, we first reformulate the calculation in the absence of overhangs in terms of a weighted, directed walk.

The interface with no overhang (see Fig. 2) may be visualized as a partially directed self avoiding walk (PDSAW) from top to bottom, where the walk is not allowed to take a step in the upward direction, but allowed to take steps in the leftward, rightward, and downward directions as long as the walk is self avoiding. We denote a downward step by $\mathbb{D}$, a leftward step of length $\Delta$ by $\mathbb{L}_{\Delta}$ and rightward step of length $\Delta$ by $\mathbb{R}_{\Delta}$. To maintain self avoidance, every set of consecutive leftward or rightward steps must be preceded (or followed) by a downward step. Thus, the different PDSAWs may be enumerated by arbitrary concatenation of substrings $\mathbb{D}$, $\mathbb{D}\mathbb{L}_{\Delta}$, and $\mathbb{D}\mathbb{R}_{\Delta}$ where $\Delta=1,2,\ldots$, and the length of the walk in the vertical direction is given by the number of $\mathbb{D}$s in the string. The formal generating function of these strings may be written as

[TABLE]

where

[TABLE]

The generating function $\mathbb{G}_{0}$ generates all possible walks of all possible lengths along the vertical direction. However, it does not assign weights to a walk.

To assign weights to each walk, we have to determine the weights $D$, $L_{\Delta}$ and $R_{\Delta}$ that correspond to the substrings $\mathbb{D}$, $\mathbb{L}_{\Delta}$, and $\mathbb{R}_{\Delta}$. To do so, we determine the weight of a PDSAW by taking the ratio of Eq. (25) and Eq. (27) to obtain

[TABLE]

where $R$ or $L$ in Eq. (40) is chosen depending on whether the step is in the rightward or leftward direction. From Eqs. (39) and (40), we immediately read out

[TABLE]

Consider now the weighted generating function

[TABLE]

where

[TABLE]

It is straightforward to see that, in terms of the interfacial tension $\sigma$, $\mathcal{G}_{0}$ may be written as

[TABLE]

Equation (46) is convergent for all $\sigma>0$, and diverges at $\sigma=0$, corresponding to the transition point. From Eq. (44), the divergence of $\mathcal{G}_{0}$ corresponds to the condition

[TABLE]

Substituting for $D$, $\widetilde{R}$ and $\widetilde{L}$ from Eqs. (41) and (45), we obtain the same equation for the phase boundary as obtained earlier in Eq. (29).

We now allow the interface to have overhangs of height one. If an overhang is present, then a horizontal line on the dual lattice will intersect the interface more than once. A schematic diagram of an interface with overhangs is shown in Fig. 7. The interface with overhangs can no longer be interpreted as a PDSAW, as the walk is now allowed to take upward steps. However, restricting the overhangs to height one implies that each upward step has to followed by a downward step before another upward step can be taken. We now generalize the PDSAW to take into account these upward steps.

An overhang will be termed as right or left overhang depending on whether it appears in the right (as shown by yellow with index A in Fig. 7) or left phase (as shown by yellow with index B in Fig. 7). We now separately determine the generating function of all possible walks with right and left overhangs of height one.

IV.2.1 Right Overhangs

There are two kinds of right overhangs depending on whether the first downward step is followed by rightward steps [see Fig. 8(a)] or leftward steps [see Fig. 8(b)]. We split each of these into two parts: the initial part shown by red and the remaining part shown by yellow in Fig. 8. This amounts to restricting to configurations where a horizontal dimer does not cross the right boundary of shaded region in Fig. 8. We denote these parts of the walk by $\mathbb{W}_{R}^{(1)}(n)$, $\mathbb{W}_{R}^{(2)}(n)$ and $\mathbb{O}_{R}(n_{1},n_{2})$ respectively, where the superscript refer to the two types of right overhangs and subscript $R$ stands for right. Thus a generic right overhang is represented by

[TABLE]

We now determine the weights for these overhangs. It is clear that

[TABLE]

where the weight of $\mathbb{W}_{R}^{(i)}(n)$ may be determined by considering an interface with only one overhang. Here $\omega(\Delta)$ is the partition function of a track of width one and length $\Delta$, and appears in the weights because the shaded region in Fig. 8 may be occupied by horizontal dimers.

The partition function $\omega(\Delta)$ is easily determined. Define the generating function

[TABLE]

where the power of $x$ represents the number of sites in the system. The recursion relation obeyed by $\widetilde{\omega}(x)$ is shown diagrammatically in Fig. 9 and may be written as

[TABLE]

with solution

[TABLE]

Expanding the denominator about its two roots $x_{\pm}$ [see Eq. (19)] and using Eq. (51), we can write

[TABLE]

where

[TABLE]

Using Eq. (54), the weights $W_{R}^{(1)}$ and $W_{R}^{(2)}$ may be rewritten as

[TABLE]

[TABLE]

Now consider the weight $O_{R}$ defined as

[TABLE]

where $O_{R}^{(i)}$ depends on the different ways the particles at the edge of the overhangs may be placed (see Fig. 10). The weight of these four kinds of overhangs may be determined in a straightforward manner by considering example of an interface with only downward step and one overhang :

[TABLE]

The sums over $n_{i}$ may be expressed in term of the generating function $G(y,0)$ [see Eq. 6]. This gives

[TABLE]

The total weight of all walks with right overhang may now be computed. Let $\mathbb{U}_{R}$ represent all possible walks with right overhang and having total weight $U_{R}$. Then

[TABLE]

where the second term in right hand side represent possible multiple overhang and the third term, the possibility of right steps after the overhang. $U_{R}$ may be rewritten as

[TABLE]

where $\widetilde{R}$ is defined in Eq. (45).

IV.2.2 Left Overhangs

The weight for left overhangs may be calculated in a manner similar to that for right overhangs. There are two kinds of left overhangs depending on whether the first downward step is followed by leftward steps [see Fig. 11 (a)] or rightward steps [see Fig. 11 (b)]. We split each of these into three parts, where the initial red line represents the first two downward steps and the intervening region, the middle yellow line represents overhangs and final red line represents leftward steps. We denote these parts symbolically by $\mathbb{W}_{L}^{(1)}(n)$, $\mathbb{W}_{L}^{(2)}(n)$, $\mathbb{O}_{L}(n_{1},n_{2})$ and ${\mathbb{L}}^{\prime}(n^{\prime})$ respectively, where the superscript (1) and (2) refer to the two kinds of left overhangs, and subscript $L$ stands for left. Thus, a generic left overhang of first kind [see Fig. 11 (a)] is represented by

[TABLE]

and that of the second kind [see Fig. 11 (b)] is represented by

[TABLE]

We now determine the weights of the different parts constituting the left overhangs. For the overhang of first kind, the total weight associated with the initial red part $\mathbb{W}_{L}^{(1)}(n)$ may be written as

[TABLE]

where $L_{\Delta}$ and $\widetilde{L}$ are given respectively in Eqs. (42) and (45). The weight $O_{L}$ associated with the left overhang is identical to that of the right overhang of similar shape, i.e.,

[TABLE]

where $O_{R}$ is given in Eq. (63).

Now, consider the left overhangs of the second kind [see Fig. 11 (b)]. The total weight associated with the first two downward steps and the first overhang may be written as

[TABLE]

where the summation over $n_{2}$ contributes $\sqrt{O_{L}}$. It is convenient to change the variable to $m=n+n_{1}$. In this new variable Eq. (70) may be rewritten as

[TABLE]

We define the functions

[TABLE]

where $i=1,2$ and $G(y,0)$ is the generating function as determined in Eq. (7). Doing the sums in Eq. (71), we obtain

[TABLE]

where

[TABLE]

with $b_{\pm}$ as defined in Eq. (55).

Finally, we calculate the total weight associated with the final set of leftward steps that may be taken after the overhangs:

[TABLE]

The total weight of all walks with left overhang may now be computed. Let $\mathbb{U}_{L}$ represent all possible walks with at least one left overhang. Let the total weight associated with these walks be $U_{L}$. Then we obtain

[TABLE]

where the second term in the right hand side represents multiple overhangs and the third term accounts for the leftward steps after the final overhang. Rewriting,

[TABLE]

IV.2.3 Phase boundary

The formal generating function of all possible walks including those with overhangs may be written as

[TABLE]

For every term in the expansion, one may uniquely identify a walk from top to bottom. The weighted generating function corresponding to $\mathbb{G}_{ov}$ may be written as

[TABLE]

As discussed earlier [see discussion following Eq. (46)], the generating function $\mathcal{G}_{ov}$ diverges at the transition point when the interfacial tension $\sigma$ vanishes, and this condition is equivalent to

[TABLE]

where $D$, $\widetilde{R}$, $\widetilde{L}$, $U_{R}$ and $U_{R}$ are given in Eqs. (41), (45), (65) and (77) respectively.

In Fig. 6, blue line represents the critical line for activity with overhang. It shows transitions at $z_{s}^{1/4}=0.733$ along SV line and at $z_{s}^{1/4}=0.642$ along SD line. The density plot of critical line with overhang is shown in the Fig. 6 by blue line. It shows transitions at $\rho_{s}=0.934$ along SV line and at $\rho_{s}=0.779$ along SD line.

IV.3 Monte Carlo simulations

In this section we determine the phase boundary numerically using grand canonical Monte Carlo simulations. Conventional algorithms that include only local evaporation and deposition moves often do not equilibrate the system, within available computer time, at high densities due to long-lived metastable states. Here, we implement an algorithm that updates two rows at a time using a transfer matrix based Monte Carlo algorithm Ramola et al. (2015). The algorithm not only succeeds in equilibrating the system at high densities, but also at full packing. This algorithm is a generalization of a cluster algorithm that is able to equilibrate systems of particles with large excluded volume at high densities Kundu et al. (2013, 2012). In this Monte Carlo scheme Ramola et al. (2015), all particles that are fully contained in a $2\times L$ track, consisting of two adjacent rows or columns of length $L$, are evaporated. The track is refilled with particles according to the correct weights in the partition function for the track, subject to the constraints induced by particles protruding into the track. The calculation of the restricted partition function is done using standard transfer matrix technique, details of which may be found in the Supplemental Information of Ref. Ramola et al. (2015).

The order parameter $Q$ is defined as

[TABLE]

where $\rho_{er}$, $\rho_{ec}$, $\rho_{or}$ and $\rho_{oc}$ are the densities of heads of particles (including both dimers and squares) in even rows, odd rows, even columns and odd columns respectively. Consider the second moment of the order parameter $Q$, defined as

[TABLE]

where $L\times L$ is the system size. Near the transition point $\chi$ obeys following scaling law

[TABLE]

where $\epsilon$ is the deviation from critical point $\epsilon=z_{s}-z_{c}$, and $f$ is the finite size scaling function, and $\gamma$ and $\nu$ are the usual critical exponents. Since the model belongs to Ashkin-Teller universality class, one knows that the exponent $\gamma/\nu=7/4$ Ramola et al. (2015), independent of the parameters. From Eq. (83), we see that $\chi/L^{7/4}$ for different $L$ should cross at $\epsilon=0$, allowing the critical point to be estimated. An example is shown in Fig. 12.

We fix the activity of dimer $z_{d}$ and vary the activity of square $z_{s}$ to get its critical value $z_{c}$. Critical value of $z_{0}$ may be obtained from the normalization condition given in Eq. (30). The data of full phase boundary for activity plane and density plane are shown by blue dots in Fig. 6 and Fig. 6 respectively.

V Discussion

To summarize, we calculated the interfacial tension between two differently ordered phases in a mixture of hard squares and dimers, within two approximation schemes. The estimates for the phase boundary between the ordered columnar phase and disordered fluid-like or power law correlated phase was obtained by setting the interfacial tension to zero. In the first calculation, we modeled the interface as having no overhangs. The estimates were improved by extending the calculations to an interface where overhangs of height one were allowed. The estimates that we obtain for critical parameters are shown in Fig. 6 and 6, and are in good agreement with results from Monte Carlo simulations. The deviation from the numerical results are largest along the fully packed square-dimer (SD) line. Along this line, we obtain the critical parameters to be $z_{s}^{1/4}=0.616$, $\rho_{s}=0.714$ for interfaces without overhang, and $z_{s}^{1/4}=0.642$, $\rho_{s}=0.779$ for interfaces with overhangs. These are to be compared with results from Monte Carlo simulations: $z_{s}^{1/4}=0.692$ and $\rho_{s}=0.843$. Along the square vacancy (SV) line, the calculation reproduces the results in Ref. Nath et al. (2016) for interfaces without overhangs, but corrects the result for interfaces with overhangs. For the latter case, it was erroneously assumed in Ref. Nath et al. (2016) that the contributions from left and right interfaces are identical.

In our calculations we assumed that the ordered phases have perfect order, thereby ignoring the presence of defects in the bulk phases. Defects may be included in a systematic manner, as was done for the case of $2\times d$ rectangles Nath et al. (2016). However, it was found that the corrections appearing from including overhangs were more dominant than that arising from including defects when $d$ was small as is the case for hard squares. At the same time, the calculations for including defects is involved and including the effect of two defects appears a formidable task. We have, therefore, ignored the corrections due to the presence of defects.

The calculations could also be improved by taking into account multiple defects using the pairwise approximation Mandal et al. (2017). Such an approach could improve the estimates for the phase boundary along the fully packed square-dimer line. This is a promising area for future study.

It will also be interesting to study mixtures of $m\times m$ squares and $m\times 1$ rods, $m=2$ being the square-dimer problem. Along the fully packed line, it would be possible to map the configurations to a vector height field. The ordered phase has a $2m$ fold symmetry. This will possibly change the nature of the transition from the ordered to disordered phases. Also, along the fully packed line, it raises the possibility of an intermediate hexatic phase. For $m>2$, the correlations even the fully packed $m$-mer problem is not known Ghosh et al. (2007). Thus, $m=3$ system (trimers+squares) would be a good starting point.

Acknowledgements.

We thank Trisha Nath for helpful discussions.