Solid-solid collapse transition in a two dimensional model molecular system

As the form of the Mercedes-Benz potential is written as a completely separable isotropic and anisotropic parts having distinct ground states (triangular and honeycomb), one would expect that the interplay between the relative strength of the isotropic to anisotropic interaction can give rise distinct solid phases – either stable or metastable. Recently, we find that the interplay between the strength of isotropic and anisotropic interactions gives rise to rich phase diagram21 (). The computed phase diagram consists of isotropic liquid and three crystalline phases – honeycomb, oblique and triangular. Interestingly, we note that the oblique phase is neither the ground state of isotropic nor the anisotropic part of the potential. It forms due to delicate balance between the two types of interactions. In order to study the solid-solid collapse transition from low density open structured honeycomb to high density oblique phase the honeycomb solid (ground state of the anisotropic part of the potential) should be metastable (as shown in Fig. 1A) at a given temperature, pressure and relative strength of isotropic interaction (${\lambda }_i$). We have constructed such a free energy surface by tuning the relative strength of the isotropic interaction (${\lambda }_i$).

The honeycomb lattice is stabilized by 3 hydrogen bonds (HBs) per particle and the oblique phase is stabilized by two HBs per particle along with LJ interactions (see Fig. 1B). Alignment of bonds in Fig. 1B indicates hydrogen bonding. Since both the phases are solid and stabilized by anisotropic interactions (either completely or partially), one would expect that entropy of two phases will not differ significantly. Only vibrational entropy (and some contribution from mobility of defects) is present and also not expected to differ significantly as both the solid phases are stabilized by anisotropic interactions, though, the fraction of directional bonding decreases on going from honeycomb to oblique phase. This indicates that entropy has very little role (compared to gas-liquid and liquid-solid transitions) in the construction of free energy surface for solid-solid transitions. Thus the temperature change should not alter significantly the free energy surface; however, it is important to overcome the required free energy barrier for the transition.
To support the proposed free energy surface (Fig. 1A) in Fig. 1C, we have shown the temperature dependence of enthalpy, energy and relative contribution of the isotropic (LJ) and anisotropic (HB) interactions to the total energy per particle of the system for isotropic to anisotropic strength ratio, ${\lambda }_i$ = 0.4. The starting configuration is honeycomb lattice at pressure $P$ = 0.10 and $T$ = 0.05, which on heating collapses to the high density oblique phases. Although, in Fig. 1C we have shown the dependence of energy and enthalpy on temperature near the transition temperature only, we find that very weak dependence persists up to much lower temperature $T$ = 0.05. Note that, contrary to the ordinary solid-liquid transition (where total energy of the system increases), in this case the total energy of the system decreases. A very weak dependence of enthalpy on temperature and decrease in the enthalpy of the system on transition supports the free energy surface proposed in Fig. 1A. Honeycomb solid is metastable at given initial pressure and temperature ($T$ = 0.05 and $P$ = 0.10). We also note that the relative contributions of LJ and HB parts of interaction potential to the total energy per particle is switching. This indicates the increased importance of isotropic interaction in the oblique solid phase.

In order to gain quantitative insight into the structures of honeycomb and oblique phases, in Fig. 2 we have shown the computed radial distribution functions and bond length distributions of nearest neighbors for both honeycomb and oblique solid phases near transition temperature. We have used the Voronoi construction in order to identify the nearest neighbors ($NN$) of each particle. Note the splitting of first peak of radial distribution function and also the bimodal bond length distribution of nearest neighbors for the oblique solid phase. This signifies that in the oblique solid all nearest neighbor particles are not equidistant. There are two length scales present in the system, which is a consequence of the two length scale potential. Particles having lower bond length are stabilized by isotropic LJ interaction and larger bond lengths are stabilized by anisotropic HB interaction. Contrary to the oblique solid, honeycomb solid is stabilized only by anisotropic HB interactions and thus the nearest neighbor distance distribution has only one peak centered at HB length and there is no splitting in the first peak of radial distribution function.
In the subsequent sections, first, we shall discuss the temperature induced transition followed by a comparative study of pressure and temperature induced transitions.

In Fig. 3, we have shown the order parameters change and their susceptibilities for temperature induced collapse transition. On increasing temperature, thermal fluctuations destabilize the honeycomb structure and collapses to a high density oblique structure with large change in density (from $\sim 0.79$ to $\sim 1.56$) (see Fig. 3A). In the collapsed transition, high density change indicates that the transition can only be first order. As due to increase in order parameter difference the surface energy cost for formation of new phase also increases. Since the barrier in the global free energy surface as function of global order parameter (density and order) is related to the surface tension, the probability of long wavelength fluctuations of the global order parameters (density fluctuation is coupled with order fluctuation) in the system is less. In Fig. 3B, we have reported the change of both the global six-fold bond orientational order (${\psi }_{6g}$) and local three-fold bond orientational order (${\psi }_{3l}$) with temperature. m-fold global bond orientational order (${\psi }_{mg}$) is defined as

