## Abstract

We report an analysis of the change in character of the fluid-to-solid transition in a quasi-two-dimensional hard-sphere colloid system as the confining wall separation changes from one to two hard-sphere diameters. Our analysis is based on a study of the bifurcation of solutions of the integral equation for the singlet density, using both direct and pair correlation function representations. Two approximations used in previous bifurcation analyses of freezing are improved in the work reported in this paper. The results of the analysis correctly predict the qualitative change in freezing as a function of the confining wall separation and density, specifically the change from a fluid-to-one-layer hexagonal solid transition to a fluid-to-two-layer square solid transition at a wall separation of 1.6 hard-sphere diameters. The numerical predictions of the theory are in semi-quantitative agreement with simulation data for the density dependence of the liquid–solid transition line for wall separations up to 1.6 hard-sphere diameters.

## 1. Introduction

Among the interesting features that differentiate confined matter from unconfined matter is the occurrence of packing structures in the former which do not exist in the latter. For example, densely packed colloid suspensions confined between smooth parallel plates only a few particle diameters apart exhibit a number of structural phase transitions as the plate separation and/or the density is varied. Specifically, as functions of the density and the distance between the two walls, a quasi-two-dimensional system undergoes the transitions (van Winkle & Murray 1986; Murray *et al*. 1990; Murray 1992)(1.1)In (1.1), the sequence of transitions displayed on the right-hand side is ordered with respect to increasing density and increasing wall separation; the triangle symbol denotes a layer with hexagonal lattice symmetry and the square symbol denotes a layer with square lattice symmetry, and the integers in front of these symbols refer to the number of layers of particles between the confining plates. The first analysis of these phase transitions of which we are aware is based on free-volume considerations. Loosely put, the successive transitions are driven by the competition between the free volume available to *n* square-packed layers, with smaller overall thickness than *n* hexagonal-packed layers, and the free volume associated with more efficient packing within a hexagonal-packed layer than in a square-packed layer.

This paper is concerned with a theoretical analysis of the dependence on wall separation and density of the fluid→1▵ and the transitions in a quasi-two-dimensional colloid system. Our primary concern in this paper is the development of an analytical representation that correctly captures the overall characteristics of the phase transition as a function of density and wall separation. For this purpose, we examine the wall separation dependence of the stability of the fluid phase with respect to hexagonal- or square-packed lattices, as monitored by the bifurcation of solutions of the first equation of the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy, and the equivalent Lovett–Mou–Buff–Wertheim (LMBW) equation (Lovett *et al*. 1976; Wertheim 1976). Our expectation is that this approach will provide a qualitatively valid interpretation of the transitions, albeit with only semi-quantitative accuracy.

Kirkwood & Monroe first suggested the identification of the bifurcation point of solutions to the first equation of the BBGKY hierarchy with the freezing transition (Kirkwood & Monroe 1941); the version of the theory we use is based on the improved implementations developed by Ryzhov & Tareeva (1981) and Cerjan *et al*. (1985). All of these implementations use several approximations that compromise the accuracy of the predicted properties of the transition, e.g. the densities of the coexisting phases. Two of the approximations commonly used, namely carrying out the bifurcation analysis without the imposition of a required consistency condition and the representation of the pair contribution to the inhomogeneous singlet density distribution, are improved in the work reported in this paper.

