Frustration of freezing in a two dimensional hard-core fluid due to particle shape anisotropy

The freezing mechanism suggested for a fluid composed of hard disks [Huerta et al., Phys. Rev. E, 2006, 74, 061106] is used here to probe the fluid-to-solid transition in a hard-dumbbell fluid composed of overlapping hard disks with a variable length between disk centers. Analyzing the trends in the shape of second maximum of the radial distribution function of the planar hard-dumbbell fluid it has been found that the type of transition could be sensitive to the length of hard-dumbbell molecules. From the ${NpT}$ Monte Carlo simulations data we show that if a hard-dumbbell length does not exceed 15% of the disk diameter, the fluid-to-solid transition scenario follows the case of a hard-disk fluid, i.e., the isotropic hard-dumbbell fluid experiences freezing. However, for a hard-dumbbell length larger than 15% of disk diameter, there is evidence that fluid-to-solid transition may change to continuous transition, i.e., such an isotropic hard-dumbbell fluid will avoid freezing.


Introduction
Since Alder and Wainwright [1] in 1962 showed by molecular dynamics simulations that two-dimensional (2D) fluid of circular hard-core species -hard disks -can freeze, a number of papers have been published studying the properties and physics behind such a phenomenon. The hard-disk fluid is one of the simplest interaction models in soft condensed matter science. This model is far from being trivial as might be thought from a first look and should be treated in the same way as some other (classical) twodimensional models such as the Ising model, X Y model, and other lattice models. Species of the hard-disk model (e.g., atoms, molecules or particles) do not experience any other interaction except being prohibited from mutual overlap. This apparent simplicity as well as a desire to understand the mechanism that enables the system with no thermal activation to become unstable, undergoing a freezing transition, are the main driving forces of the interest to this model fluid. The most recent advances in understanding the physics behind the freezing transition in a hard-disk fluid concern the large scale computer simulation data reported by Zollweg and Chester [2], Mak [3] and Bernard and Krauth [4]. According to these simulation studies, the hard-disk fluid becomes solid in two steps: (i) by means of the first-order type of transition from an isotropic fluid phase to a hexatic phase and (ii) continuously from a hexatic phase to a solid.
The monolayers of colloidal particles that very often serve as a prototype of two-dimensional and quasi-two-dimensional systems (e.g., see review article [5] and references therein) still in various ways (the size polydispersity, shape anisotropy, extra interactions, etc) differ from the "ideal" hard-disk system. Some of these features, namely, the extra out-of-core interaction between disks as well as the disk size bidispersity, have been investigated [6,7] and their impact on freezing transition has been revealed. In particular, it has been shown [7] that in a binary equimolar mixture of hard disks, there is a limiting large-to-small disk diameter ratio of around 1.2 when freezing transition is still localized, while the same binary hard-disk mixture having the diameter ratio of around 1.4 does not exhibit the freezing behavior.
In the present study we aim to explore how the freezing transition in a hard-disk fluid might be affected by the shape anisotropy imposed on hard-disk species. In the following section 2 we start with the introduction of the shape anisotropy into a hard-disk model and the description of the necessary information regarding the details of computer simulations. The full set of the results for the present study and discussions are collected in section 3, while conclusions are presented in a final section 4.