and m-fold local bond orientational order (${\psi }_{ml}$) is defined as

where m = 6 for six-fold and m = 3 for three-fold symmetry. NN is the number of nearest neighbors, ${\theta }_{ij}$ is the angle of the bond that the j${}^{th}$ neighbor makes with the tagged i${}^{th}$ particle. Significantly large values of global order parameter,${\psi }_{6g}$, in both phases indicate that there is a global translation of either six-fold or three-fold bond orientational order. As ${\psi }_{6g}$ does not distinguish between the six-fold and three-fold bond orientation order, we have also plotted the local three-fold bond orientational order (note that global three-fold bond orientational order is zero for both six-fold and three-fold symmetry). However, the later does not guarantee the long range translation of local bond orientational order. Sharp transition in ${\psi }_{3l}$ and significantly large values of ${\psi }_{6g}$ in both the phases indicate a one step transition from three-fold honeycomb lattice to six-fold (distorted) oblique solid phase.

In Fig. 3C and 3D, we have plotted the susceptibilities along density (isothermal compressibility, ${\kappa }_T$) and global six-fold bond orientational order parameter $\left({\chi }_6\right)$. Six-fold bond orientational order susceptibility is defined as ${\chi }_6=⟨|{\psi }_{6g}{|}^2⟩-⟨|{\psi }_{6g}|{⟩}^2$ and and isothermal compressibility is defined in terms of volume fluctuations as ${\kappa }_T=\left(⟨V^2⟩-⟨V{⟩}^2\right)/T⟨V⟩$. Weak temperature dependence and finite jump in both ${\kappa }_T$ and ${\chi }_6$ suggest the strongly 1${}^{st}$ order transition.

In order to find the pathway of transition first we need to distinguish the particles belonging to new phase in the sea of metastable bulk honeycomb solid phase. To distinguish the nucleating and growing new oblique-like particles from the parent metastable honeycomb phase, we computed the distribution of local 3-fold bond orientational order (${\psi }_{3i}$), defined as ${\psi }_{3i}=\left|{\sum }_{j=1}^{NN}{e}^{i3{\theta }_{ij}}/{\sum }_{j=1}^{NN}{e}^{i3{\theta }_{ij}}NNNN\right|$, for both the solid phases. As shown in Fig. 4, the complete separation of two distributions

enables us to separate the growing new phase from the parent honeycomb solid phase. Due to significant overlap in the distribution of local six-fold bond orientational order we have not used it as a parameter to distinguish the two phases. Note that, though, ${\psi }_{3i}$ is able to separate the oblique solid from the honeycomb solid phase, it does not distinguish the oblique solid-like and liquid-like particles (as the peak of distribution is centered at very low value of ${\psi }_{3i}$). A two order parameter criteria using the values of both ${\psi }_{3i}$ and ${\psi }_{6i}$ is required to separate the oblique phase with liquid-like particles as the oblique solid is characterized with low ${\psi }_{3i}$ and high ${\psi }_{6i}$ values (see Fig. 3B).
In Fig. 5A - 5D, we have shown the snapshots of the system at different stages during the collapse transition. Depending on ${\psi }_{3i}$ value we have defined three types of particles in the system. The particles having ${\psi }_{3i}$ value greater than 0.75 is considered as honeycomb type and indicated with black color, red color particles have intermediate order with ${\psi }_{3i}$ value between 0.25 and 0.75. The particles having ${\psi }_{3i}$ value less than 0.25 and ${\psi }_{6i}$ value greater than 0.2 (lower cut-off value of ${\psi }_{6i}$ below which the distribution of ${\psi }_{6i}$ $\left(P\left({\psi }_{6i}\right)\right)$ in oblique phase is zero) are considered as oblique-like. Note that these criteria do not distinguish oblique-like and triangular solid-like particles. However, simplicity of 2D systems allows us to distinguish these two types of orders by visual inspection of snapshots. Recently, Onuki et al. developed an elegant approach based on local bond orientational order parameter to quantify different phases, stacking faults and grain boundaries in both two and three dimensional systems22 (); 23 ().