Since we are interested in how the symmetry of the phase, with respect to which the quasi-two-dimensional liquid becomes unstable, changes as the thickness of the liquid increases, it is pertinent to consider the freezing transition in the limit that the fluid is confined to two dimensions. Considerable attention has been focused on understanding the phase diagrams of two-dimensional systems. At present, the Kosterlitz–Thouless–Halperin–Nelson–Young (KTHNY) theory (Kosterlitz & Thouless 1972, 1973; Halperin & Nelson 1978; Nelson & Halperin 1979; Young 1979) is the most widely accepted description of two-dimensional melting/freezing. This theory is based on a characterization of the two-dimensional solid as a deformable medium, with the inclusion of the two classes of point topological defects with smallest excitation energy to mediate structural changes; it relates the melting process to the mechanical instability of the two-dimensional solid. Although the theory allows for other possibilities, it is commonly taken to predict that a two-dimensional system that supports only one ordered solid phase melts via sequential continuous phase transitions. The first transition is from the solid with quasi-long-range positional order and long-range bond-orientation order to a phase with short-range positional order and quasi-long-range bond-orientation order, the hexatic phase. This transition is driven by the dissociation of bound dislocation pairs in the solid. The second transition transforms the hexatic phase into the liquid phase in which both positional and bond-orientation orders have short range; it is driven by the dissociation of individual dislocations to form disclinations. The formalism of the theory obscures the influence of the intermolecular interaction on the melting transition.

In its most common form, the density functional approach to the theory of two-dimensional melting cannot be used to determine whether a hexatic phase intervenes between the liquid and ordered solid phases of the system, owing to the mean-field character of the trial density fluctuations examined. Ryzhov & Tareyeva (Ryzhov & Tareyeva 1995; Pomirichi *et al*. 2002) have proposed a modification of the density functional approach to two-dimensional melting that combines macroscopic and microscopic elements, and which allows the two-particle distribution function to become anisotropic due to long-range bond-orientation correlation. They allow the coefficients in the Fourier expansion of the single-particle density in terms of the reciprocal lattice vectors of the ordered phase to vary slowly in amplitude and rapidly in phase. The latter variation, which depends on displacements of the particles from the lattice, is decomposed into two parts: one corresponding to a smoothly varying phonon field and the other a singular contribution associated with dislocations. The effect of dislocations is treated as in the KTHNY theory, i.e. using an elastic continuum model. The results of this analysis predict that the two-dimensional electron solid melts via two separate transitions with an intermediate hexatic phase, but that the two-dimensional hard-disc system melts via a single first-order transition without intervention of a hexatic phase.

Experimental verification of the KTHNY theory must contend with the fact that all physical realizations are at best quasi-two-dimensional. Alba-Simionescu *et al*. (2006) have reviewed the evidence for the influence of confinement on freezing and melting in molecular systems, with a brief mention of the colloid systems that are of interest to us. The available studies show that colloid particles that have magnetic dipole interactions are supported on a liquid surface freeze, as predicted by the KTHNY theory (Zahn *et al*. 1999). The situation is less clear for colloid particles with short-range interactions. Studies of those systems show the existence of a hexatic phase, but not the continuous fluid-to-hexatic and hexatic-to-solid transitions (Marcus & Rice 1996).

Verification of the KTHNY two-dimensional melting scenario via computer simulation has proven to be very difficult, primarily because both the orientation and translation correlation lengths become very large when the system is close to a continuous transition. For a system that contains hard discs, the best available evidence from computer simulation studies suggests that the liquid-to-hexatic transition is either continuous or very weakly first order, and that the hexatic-to-solid transition is continuous (Mak 2006). These conclusions are consistent with the results of a study of the convergence of the virial expansion for the density as a function of the fugacity, which imply that the hard-disc liquid-to-ordered solid transition is continuous, and with the results of a study of the elastic moduli of the hard-disc system, which suggest that the KTHNY continuous transition from the solid to a hexatic phase just barely pre-empts the solid-to-liquid first-order transition (Eisenberg & Baram 2006).