Modelling and details of simulations
A simple way to introduce shape anisotropy to the planar hard-disk species is to consider the planar hard dumbbell-like objects that consist of two fused hard disks of diameter σ and their centers at a distance d σ. Such an approach permits to have the hard-core species elongated in one direction with aspect ratio (σ + d σ)/σ = 1 + d , where d plays the role of the anisotropy parameter, assuming values 0 d 1. When d = 0, the hard-disk fluid is recovered, while the largest anisotropy that can be reached for this hard-core model, d = 1, corresponds to the case of tangent hard disks and will be referred to as a harddimer fluid. Besides the anisotropy parameter, the system is characterized by packing fraction η defined as η = ρ A d , where ρ = N /A and A is the area of the system, A d is the area of a hard-dumbbell molecule that depends on the anisotropy parameter d , while N is the number of hard-dumbbell species. Figure 1 (the left-hand-side frame) shows the hard-dumbbell modelling as well as the phase diagram in anisotropy-density coordinates that was suggested by Wojciechowski [8] based on heuristic reasoning   [8], and schematic of a single hard-disk molecule of diameter σ, and hard-dumbbell molecules that consists of two fused hard disks of diameter σ and centers at a distance d σ. Right: Monte Carlo simulations data for the radial distribution function of a hard-disk system at two packing fractions η: before freezing (dashed line) and before melting (solid line). Inset: Sketch of the disk regular hexagonal configurations where four disks filled with black and two empty disks drawn by the dashed line are the first coordination shell neighbors of the central disk A, while those six filled with gray and those six empty disks drawn by the solid line belong to the second shell neighbors. using free-volume arguments. According to this diagram, any hard-dumbbell system will always undergo a first-order transition from an isotropic fluid phase to a solid phase. The latter could be of various types such as an ordinary crystal, a plastic crystal (PC) or a disordered crystal (DC) depending on the magnitude of the anisotropy parameter d . As it has been already mentioned in the Introduction, a scenario that involves first-order transition indeed takes place in the limit d = 0 case, i.e., for a hard-disk fluid [2][3][4].
The transition with two coexisting phases -isotropic fluid and disordered crystal phase -has been also confirmed for another limiting case d = 1, i.e., for a hard-dimer system [9], as well as for a hard-dumbbell fluid with d < 1, namely, in the case of d = 0.924 [10].
In the present study, we consider the features of fluid-to-solid transition in a hard-dumbbell fluid that is characterized by a small deviation of molecule shape anisotropy from a hard-disk fluid, by assuming the values of the anisotropy parameter within the range 0 d < 0.28. The latter was chosen by means of the anisotropy-density phase diagram in figure 1 that associates this range of a hard-dumbbell fluid anisotropy with the coexisting plastic crystal solid.
We have used Monte Carlo simulations technique to study the properties related to freezing transition in a hard-dumbbell fluid. All simulations were performed in a rectangular box with an aspect ratio of 3/2 with the usual periodic boundary conditions. Two different thermodynamic ensembles have been employed. Namely, the constant-density NV T ensemble was used to evaluate the radial distribution function and global bond-orientational order parameter of the hard-dumbbell fluid, while the constant-pressure N pT ensemble was used to discriminate between a first-order and continuous transition in the same system. The symbols N ,V, T and p specify the number of particles, volume (area in two dimensions), temperature and pressure used in the simulations, respectively.
In the course of our simulations, the single molecular move (in what follows is referred to as the iterative step) consists of a random positional displacement of molecular mass center accompanied by a rotation of molecular axes. For denser states, we also attempted a π/3 orientational displacement of molecular axes that may be commensurable with a crystalline solid phase. The N pT simulations additionally consist in the change of the area of the simulation box. The acceptance ratio was maintained between 20% and 30% by adjusting the maximal size of positional displacement and maximal amplitude for the rotation in the case of NV T ensemble, and between 20% and 50% by adjusting the maximal size of positional displacement, maximal amplitude for the rotation and maximal variation of the box size in the case of N pT ensemble.
A standard Metropolis algorithm was used to obtain the ensemble averages. For NV T simulations, each equilibration run was relaxed for at least 10 4 iterative steps. The resulting data for radial distribution function and global bond-orientational order parameter were obtained by averaging at least over 600 configurations, each being relaxed by 10 3 iterative steps. Whereas for the N pT simulations, each equilibration run was relaxed for at least 10 6 iterative steps, the distribution of densities was accumulated after at least 2 × 10 7 iterative steps to obtain a good sampling. Various sizes of system were investigated, ranging from 100 to 1000 molecules. Similar to the case of a hard-disk fluid [11,12], no systematic system size effects were observed in the NV T data for the second maximum of the radial distribution function, and subsequent calculations were made using N = 400 hard-dumbbell molecules. These simulations have been used to collect data for the global bondorientational order parameter. However, in this case it is well documented [13] that N = 400 is too small to identify the boundary of freezing and melting densities. Therefore, the data for global bondorientational order parameter could be used for qualitative purposes only.
The Monte Carlo simulations in the N pT ensemble are efficient at identifying the type of fluid-tosolid transition. According to the studies by Lee and Strandburg [14] on the application of the isobaric Monte Carlo simulations, the system size 100 N 400 is sufficiently large to identify the fluid-to-solid transition in a hard-disk system as first order transition. Therefore, the N pT calculations of the equation of state were made using N = 100 hard-dumbbell molecules. The consequences of such a system size for pressure magnitude of the hard-disk fluid in the vicinity of transition region are already discussed in the literature [15].