Contrary to the common believe and classical nucleation theories that the nuclei of new phase are circular in 2D) and have same order and density as of bulk stable phase, interestingly, we find highly anisotropic growth of new phase. At the early stage of growth, the cluster grows as linear strip (having significantly different order than the bulk oblique phase) followed by branched and ring-like strip and finally engulfs the whole system. The geometry of nucleating and growing cluster (shown in Fig. 5E) is a consequence of delicate balance between the two (isotropic and anisotropic) types of interactions which enables the dominance of stabilizing energy over the destabilizing surface energy. Sharing of the H-bonds (see Fig. 5E) between the growing cluster and the bulk metastable honeycomb phase minimizes the surface energy cost for formation of new phase. We also observe that the structure of the growing cluster is different from the bulk oblique solid phase. As shown in Fig. 5B and 5C, the oblique-like phase starts appearing inside the high density intermediate order cluster. There is no clear dividing surface between the oblique and the honeycomb and hence it is difficult to assign unambiguously the identity of particles as honeycomb or oblique-like. Wetting of the oblique-like phase by intermediate order particles reduces the surface free energy cost for the formation of stable oblique phase.
Thus, the solid-solid transition in a molecular system is a two step process. In the first step small anisotropic cluster having significantly different order nucleates and grows. Anisotropic shape enables the cluster to maximize the energy gain due to strong specific interactions among particles inside the cluster as well as with the bulk parent honeycomb phase. In the second step, due to internal rearrangement inside the dense (intermediate ordered) anisotropic cluster oblique-like phases start appearing. Wetting of nuclei by intermediate order phase is also observed in crystal nucleation from metastable liquid24 () where nucleation preferentially takes place in regions of high structural order. Note that clear dividing surface with no sharing of particles between the two phases (cluster of stable phase and bulk metastable phase) creates a density/order deficiency and thus increases the surface energy25 (). Wetting is one of the factors responsible for the observed discrepancy between the predictions of classical nucleation theory and experimental/simulation studies.
In Fig. 6A, we have shown the evolution of the average order and average coordination number of a growing cluster. We consider a particle not as honeycomb like if it has ${\psi }_{3i}$ value less than 0.75. Particles having ${\psi }_{3i}$ value less than 0.75 and are connected by neighborhood form a cluster of new phase. The order of the cluster (${\psi }_{3lc}$) is defined as, ${\psi }_{3lc}={\sum }_{i=1}^{n_{lc}}{\psi }_{3i}/{\sum }_{i=1}^{n_{lc}}{\psi }_{3i}{n}_{lc}{n}_{lc}$, where $n_{lc}$ is the size of the largest cluster and summation is over all the particles that belong to the largest cluster. Average coordination number is scaled with the average coordination number of particles in the oblique solid phase. Since only one cluster grows in the system the largest cluster contains the maximum fraction of particles belonging to the new phase. The cluster of new phase has significantly different order from the bulk stable phase. During the growth, both the order and average coordination number of the cluster gradually evolve and finally merge to the value corresponding the bulk oblique solid phase. As shown in the figure, there is a crossover (indicated by blue dotted line in Fig. 6A) from weakly dependence to non-linear dependence of order as well as average coordination number on size of the cluster. At early stage of growth, weak dependence of the order and average coordination number on size of the cluster indicates the linear strip-like growth, where surface to volume ratio does not change appreciably. The non-linear dependence of order/average coordination number on cluster size is a consequence of two factors: (i) the change of surface to bulk ratio due to branching and swelling of the growing cluster and (ii) the rearrangement of the particles inside the cluster towards oblique-like structure.
In Fig. 6B, we have shown the evolution of both global six-fold bond orientational order (${\psi }_{6g}$) and local six-fold bond orientational order (${\psi }_{6l}$) of system with density during collapse transition. We observe an almost monotonic decrease of ${\psi }_{6l}$ during transition from the value in the honeycomb to its value in the oblique solid phase. However, ${\psi }_{6g}$ shows strongly non-monotonic evolution from its value in the honeycomb to the oblique solid phase. This non-monotonic evolution of ${\psi }_{6g}$ and monotonic (or very weakly non-monotonic) evolution of ${\psi }_{6l}$ arises due to cancellation of components of global bond orientational order at intermediate stage of growth and suggests the presence of domains (separated by grain boundaries) having different symmetry axes. Presence of domains having distinct symmetry axes does not alter the values of local order parameters. This scenario is also evident from the snapshots shown in Fig. 5. In the late stage of growth due to relaxation of grain boundaries the global order parameter value again increases and finally merges to the value corresponding to the bulk oblique phase.