Now consider the case that the colloid system allows the particles limited out-of-plane motion. Does the change from two-dimensional to quasi-two-dimensional confinement alter the character of the solid-to-hexatic transition in a one-layer system? Moreover, as the plate separation defining the confined system increases further, when, in the transition from two to three dimensions, does the melting/freezing transition become first order? With respect to the first question, the issue to be resolved is whether small amplitude out-of-plane motion in a quasi-two-dimensional system generates only small quantitative corrections to the phase diagram predicted under the assumption that the motion is strictly two-dimensional, or whether it generates qualitative changes to that phase diagram. Zangi & Rice (1998) and Frydel & Rice (2003) have reported the results of extensive simulations of the several phase transitions in a quasi-two-dimensional system designed to mimic a real colloid assembly. The colloid–colloid interaction used in those simulations has an attractive component, with a range that is only a few per cent of the molecular diameter, a soft repulsive component of comparable range that connects the attractive well and a very steep repulsive interaction. The results of the simulations establish the occurrence of first-order liquid-to-hexatic and hexatic-to-solid transitions. At higher densities, they also reveal, consistent with the form of the colloid–colloid interaction, an isostructural solid I-to-solid II transition and a buckling transition, both of which are continuous. The dislocation pair, free dislocation and free disclination pair concentrations found in the simulations do not satisfy the predictions of the KTHNY theory. A comparison of the phase diagrams of the quasi-two-dimensional and two-dimensional systems and systems with the same interaction reveals the influence of out-of-plane motion on the locations of the phase boundaries and the ability of the system to support particular ordered arrangements of the particles. Overall, the most apparent macroscopic consequence of permitting out-of-plane motion is a shift of some of the phase boundaries to higher density. These shifts are most visible for the low-temperature–low-density portion of the liquid-to-solid I transition line, and for the low-temperature–low-density portion of the solid I-to-solid II transition line. As expected, the introduction of the out-of-plane motion increases the density range in which the liquid is stable. Simply put, at any selected density at which the liquid is stable, its entropy is increased by allowing out-of-plane motion. Similarly, at low temperatures, allowing out-of-plane fluctuations of the solid I structure increases its entropy and thereby extends the density range over which it is stable. Eventually, the density becomes so large that the puckered structure of solid II becomes stable. At higher temperatures, the liquid–solid I coexistence region is very nearly the same for the two-dimensional and quasi-two-dimensional systems. Moreover, a comparison of the density dependencies of the renormalized elastic moduli of the two-dimensional and quasi-two-dimensional systems reveals the following. For the two-dimensional system, at low temperatures, the KTHNY transition barely pre-empts the first-order melting transition from solid-to-hexatic, and, at higher temperatures, the ordinary thermodynamic first-order transition from the solid to the hexatic phase pre-empts the KTHNY transition between those phases. On the other hand, for the quasi-two-dimensional system at both low and higher temperatures, the first-order transition from the solid to the hexatic phase pre-empts the KTHNY transition between those phases. Clearly, the change from two-dimensional to quasi-two-dimensional geometry has an important effect on the character of the phase transitions in the system studied.

## 2. Theoretical background

We consider a one-component system. The first equation of the BBGKY hierarchy (Hansen & McDonald 1976*a*), relating the one-particle distribution function in the system to the pair interaction, is(2.1)and the equivalent LMBW equation is (Baus 1984*a*)(2.2)In equations (2.1) and (2.2), *ρ*(*r*_{1}) is the one-particle distribution function, *g*_{2}(*r*_{1}, *r*_{2}) is the pair correlation function, *c*(*r*_{1}, *r*_{2}) is the direct correlation function, and *v*(*r*_{1}, *r*_{2}) is the pair interaction. We also use the customary notation *β*=(*k*_{B}*T*)^{−1}. Equations (2.1) and (2.2) support solutions with uniform one-particle density distributions when the system density is less than the liquid density at the freezing point, *ρ*_{L}; for densities greater than the density of the crystal at the freezing point, *ρ*_{S}, the one-particle distribution is non-uniform. In the crystalline phase, the non-uniform one-particle distribution function is conveniently represented in the form(2.3)In equation (2.3), {** G**} is the reciprocal lattice vector set for the specified crystal. Baus has shown (Baus 1984

*b*) that to terms linear in Δ

*ρ*(

**), the two equivalent equations (2.1) and (2.2) take the following forms:(2.4)(2.5)in which(2.6a)(2.6b)In a uniform system, the pair correlation function can be represented in the form(2.7)in which is the sum of all distinct, at least doubly connected, diagrams consisting of two white 1-circles and**

*r**n*black