Results and discussions
Following the observation by Truskett et al. [11], the fluid composed of planar hard disks (the case of the anisotropy parameter d = 0) exhibits a structural precursor to freezing transition that manifests itself as a shoulder in the second maximum of the disk-disk radial distribution function. In fact, an increase of the density in a planar array of N hard disks of diameter σ, that are uniformly spread on area A, is accompanied by formation of the local quasi-regular hexagonal arrangement. By quasi-regular hexagons we mean the sixfold configurations (see the right-hand side frame of figure 1 adopted for details from reference [12]) that, at certain disk density or packing fraction η = πN σ 2 /(4A), is characterized by a constant gap width between each pair of the neighbor disks [12]. Here, η cp = π/(2 3) ≈ 0.907 is the disk close packing (CP) fraction. Such a sixfold formation implies that center-to-center distances r between disks and, consequently, the gap width ∆ shortens under the density increasing.
By analyzing the schematic drawing in the inset of figure 1, Huerta et al. [12] pointed out that there exists a packing fraction η cage = π 3/8 = 0.680 and, consequently, a gap width ∆/σ = 2/ 3 − 1 ≈ 0.155, which in transparent way demonstrates the existence of two kinds of the second shell neighbors discerned by means of the distance criterion. Namely, the second shell neighbors on the distances (e.g., AD, . . . ) that are always larger than 2σ, and those on the distances (e.g., AB, AC, . . . ) that could be shorter than 2σ. In fact, the existence of these two kinds (larger and shorter) of the second shell neighbors is the origin, initially of a shoulder and then, with an increase of the density, a split of the second peak of the disk-disk radial distribution function g (r ) for distances around r /σ = 2.
On the other hand, it implies that the center-to-center distances r between the first-shell alternating neighbors of any disk in a fluid system become, on average, shorter than 2σ, i.e., the set of only three firstshell alternating neighbors is capable of forming a cage, preventing the central disk to wander. Thus, the hard-disk fluid will exhibit a tendency to freeze. The corresponding hard-disk fluid caging density, η cage = 0.680, is lower than the freezing transition density, η f = 0.700, that follows from large scale computer simulation data [4]. Actually it means that in practice disks are not spread uniformly and some of the center-to-center distances (AB, AC, . . . ) become shorter than 2σ already for the disk packing fractions η 0.6, initiating deformation of the smooth second maximum of the disk-disk radial distribution function g (r ).
Therefore, a starting point to discuss the effect of shape anisotropy on the freezing transition in a fluid composed of planar hard dumbbells will be to analyze the trends in the second maximum of the center-center radial distribution functions shown in figure 1 for packing fractions η = 0.686 and 0.723 upon increasing the value of anisotropy parameter d . As expected, and in agreement with a prediction of the density-anisotropy phase diagram in figure 1, for small values of the anisotropy parameter, d = 0.05, the freezing scenario seems to be quite similar to the case of a hard-disk fluid. Namely, for the case of d = 0.05 we are observing practically the same shape of the second maxima for both densities with a tiny shift in the positions to slightly larger center-to-center distances r . As the shape anisotropy increases to d = 0.15, the shift in the positions of maxima and minima of the g (r ) becomes more pronounced. Although for packing fraction η = 0.686, the shoulder in the second maxima is still noticed but it is less evident. However, since the split of the second maxima is definitely

23605-4
Frustration of freezing in a two dimensional hard-core fluid present, one could expect that the freezing mechanism that was discussed for a hard-disk fluid is still at work.
However, the case of the shape anisotropy parameter d = 0.27 seems to be crucial from the point of view of the freezing transition. The resulting radial distribution function in such a hard-dumbbell fluid, shown in figure 4, is totally different from the one of the hard-disk fluid for both packing fractions considered. Most important is that there are no indications either of a shoulder at lower packing fraction or of a split of the second maximum at higher packing fraction. When discussing figure 1, we have already pointed out that besides serving as a precursor towards approaching the freezing transition, the presence of these features is signaling about the formation of the local quasi-regular hexagonal arrangement.
A quantitative measure of hexagonal order in the system is provided by the so-called global bondorientational order parameter Ψ 6 , that was evaluated during the NV T Monte Carlo simulation runs using definition where j runs over all disks in the system, k runs over all geometric nearest neighbors (nn) of disk j , each obtained through the Voronoi analysis and N nn is the total number of such nearest neighbors in the system. The angle θ j k is defined between some fixed reference axis in the system and the vectors bonds connecting nearest neighbors j and k. As one can see from figure 5, indeed, the global bond-orientational   show the N pT Monte Carlo data for pressure βp, where β = 1/k B T and k B is the Boltzmann constant, in a wide range of packing fraction η (the left-hand-side frame of figure 6). The hard-dumbbell systems presented in figure 6 are the same as have been discussed in figures 2-5 to learn about the trends in the features of the second maximum of the radial distribution function. As expected, those features that are sensitive to the shape anisotropy in the density region that precedes fluid-to-solid transition; for better visualization, this region is shown separately on the right-hand-side frame of figure 6. On the other hand, from the results in the left-hand-side frame of figure 6 it follows that in the limiting case of low (η < 0.3) as well as at high (η ≈ 0.8) packing fractions, all pressure isotherms tend to coincide, showing no dependence on the anisotropy parameter.
Although quantitatively the results for pressure p and for the freezing and melting transition densities could be affected by the system size effects, which is rather evident in the case of a hard-disk fluid,   [4]; the filled squares at the top represent computer simulation data for a hard-disk fluid that are corrected for finite-size effects [2,15,16].
the qualitative picture must be correct. Namely, like in the previous N pT Monte Carlo simulations [14] as well as in quite recent large-scale NV T Monte Carlo simulations [4,17,18]  Regarding the issue of the system size effects in the case of a hard-dumbbell fluid, we note that contrary to a hard-disk system, all the existing computer simulation data for a hard-dumbbell system have been obtained with N = 112 molecules [8][9][10].
In fact, the statement regarding the frustration of freezing in a hard-dumbbell fluid with the anisotropy parameter d = 0.27 is consistent with the information that follows from the data presented in figure 7 where density distributions obtained in the N pT Monte Carlo simulation are shown. Each curve