Temperature plays two roles in the phase transition – (i) it couples with entropy and perturbs the entire free energy landscape and (ii) it provides sufficient thermal energy to overcome the free energy barrier for phase transition. In the case of condensation and crystallization, where the entropy difference between the two phases is large, the dependence of free energy barrier on temperature is stronger than the case of solid-solid phase transitions. As discussed earlier, in the case of solid-solid transition where specific interactions play important role in the stability of both the phases we can assume that the major role of temperature is to provide sufficient thermal energy to overcome the free energy barrier for transition. As the honeycomb solid is metastable at a given initial pressure and temperature with sufficiently high free energy barrier, transition can take place in two ways – (i) increase the temperature until the thermal energy is sufficient to overcome the free energy barrier (this is called temperature induced transition and discussed earlier) and (ii) increase the pressure until the given initial thermal energy is sufficient to overcome the barrier. Increase of pressure favors high density phase and thus decreases the free energy barrier for transition. In this section we shall compare these two scenarios and discuss how the pathway of pressure induced transition at low temperature differ from temperature induced transition.

In Fig. 7 we have compared the pathways for pressure and temperature induced transitions in two dimensional order parameter (density and six-fold bond orientational order) space. For temperature induced transition the transition temperature is $T$ = 0.105 at $P$ = 0.10. In the case of pressure induced transition we have increased the pressure of the system (details are mentioned in Section II) at constant temperature $T$ = 0.05 and the observed transition pressure is $P\sim 0.45$. On increasing pressure transition temperature shifts towards lower temperature. In the case of pressure induced transition, we find that different trajectories collapse to different metastable disordered states (with slightly different density and order) and transition pressure also vary in orders of 0.01. As shown in the figure, contrary to the temperature induced transition where the honeycomb solid collapses to high density crystalline oblique solid phase, in the case of pressure induced transition the honeycomb solid collapses to a long lived disordered metastable state having very small global order parameter value (see Fig. 7A). Disordered metastable phase evolves towards the stable state with an extremely slow rate. In spite of the absence of the global six-fold bond orientational order, significantly large value of the local six-fold bond orientational order (see Fig. 7B) suggests the heterogeneous structure of the metastable phase where significant local six-fold symmetry is present but there is no global translation of the local order.

In order to gain more insight in the microscopic structure of the disordered solid phase, in Fig. 8, we have shown the distributions of local six-fold bond orientational order (${\psi }_{6i}$) for both the oblique and the disordered solid phase. The emergence of the shoulder in the distribution corresponding to the disordered solid phase indicates the presence of heterogeneous local environments. To find the structures of different local environments responsible for the broadening of ${\psi }_{6i}$ distribution towards the lower value, in inset we have shown the distributions of ${\psi }_{mi}$ $\left(P\right({\psi }_{mi}\left)\right)$ for the particles having coordination number m. That is we have computed the distribution of ${\psi }_{6i}$ for only those particles which have coordination number six and so on. The distribution is normalized with the total number of particles so that area under curve of $P({\psi }_{mi}$) will give the fraction of particles having coordination number m. We note the presence of significant five-cordinated ($\sim 30\%$) and seven-coordinated ($\sim 18\%$) along with six-coordinated ($\sim 50\%$) particles with respective symmetries. On heating the metastable disordered solid phase undergoes a transition to the stable oblique phase having unimodal ${\psi }_{6i}$ distribution. On increasing pressure the peak of the distribution of ${\psi }_{6i}$ in the oblique solid phase shifts towards lower ${\psi }_{6i}$ value.
Recently, Kadau et al.5 () used molecular dynamics simulations to investigate the shock-induced phase transformation of solid iron. Above a critical shock strength the nucleation and growth become collective as many small grains start growing simultaneously. However, we have observed a single cluster growth scenario in both pressure and temperature induced transitions. Collective nucleation and growth scenario indicates the absence of free energy barrier for nucleation. Note that shock waves increase both the pressure as well as temperature and increment in both decreases the free energy barrier for transition. Above a critical shock strength the free energy barrier might be negligibly small and hence many clusters start growing simultaneously. This transition closely resembles with the gas-liquid transition at large metastabilty (above kinetic spinodal) where one observes similar collective nature of nucleation and growth26 (); 27 (); 28 (). In the case of pressure induced transition (or low temperature transition as on increasing pressure transition point shifts towards the lower temperature) trapping of the system in a long lived local minima indicates that the free energy landscape is rugged. Although one might observe the rugged energy landscape in atomic systems also, inclusion of specificity in interaction potential enhances the ruggedness of the landscape. The origin of enhanced ruggedness lies in the competition between specific and non-specific interactions where local regions might be energetically stabilized by a particular type of interaction but system as a whole might not be in global minimum. These local traps in the landscape decreases significantly the dynamics of the transition. On increasing temperature the disordered metastable solid phase undergoes a transition to the crystalline oblique phase as on increasing temperature system attains sufficient thermal energy to overcome local energy barriers and finally undergoes a transition to global free energy minimum.