*ρ*-circles and

*f*-bonds (

*f*=e

^{−βv}−1); and

*δg*

_{2}(

*r*_{1},

*r*_{3})/

*δρ*(

*r*_{2}) is the sum of all distinct diagrams obtained by replacing a black

*ρ*-circle of

*g*

_{2}(

*r*_{1},

*r*_{3}) by a white 1-circle (Hansen & McDonald 1976

*b*). Substitution of equation (2.3) into equations (2.4) and (2.5) yields equations with the same functional form, namely(2.8)(2.9)with the definitions(2.10)The definitions of

*σ*

_{0}and

*σ*

_{G}are slightly different for the direct and the pair correlation function-based theories, respectively. In the former case,(2.11)whereas in the latter case(2.12)In both the direct and the pair correlation function-based theories, the representation used requires a sum over all reciprocal lattice vectors. Since we are primarily concerned with capturing the qualitative features of the alteration of lattice packing with plate separation in the confined system, albeit as accurately as possible, we retain in the sum only the smallest reciprocal lattice vectors. Although the direct and pair correlation function analyses are equivalent up to terms linear in Δ

*ρ*(

**), one or the other can be more difficult to implement in specific applications. For example, the direct correlation function is available for fewer systems than the pair correlation function, while the computation of**

*r**H*

_{L}(

*r*_{12}) requires the evaluation of an expansion in powers of the density of which only the first terms can be calculated.

In previous studies of freezing, the bifurcation point was found subject to the conditions(2.13)The first stated condition looks for a bifurcation point of equation (2.9), while the second condition looks for the smallest possible density for which the first condition is valid. Thus, only equation (2.9) is solved and one particular solution to (2.9), the one that has the smallest density value, is chosen. This procedure does not involve the use of equation (2.8), and it results in a large discrepancy between the theoretical and experimental values of *ϕ*_{0}. In the work reported in this paper, we find the bifurcation point of equation (2.9) that simultaneously satisfies equation (2.8); the results obtained are superior to those found using only equation (2.13).

## 3. Two-dimensional hard-disc freezing

To provide a reference for the change in lattice type with spacing of the plates in a quasi-two-dimensional system, we consider first the freezing of a two-dimensional hard-disc fluid (table 1). We have already noted that the latest and largest simulation has been interpreted as showing that the fluid-to-hexatic transition is weakly first order and the hexatic-to-crystal transition is continuous. The two-dimensional system isotherm appears to be converging to a limit with a flat liquid–hexatic coexistence region; the width of this coexistence region is very small (0.90<*ρ*<0.92; Mak 2006). We have also noted that the density functional approach to the theory of two-dimensional melting cannot be used to determine whether a hexatic phase intervenes between the liquid and ordered solid phases of the system owing to the mean-field character of the trial density fluctuations that are examined. A few words about those density fluctuations are in order at this point. Given the small density difference between the phases, the packing in the first few shells surrounding a particle in a hexatic phase is very similar to that in a crystal phase. Although the translation correlation function decays exponentially in the former case and algebraically in the latter case, at high density the difference is scarcely noticeable on the scale length of the first few shells. Moreover, the difference between the orientation correlation functions of the hexatic and hexagonal solid phases on this scale length is also very small. Accordingly, an analysis such as we present, which uses only very few reciprocal lattice vectors, will not distinguish between the hexatic and hexagonal solid phases, and we examine the stability of the hard-disc fluid with respect to the hexagonal solid. It is worth noting that a liquid phase that is unstable with respect to a hexatic phase is also unstable with respect to a hexagonal solid phase since the latter has a lower chemical potential than the former.

We first consider the predicted behaviour using the direct correlation function-based bifurcation theory. Baus & Colot (1987) have reported a very good approximation for the direct correlation function of a two-dimensional hard-disc liquid. It has the form(3.1)with(3.2)and *a*=*a*(*η*) obtained from(3.3)where is the Heaviside step function. Retaining only the first reciprocal lattice vector in the expansion (2.3), we find (Bagchi *et al*. 1983)(3.4)where ; ; *σ*_{G}=*ρ*_{S}*c*_{L}(** G**); and