23605-7
for a particular system in figure 7 corresponds to a constant pressure p, and the position of its maximum determines the density (packing fraction η) of the stable phase of this system in this thermodynamic state. The curves with two maxima in the case of d = 0 and d = 0.15 each tells us that at pressure p, that corresponds to this curve, there is coexistence of two phases and, consequently, the fluid-to-solid transition is discontinuous or of first order. By contrast, in the third case of hard-dumbbell fluid with d = 0.27, all curves show only one maximum, and under pressure increasing, the position of this maximum is moving towards higher densities (packing fractions η). Thus, the fluid-to-solid transition in this case is continuous.

Conclusion
Contrary to the hard-disk system, where fluid-to-solid transition has been the subject of a long-standing debate [19], the hard-dumbbell system has not been explored in full yet. This fact looks rather strange in view that almost three decades ago the theoretically predicted phase diagram of the harddumbbell system in full range of the anisotropy parameter, 0 d 1, was reported based on heuristic reasoning and using free-volume arguments [8]. The main conclusion of that study was that in the full range of the anisotropy parameter values, fluid-to-solid transition is of the first order. Up to date, this result was confirmed by computer simulations for two particular values of the anisotropy parameter, i.e., d = 1 [9] and 0.924 [10]. Obviously, theoretical predictions are valid for the case of d = 0, i.e., for a hard-disk fluid [3,4,13,17,18].
In this study, we turned our attention to the case of anisotropy parameter values range 0 d 0.28 that was indicated in that theoretical analysis [8] as a unique one where the isotropic fluid phase freezes by coexisting with the plastic crystal solid. First of all, using the NV T Monte Carlo simulations data we have shown that such a feature of the radial distribution function as the shoulder in its second maximum, that was suggested by Truskett et al. [11] as the structural precursor to freezing transition in the case of a hard-disk fluid, is still preserved for a hard-dumbbell fluid with anisotropy parameter values up to d 0.15. However, in the case of a hard-dumbbell system with anisotropy parameter d = 0.27, the radial distribution functions do not show any sign of the shoulder in the range of distances that correspond to the second maximum. Secondly, by performing the N pT Monte Carlo simulations we obtained that indeed first-order freezing transition from an isotropic fluid to denser highly ordered phase persists in the case of a hard-dumbbell fluid with an anisotropy parameter d 0.15 while two-phase coexistence is absent in the case of d = 0.27.
Although in the present study we did not consider this question, the coexisting denser phase in the case of a hard-dumbbell fluid, with anisotropy parameter d 0.15, could be associated with a plastic crystal, as it was predicted [8]. We remind the reader that the plastic crystal phase represents one where densely packed molecules are still free to rotate. In fact, close examination of the hard-disk fluid configurations in the vicinity of freezing transition revealed the formation of a local quasi-regular hexagonal arrangement of hard disks [12]. Under increasing density, the average center-to-center distance between neighboring disks continuously shortens and, at the point of freezing. This results in a disk spacing of approximately 15% of disk diameter. Evidently, any disturbance that will tend to violate such a spacing, e.g., due to the attractive or repulsive interaction between disks, or the degree of disk size polydispersity, etc., could potentially result in the frustration of the freezing phenomenon. This is what we have observed in the present study: if the disk elongation starts to exceed 15% of the disk diameter (the anisotropy parameter d = 0.27), the freezing transition cannot survive and under an increase of the density, the disordered isotropic hard-dumbbell fluid undergoes a continuous transition to a denser and more ordered phase. Such a mechanism has already proved to work in the case of hard disks with square-well attraction [6] and for an equimolar binary hard-disk mixture [7].