*a*is the edge of a unit cell.

We then search for the bifurcation point,(3.5)of equation (2.9) and look for a particular solution that satisfies equation (2.8), which, in this case, becomes(3.6)The bifurcation point, determined by the same method as used in Cerjan *et al*. (1985), is found to be *σ*(*ρ*_{L}, *ρ*_{S})=0.84. This result and equation (3.6) define two equations for the two unknown quantities *ρ*_{L} and *ρ*_{S}. We used the Mathematica procedure *FindRoot,* which is an implementation of Newton's method, the secant method and Brent's method, to solve the simultaneous equations. We find that the freezing transition occurs when *ρ*_{L}=0.931 and that the solid density is *ρ*_{S}=0.934. The difference between liquid and solid densities is very small, so the predicted transition is almost continuous. The predicted density difference is somewhat smaller than that estimated from the simulations by Mak (2006), namely *ρ*_{S}−*ρ*_{L}=0.92−0.90=0.02. We have checked the adequacy of the one-reciprocal-lattice-vector representation by carrying out calculations with up to seven reciprocal lattice vectors. At that level of representation, we find *ρ*_{L}=0.902 and *ρ*_{S}=0.910, giving *ρ*_{S}−*ρ*_{L}=0.008. Given the simplicity introduced by the truncation, we consider the one-reciprocal-lattice-vector representation adequately accurate for our purposes.

We now consider the predicted behaviour using the pair correlation function-based bifurcation theory; it starts from equation (2.4). To the zeroth order of *ρ*, *H*_{L} can be dropped and *K*_{L} can be evaluated as follows. Let . Noting that , we have(3.7)Again, retaining only the terms up to the first reciprocal lattice vector, we findwhere *g*_{2}(1, *ρ*_{L}) is the liquid-state pair correlation function evaluated at hard-disc contact. We have calculated *g*_{2}(1, *ρ*_{L}) using the Pade approximant to the equation of state proposed by Ree & Hoover (1963). The fluid-to-solid transition is found to occur when *ρ*_{L}=0.655 and *ρ*_{S}=0.686.

We now evaluate *H*_{L} using equation (2.6*b*) and the leading two terms in equation (2.7). Setting *Φ*=arccos(*R*_{12}/2), the term linear in *ρ* is(3.8)which leads to(3.9)The term that is second order in *ρ* consists of five diagrams (Hansen & McDonald 1976*c*). Writing , they have the forms(3.10)The terms *H*_{L2D} and *H*_{L2E} were evaluated numerically. The sum of the other three diagrams iswhereUsing the terms linear and quadratic in *ρ*, we construct a Pade approximation to *H*_{L},(3.11)Approximating *H*_{L} with *H*_{L1} (first order in *ρ*), *H*_{L1}+*H*_{L2} (second order in *ρ*) and *H*_{L,Pade} (Pade approximation), respectively, we search for the fluid-to-solid transition. The first- and second-order approximations have little effect on shifting the liquid and solid transition densities. When equation (3.11) is used, we find that the transition occurs when *ρ*_{L}=0.85 and *ρ*_{S}=0.87. These liquid and solid densities are not as close to those estimated from the simulations by Mak as are the densities derived from the direct correlation function theory, but they share the characteristic that *ρ*_{S}−*ρ*_{L} is very small.

## 4. Thickness dependence of quasi-two-dimensional freezing

As described in §1, a hard particle system confined between two smooth parallel plates with a small separation undergoes the transitions as the plate separation increases and the density is varied at a fixed plate separation. We now examine the thickness dependence of the transition from the liquid to the solid, in particular the change in packing structure of the solid with respect to which the liquid is unstable. The input required for the bifurcation analysis is the thickness dependence of the direct or the pair correlation function. We calculate approximations to these functions by mapping the system confined to a slab with height *Hσ*=(1+*h*)*σ* onto the two-dimensional system, i.e. we define an effective diameter and use that diameter in the known two-dimensional functions (Schmidt & Lowen 1997). For small *h*, this effective diameter is *σ*_{eff}/*σ*=1−(*h*^{2}/6). Thus, the *h*-dependence of the direct correlation function is taken to be(4.1)and the *h*-dependence of the pair correlation function is taken to be(4.2)

It is intuitively evident that the approximation described is plausible when *h* is very small, but must become increasingly inaccurate as *h* increases. With this observation in mind, we first calculate the shape of the boundary between a fluid and a single-layer hexagonal crystal as a function of *h*. Equation (4.1) is inserted in equation (2.11) and we search for the bifurcation point as a function of *h*; the results are displayed in figure 1. The shapes of the predicted phase boundaries are very similar to those found in the simulations, but they are displaced to a higher density. Clearly, the theory gives a qualitatively sound description of the thickness dependence of the one-layer freezing transition, with semi-quantitative accuracy. Given the approximations employed, we consider the agreement between theory and simulations to be good.

With the same approximations to *H*_{L} as described above, equation (4.2) is inserted in equations (2.6*a*), (2.6*b*) and (2.12) and we search for the bifurcation point as a function of *h*; the results are displayed in the lower parts of figures 2 and 3. In the lower parts of these figures, the solid lines refer to the fluid→1▵ transition and the dotted lines to the transition. The theory, either with second order in density corrections or a Pade approximant, correctly predicts that the fluid→1▵ transition pre-empts the transition for *h*<0.6. Although the shapes of the phase boundary predicted using the Pade approximant are reasonably accurate, those predicted using the theory with second order in density corrections are not as good as those obtained from the direct correlation function analysis. We note that the predicted phase boundaries using the Pade approximant are displaced to a lower density than the simulation results by a small amount, whereas those predicted using the second order in density corrected theory are shifted to a lower density by a considerable amount. All in all, the theory based on the pair correlation function does capture the essential features of the *h*-dependence of the fluid→1▵ transition.

When considering the liquid-to-two-layer crystal phase transition, our approximation requires the use of two effective diameters: one for the interaction between particles in the same layer (*σ*_{same}=*σ*) and the other for the interaction of particles in different layers . This very crude approximation does not work when used with the direct correlation function theory; no meaningful bifurcation point can be found. In contrast, when used with the pair correlation function theory, we obtain qualitatively valid predictions that capture the essence of the physics of the liquid-to-two-layer-ordered solid transition.

Suppose the ordered solid phase has two layers of width 2*ϵ*, where *ϵ* is small. To the zeroth order in *ρ*, we need only *K*_{L}. We assume that, for particles 1 and 2 within the same layer,(4.3)where *g*_{2,same}(1) is the pair correlation function at contact at the effective diameter for particles in the same layer, , and −*ϵ*<*z*<*ϵ*. When particles 1 and 2 are in different layers,(4.4)where *g*_{2,diff}(1) is the pair correlation function at contact for particles in different layers. We now expand the singlet distribution function in a Fourier series keeping only one reciprocal lattice vector. The result is(4.5)and the nonlinear equation takes the form as follows.

For particles 1 and 2 in the upper layer and particle 3 in the lower layer,

(4.6)

For particles 1 and 3 in the lower layer and particle 2 in the upper layer,

(4.7)Thus, we have(4.8)where(4.9)Then,(4.10)and(4.11)

Define *Π*_{trans}=−*A*^{−1}(∂*F*/∂*h*) and *Π*_{lat}=−*h*^{−1}(∂*F*/∂*A*), where *A* is the area of the plate and *F* is the free energy. The *h*-dependence of the pair correlation function is , as we now have two layers. The transverse compressibility function and the lateral compressibility function can be written as(4.12)(4.13)That is(4.14)(4.15)Substitution of equations (4.14) and (4.15) into equations (4.10) and (4.11), followed by a search for the bifurcation point as a function of *h* with the same approximations to *H*_{L} as described above, leads to the results displayed in the upper parts of figures 2 and 3. In the upper parts of these figures, the solid lines refer to the transition and the dotted lines to the fluid→1▵ transition, i.e. the dotted lines are the extensions of the fluid→1▵ transition line for *h*>0.6. The theory correctly predicts that the transition pre-empts the transition for *h*>0.6. Although the shapes and positions of the predicted phase boundaries are not in good agreement with the experiment, the theory does capture the essential feature of the simulation results, namely that the symmetry of the solid phase to which the fluid freezes changes from hexagonal to square and the number of layers from one to two when *h*>0.6.

## 5. Discussion

The results reported in this paper demonstrate that the analysis of the bifurcation of solutions of the first equation of the BBGKY hierarchy or the equivalent LMBW equation provides a qualitatively satisfactory interpretation of the dependence on wall separation and density of the fluid→1▵ and the transitions in a quasi-two-dimensional colloid system. Our analysis cannot answer whether or not a hexatic phase intervenes between the two-dimensional liquid and ordered solid phases or whether the said hexatic phase persists as the two-dimensional system is transformed to a quasi-two-dimensional system by increasing *h*. The analysis also does not reveal the mechanism by which the one-layer hexagonal lattice rearranges to a two-layer square lattice.

Some insight into this dramatic change in lattice symmetry is provided by the results of simulations reported by Zangi & Rice (2000), for the case of a two-layer quasi-two-dimensional suspension of hard spheres constrained to have thickness *h*=0.8. They studied the transition over the density range of 0.9300<*ρ*<1.0300. When *ρ*=1.0300, the two-layer-ordered solid is hexagonal, and when *ρ*=0.9300 the ordered solid is square. However, when *ρ*=1.0200, there is a rhombic phase whose smaller lattice angle, *θ*_{1}, increases with decreasing density until *ρ*=0.9700, where *θ*_{1}=75°, signalling the onset of coexistence between the rhombic phase and the phase. Both the linear and zigzag rhombic phases are observed to occur. The lateral pressure as a function of the one-layer two-dimensional number density displays a van der Waals loop in the density range of 0.9500<*ρ*<0.9800; this density range matches that for coexistence between the two-layer-rhombic and the phases, which clearly identifies the first-order character of that transition. Zangi & Rice could not clearly determine the character of the 2▵→two-layer-rhombic transition, which spans a very small density range; they suggest that it is either weakly first order or continuous.

If we suppose the rhombic phase to be generated by the distortion of the hexagonal lattice via particle displacement, the structure of the distortion can be realized by connecting the six closest neighbours of a central particle. The structure of the distortion for the zigzag rhombic phase is shown in the lower right of figure 4, reproduced from the paper by Zangi & Rice. Note that two of the six neighbours are no longer nearest neighbours, and the nearest neighbour configuration is kite shaped (lower left). The corresponding distortion for the linear rhombic phase is shown in the upper part of figure 4, where the shape of the arrangement of the four nearest neighbours is rectangular. The results obtained from the Zangi–Rice simulations show that for the pure rhombic phase, 1.0200<*ρ*<0.9800, the lateral pressure isotherm does not show a van der Waals loop and there is only one set of values for *θ*_{1} for any density; hence, in agreement with the suggestion by Schmidt & Lowen, the linear and zigzag rhombic phases are degenerate and have the same lattice angle.

We conjecture that the change in lattice symmetry from 1▵ to has a similar genesis. Noting that as *h* increases, the density of the 1▵ lattice decreases, distortions of the hexagonal packing of the type just described become easier, and ultimately the competition between the free volume available to two square-packed layers with smaller overall thickness than one hexagonal-packed layer, and the free volume associated with more efficient packing within the hexagonal-packed layer than in a square-packed layer, favours the former.

## Acknowledgments

This work was supported by the National Science Foundation funded MRSEC Laboratory at The University of Chicago.

## Footnotes

- Received July 3, 2007.
- Accepted September 19, 2007.

- © 2007 The Royal Society