# Defect-unbinding transitions and inherent structures in two dimensions

## I Introduction

Some years ago, Stillinger and Weber  introduced a theory of liquids, based on the partitioning of the configuration space into potential-energy (PE) basins. Each of these basins contains a single PE minimum, to which all other points within the basin are connected via steepest-descent paths. The PE minima were coined “inherent structures” (IS); all other configurations are taken to be vibrational excitations of them. This approach allows for the decomposition of the configurational partition function into a sum, over PE basins, of intrabasin terms. The resulting partition function may be approximated as follows:

$$Q=\underset{\alpha }{}Q_\alpha G(𝐩)Q_𝐩𝑑𝐩G(𝐩^{})Q_𝐩^{}.$$
(1)   
In the first step of this transformation, one splits the partition function into a sum over “basin partition functions” (the usual Boltzmann integral, limited to configurations within a given basin), $`\alpha `$ being the basin index. This step is exact in principle. The second step transforms this sum into an integral, via the introduction of the (generally vector-valued) structural parameter $`𝐩`$, characterizing the IS. Typically, this parameter would include such information as average coordination numbers, densities and spatial distributions of defects, etc. $`G(𝐩)`$ is a density-of-states function, enumerating the basins having a given value of $`𝐩`$; $`Q_𝐩`$ is then the corresponding constrained partition function. This step in the transformation inevitably loses some information through the necessity to make $`𝐩`$ finite-dimensional. The final step of the transformation is a result of the fact that, in the thermodynamic limit, the density-of-states function is essentially exponential in the number of particles:

$$G(N,𝐩)\mathrm{exp}[Ng(𝐩)],$$
(2)   
which permits a maximum-integrand evaluation of the integral over $`𝐩`$. $`𝐩^{}`$ is the value of $`𝐩`$ which maximizes the integrand for a given set of thermodynamic conditions (e.g. volume and temperature).

The partition function may be further transformed by writing the potential energy as

$$\mathrm{\Phi }(𝐫)=\mathrm{\Phi }_\alpha +\mathrm{\Delta }_\alpha \mathrm{\Phi }(𝐫),$$
(3)   
where $`\mathrm{\Phi }_\alpha `$ is the “structural energy” (the PE of the IS of the occupied basin) and $`\mathrm{\Delta }_\alpha \mathrm{\Phi }(𝐫)`$ is the “vibrational energy” (the difference between the total PE and the structural energy). This allows the generic basin partition function to be written as

$$Q_\alpha =\mathrm{exp}(\mathrm{\Phi }_\alpha /k_BT)Q_\alpha ^{vib},$$
(4)   
where

$$Q_\alpha ^{vib}=_{R(\alpha )}\mathrm{exp}(\mathrm{\Delta }_\alpha \mathrm{\Phi }/k_BT)𝑑𝐫.$$
(5)   
Here, $`R(\alpha )`$ limits the integration to basin $`\alpha `$, and $`k_B`$ is Boltzmann’s constant. Hence, the total partition function, as given by Eq. (1.1), becomes

$$Q=Q_𝐩^{}^{struct}Q_𝐩^{}^{vib},$$
(6)   
where,

$$Q_𝐩^{}^{struct}=G(𝐩^{})\mathrm{exp}(\mathrm{\Phi }_𝐩^{}/k_BT).$$
(7)   
Phase transitions are defined by singularities in the free energy, $`F_{conf}=k_BT\mathrm{ln}Q`$. Thus, in order for a phase transition to occur, there must be singularities in $`Q_𝐩^{}^{struct}`$, $`Q_𝐩^{}^{vib}`$, or both. The evident unlikelihood of such singularities, without discontinuities in $`𝐩^{}`$ (that is to say in the types of basins occupied), implies that the existence of marked differences between the IS associated with different phases is a practical requirement for the applicability of the inherent-structures theory (IST) to systems exhibiting phase transitions. This requirement has been shown to be satisfied for certain three-dimensional systems . Also, some limited results  have been obtained for 2D systems.

Recently, we  have performed extensive numerical studies of inherent structures in simple, single-component fluids in 2D. An early study  extended the range of such studies to $`N=4096`$ particles, but found (as expected) only evidence for two phases (solid and liquid). More recently  we extended these studies to $`N=36864`$ particles. This choice of system size was motivated by earlier molecular-dynamics studies  giving strong—but not conclusive—evidence for an intermediate, hexatic phase for systems of this size and larger (see also Ref. ). Ref.  gave a brief report of the principal results reported here. In this paper, we offer a detailed discussion of our computational methods and results. We also provide a clear picture of the disclination-unbinding “transition” in the IS—a result which was not clear in —and some calculations of the disclination charge-charge correlation function in the equilibrium fluids. These latter calculations provide further evidence for disclination unbinding at the hexatic/liquid transition.

Roughly contemporaneously with the development of IST, Halperin and Nelson (HN) , following work on the melting of 2D solids by Kosterlitz and Thouless , predicted a two-stage melting mechanism for 2D systems. A number of results on the first stage of melting were obtained independently by Young  (see also Nelson ). In the resulting picture of two-stage melting—commonly called the ‘KTHNY theory’—each successive phase (in order of increasing energy) is characterized by the presence of an additional type of defect: the solid contains only dislocations, bound in pairs; the intermediate, hexatic phase adds unbound dislocations; and the liquid further adds unbound disclinations. Attendant to this progression of defects are differences in the bond-orientational correlation function, which exhibits long-range order, power-law decay, and exponential decay for the solid, hexatic, and liquid phases, respectively. (For a detailed review of the KTHNY theory and the defects involved, see the review of Strandburg .) Since KTHNY predicts the existence of three condensed phases—for which IST requires an equal number of distinct classes of IS—2D would seem to offer an ideal testing ground for IST. Furthermore, it seems that a study of the IS underlying the different phases in 2D systems might provide novel and useful microscopic evidence for the defect-unbinding transitions expected from the theory—assuming that such defects can be ‘trapped’ by the quenching procedure (which yields the mechanically stable inherent structure from a snapshot configuration at thermal equilibrium). Indeed, the defects present in each equilibrium phase should also show up in the IS underlying that phase—but much more clearly, due to the attendant removal of the vibrational distortions present at equilibrium—if these defects are mechanically stable. Before discussing our numerical results, then, we will discuss the question of the mechanical stability of the defects—dislocations and disclinations.

The interaction energies of the defects are typically calculated by well known methods of linear elasticity theory . To make an elastic defect, one takes an unstrained elastic medium and introduces one or more topological changes—yielding one or more topological defects—by cutting, shifting, and gluing. The next step is to invoke the “equations of equilibrium” (EOE)—which require that there be no net force at any point in the medium—at constant topology. This then gives the strain field which, in turn, gives the self-energies and energies of interaction for the defects. The assumption of constant topology is crucial, because it is the topology that defines the elastic problem to be solved. However, it is important to note that this can only be an assumption in the continuum theory: in the absence of a microscopic atomic structure in the fluid, there is no reason to expect any dislocation to be pinned (against the calculated forces of attraction or repulsion) at any point in the medium, except ‘by hand’.

Thus we cannot expect to see such structures as free dislocations in any mechanically stable configuration, unless some justification can be given for this assumption of constant topology. For the case of dislocations, this justification comes in the form of the Peierls-Nabarro potential . This is a periodic “corrugation” in the inter-dislocation potential, arising from the underlying microscopic structure of the material. This potential is well known to be capable of pinning dislocations, such that arrays of dislocations may be rendered mechanically stable. This leads us to anticipate that the configurations of dislocations trapped in our numerically-obtained IS may in fact provide useful insight into the equilibrium defect structures, without the almost overwhelming ‘noise’ associated with the vibrations about the IS, occurring in thermal equilibrium.

There is less justification for this assumption, as applied to disclinations. In fact, there is reason to doubt the mechanical stability—and hence the presence in IS—of free disclinations . The question is then, can the inherent-structures idea, invoking as it does a qualitative difference in IS between different thermodynamic phases, be reconciled with the KTHNY picture of melting (hexatic $``$ liquid) by the unbinding of disclinations—even when there is good reason to expect that no free disclinations can be seen in mechanical equilibrium (i.e., in any IS)? We provide a conclusive (‘yes’) answer to this question, below, while at the same time failing to find any evidence for mechanical stability of free disclinations.

Our results reveal an extremely clean correspondence between the predictions of the KTHNY theory of two-stage melting and the inherent structures associated with each thermodynamic phase. Previously, the principal barrier to this sort of study has been the difficulty associated with finding the hexatic phase in simulations , in part due to limitations in system size. Boundary conditions and equilibration methods may also play an important role. The simulations in Refs.  ($`N>100000`$) and  ($`N65000`$) gave some compelling evidence for the hexatic phase. However this phase was found to be only metastable thermodynamically in ; and the differences in method and boundary conditions in the two studies leaves some room for controversy. Our own studies use a system size ($`N36000`$) for which a metastable hexatic phase appeared in the study of Ref. , with quenches from the ‘equilibrium’ hexatic phase being taken from snapshots in this metastable thermodynamic state. We believe that our results, revealing as they do three classes (one clearly ‘hexatic’) of IS for these fluids—classes which are expected to persist in the thermodynamic limit—provide further support for the hypothesis that the hexatic phase is thermodynamically stable, in some region of the phase diagram, for $`N\mathrm{}`$. This evidence is distinct from, and complementary to, that obtained from equilibrium studies. Our results also further strengthen the basic premise of inherent-structures theory: that distinct thermodynamic phases are characterized by qualitatively distinct inherent structures, such that singularities in thermodynamic functions may be ascribed to discontinuities in the occupation probabilities of potential-energy basins, at a phase boundary.

## II Computational procedure

In order to study inherent structures, one first needs starting configurations, taken as snapshots from thermal equilibrium. In the present study, these were obtained directly from the molecular-dynamics (MD) simulations described in Ref. , which made use of a computational framework described by Melchionna, Ciccotti, and Holian (MCH)—specifically, a constant-NPT molecular dynamics simulation, using the Parrinello-Rahman  shape-varying box with periodic boundary conditions (BC). Let us, then, give a brief summary of those aspects of this framework which are most relevant to the present study, the details being available in Ref. .

For the sake of simplicity, we will outline the method for the case of isotropic volume fluctuations and then state the changes necessary to account for the shape-varying box. The MCH equations of motion, for the case of isotropic volume fluctuations, are

$`\dot{𝐫}_i={\displaystyle \frac{𝐩_i}{m_i}}+\eta (𝐫_i𝐑_0)`$ (8)   
$`\dot{𝐩}_i=𝐅_i(\eta +\zeta )𝐩_i`$ (9)   
$`\dot{\zeta }=\nu _T^2\left[{\displaystyle \frac{T(t)}{T_{ext}}}1\right]`$ (10)   
$`\dot{\eta }={\displaystyle \frac{\nu _P^2}{Nk_BT_{ext}}}V[P(t)P_{ext}]`$ (11)   
$`\dot{V}=dV\eta .`$ (12)   
Here, $`𝐫_i`$, $`𝐩_i`$, and $`m_i`$ are the position, momentum, and mass, respectively, of particle $`i`$. $`𝐅_i`$ is the instantaneous force acting on particle $`i`$, and $`𝐑_0`$ is the center of mass of the system. $`\eta `$ is a barostating factor which tends to restore the instantaneous pressure, $`P(t)`$, to the preset value $`P_{ext}`$. It is modulated by the adjustable parameter $`\nu _P`$, which is termed the “barostating rate”. Similarly, $`\zeta `$ serves to equilibrate the temperature and is tuned by way of the “thermostating rate”, $`\nu _T`$. $`V`$ is the volume and $`d`$ is the dimensionality. The main alterations to Eq. (10), needed to accommodate a shape-varying box, are to change the scalars $`\eta `$ and $`P`$ to tensors and to replace $`V`$ with a “box matrix” whose columns are the basis vectors of the box. The basic form of the equations remains that of Eq. (10). The potential used was of the shifted Lennard-Jones form,

$$V(r)=\{\begin{array}{ccccccccc}4ϵ[(\sigma /r)^{12}(\sigma /r)^6]+V_c,\hfill & & & & & & & & r<r_c,\\ 0,\hfill & & & & & & & & rr_c,\end{array}$$
(13)   
where $`ϵ`$ and $`\sigma `$ are parameters, $`r`$ is the interatomic distance, $`r_c`$ is the cutoff radius, and

$$V_c=4ϵ[(\sigma /r_c)^{12}(\sigma /r_c)^6].$$
(14)   
In units in which $`ϵ`$, $`\sigma `$, $`k_B`$, and $`m`$ (the atomic mass) are all equal to 1—units which we use throughout this paper—the parameters used in the equilibrations were as follows:

$$\begin{array}{c}N=36864\hfill \\ P_{ext}=20\hfill \\ T_{ext}\{\begin{array}{ccc}2.15& & (\mathrm{𝖼𝗋𝗒𝗌𝗍𝖺𝗅})\hfill \\ =2.154& & (\mathrm{𝗁𝖾𝗑𝖺𝗍𝗂𝖼})\hfill \\ 2.17& & (\mathrm{𝗅𝗂𝗊𝗎𝗂𝖽})\hfill \end{array}\hfill \\ r_c=4\hfill \\ \mathrm{\Delta }t=.0005,\hfill \end{array}$$
(15)   
where $`N`$ is, of course, the number of particles.

Having obtained an equilibrium configuration (in the hexatic case, this is only a thermodynamically metastable equilibrium ), in order to find the associated IS one must perform a steepest-descent minimization of the total potential energy—which, for the current constant-$`P`$ case, consists of the total interparticle potential plus $`P_{ext}V`$. In practice, this is too time consuming to be practical  for the large system size of the current study—indeed, it is impractical for all but the smallest systems. Instead, we make use of the above-outlined MD method. Beginning with an equilibrium configuration, we first zero all particle velocities and the velocities of the “box walls”. We then run the MD simulation at very low temperature ($`T_{ext}`$ something like $`10^6`$ to $`10^4`$ of the equilibrium temperature), carefully adjusting the thermostating and barostating rates, such that the instantaneous temperature remains very close to $`T_{ext}`$ and the PE smoothly decreases with the time. Furthermore, the PE is checked at each time step, and if an increase is found, we go back to a previous configuration (saved before the occurrence of the increase), re-zero the velocities, and restart the simulation. This process is continued until the duration of the MD runs (i.e. before a PE increase occurs) becomes only a few time steps. At this point, we run the MD simulation at the same, very low temperature (and without, of course, the requirement of a strictly decreasing PE), long enough (typically something like $`10^5`$ time steps) to ensure that the system is, indeed, vibrating about a PE minimum. If this test is successful, we have an IS. (At such low temperatures, the vibrations are of small enough amplitude to be negligible for structural considerations.) Otherwise, we continue the minimization procedure, until a minimum is found which does pass our test.

Of course, the structures obtained from the above-outlined minimization procedure (hereafter referred to as “quenching”) will, in general, differ from the ‘true’ IS, connected by a steepest-descent path to the starting equilibrium configuration. However, the very low temperatures—which is to say very small particle momenta—maintained throughout the quenching process ensure that the system trajectory nearly  follows the steepest-descent path prescribed by IST. Also, according to IST, fluid systems at equilibrium transit frequently and exclusively (at least for $`N\mathrm{}`$) between thermodynamically equivalent PE basins  (i.e. basins having the same value of the structural parameter, $`𝐩`$). Hence we believe that the structures obtained by our quenching procedure will be thermodynamically equivalent to the actual IS of the starting configuration. That is, they should be representative of the set of IS associated with the thermodynamic conditions of the equilibrium configuration.

## III Results

In earlier work  we performed a large number of quenches of equilibrium systems at $`N=4096`$. In the present work we increase $`N`$ to 36864 particles, in order to be able to quench from all three thermal phases: solid, hexatic, and liquid. At this system size the quenches are very intensive computationally. We have performed three quenches—following the procedure outlined above—for each phase, and several additional quenches which did not strictly enforce the requirement of monotonic decrease of the PE. The results were qualitatively the same for all quenches derived from the same starting phase.

Fig. 1(a) shows an equilibrium snapshot for the solid phase. We plot only those atoms which are not 6-fold coordinated according to a Voronoi construction, labeling all such atoms by their coordination number (mostly 5 or 7). Fig. 1(b) then shows the configuration of Fig. 1(a), when relaxed to mechanical equilibrium by our quenching procedure. Although Fig. 1(b) is mostly white space (i.e., 6-fold coordinated atoms), we include it here to illustrate the dramatic reduction in defect number as a result of quenching from thermal equilibrium to mechanical equilibrium. It is apparent to the eye that there are no free dislocations in the solid-derived IS: every dislocation is closely bound in a (vector-)charge-neutral composite. We have also computed the bond-orientational correlation functions (BOCF) for the various quenches. Obviously, for the solid-derived IS, the BOCF has long-ranged order.

On quenching from the (metastable) equilibrium hexatic phase, we obtain structures such as that shown in Figure 2. (We do not show the equilibrium defect configuration as the defects are very dense.) Clearly, there is a large density of defects, even in mechanical equilibrium, for this case. Of course, there are still some bound dislocations, some of which compose large-angle grain boundaries identifiable as chains of very closely spaced dislocations. In addition to these—and in contrast to the crystal IS—there are many dislocations which do not have any ‘canceling’ dislocations within several lattice spacings, some of which show a clear tendency  to arrange themselves into small-angle grain boundaries. We term these the “free” dislocations for the IS; reasoning from the existence of the Peierls-Nabarro potential, we argued above that free dislocations, if present in the equilibrium snapshot, will survive the quench and thus appear in the IS. Here, we claim that a comparison of Figs. 1 and 2 graphically reveals the dislocation-unbinding transition in the inherent structures.

We can also test this idea with the BOCF. While networks of large-angle grain boundaries are capable of destroying the quasi-long-range orientational order characteristic of the hexatic phase, those present in our hexatic quenches are relatively small and isolated, so that this order is in fact preserved. Log-log plots of the BOCF \[$`g_6(r)`$\] for a typical hexatic “MD snapshot”  and its associated quenched structure  reveal that both obey a power-law behavior, with the IS showing a smaller exponent (i.e. a slower rate of algebraic decay of the orientational order). This may be attributed to the removal, on quenching, of long-wavelength torsional phonons, which are supported by the hexatic phase’s finite Frank constant.

An IS of this nature has not been seen in any smaller system. In fact—if we use the word “glass” as shorthand for structural glass, i.e. an atomic configuration in mechanical (but not thermal) equilibrium—then Fig. 2 shows a hexatic glass. Two-dimensional glasses have mostly been studied using two or more atomic species , since the ‘frustration’ in 2D monatomic fluids is very small (it is zero for the 2D packing problem) . Hexatic glasses have, to our knowledge, only been seen before in two-component systems—in simulations , and in ball-bearing experiments .

With the isotropic liquid as the starting point, our quenching procedure results in structures such as that shown in Figure 3. The two-stage melting theory predicts that the transition from hexatic to isotropic liquid takes place via the unbinding of disclinations (our 5’s, 7’s etc.). However there are no free disclinations in any of our liquid-derived IS. Rather, the only additional defects, as compared to the hexatic phase, are—as is clear from Fig. 3—percolating networks of large-angle grain boundaries. These in themselves can destroy the quasi-long-ranged orientational ordering, as may be verified by calculating the BOCF for the IS of Fig. 3. In fact, both the equilibrium snapshot and its IS show an exponential decay of orientational order, with roughly the same exponent .

It is thus tempting to suppose that the equilibrium liquid is also characterized by percolation of grain boundaries—i.e., that the hexatic $``$ liquid transition takes place by grain-boundary melting . Certainly \[compare Figs. 2(a) and 3\] the transition appears in the IS as a percolation of grain boundaries. However, we believe that we can see the unbinding of disclinations in our IS, with a bit more effort.

In Fig. 4 we show the quenching of an artificial starting condition whose only defects are four widely spaced disclinations (two positive and two negative). The corresponding quenched (mechanically stable) structure is a roughly square grain-boundary network, whose nodes correspond closely to the positions of the original disclinations. The “free” disclinations of Fig. 4(a) are (as expected) not mechanically stable; and they relax upon quenching to a network of grain boundaries, which serves as a “fossil relic” of the free disclinations in the starting configuration. This suggests that there should be a strong correlation between the average separation of free disclinations in an equilibrium configuration, and the average grain size in the corresponding quenched structure. We find further support for this idea from other quenches like that shown in Fig. 4: above a threshold separation distance, the grain size in the quenched structure closely reflects the spacing of the original disclinations. (For disclinations closer than the threshold distance, the relaxed structure is a single grain.) Thus, the fact that the equilibrium and quenched configurations have nearly identical orientational correlation lengths  is consistent with the presence of free disclinations in the equilibrium liquid.

It is also interesting to note that, while the transformation shown in Fig. 4 is quite dramatic, the final configuration \[Fig. 4(b)\] still, in a sense, contains “free” disclinations. These are at the nodes of the grain-boundary network, which are near the positions of the original disclinations of Fig. 4(a). Disclinations are defined as centers of lattice rotation; that is, tracking the local lattice orientation, while making a closed circuit around a disclination, will show a net rotation. For the present case, this means that if we track the orientation of six-coordinated cells, as we make a closed circuit around a disclination, we will find a net rotation of some integral multiple of $`\pi /3`$. The disclinations centered on the five- and seven-coordinated atoms of Fig. 4(a) give rotations of $`+\pi /3`$ and $`\pi /3`$, respectively. Similarly, on making a circuit enclosing a set ($`A`$) of atoms whose average coordination number is different from six, we will find a net lattice rotation—specifically,

$`\theta _{rot.}={\displaystyle \frac{\pi }{3}}q_A`$ (16)   
$`q_A={\displaystyle \underset{iA}{}}(6z_i).`$ (17)   
Here, $`q_A`$ is the net “disclination charge” in $`A`$, and $`z_i`$ is the coordination number of atom $`i`$. Of course, in order to properly define this lattice rotation, we need a circuit consisting solely of “good crystal” (i.e. six-coordinated cells). Such circuits do exist around the grain-boundary nodes of Fig. 4(b) (as shown in closeup in Figure 5). Hence we see that there are, indeed, “net-sevens” ($`q_A=1`$) and “net-fives” ($`q_A=1`$) at the grain boundary nodes, near the positions of the original negative and positive disclinations, respectively.

We next examine the distribution of such “net disclinicity” in our IS. In doing this, the disclinations are identified as outlined above—by identifying groups of atoms whose average coordination number is different from six, and which can be enclosed by a path consisting entirely of six-coordinated atoms. If, for a given disclination, the smallest such path encloses more than one atom, the location of the disclination is somewhat arbitrary (except, of course, that it should be somewhere within the enclosed area). For the purpose of illustration in the present work, we have used the following rule: if the disclination is near a grain-boundary node, we assign it to the appropriately coordinated atom ($`z=5`$ for a positive disclination, or $`z=7`$ for a negative disclination) nearest the node; otherwise, we assign it to the appropriately coordinated atom nearest to the closest, opposite-signed disclination. Application of this procedure to the quenched structures of Figs. 2 and 3 results in the disclination arrays shown in Figure 6. Comparison of these arrays with the structures from which they were obtained reveals that the presence of net disclinations is restricted to the vicinity of large-angle grain boundaries. That is, Figs. 2, 3, and 6, considered together, are entirely consistent with the notion obtained from Fig. 4: that free disclinations, upon quenching, relax to mechanically stable IS in which a network of grain boundaries marks the extent, and even positions, of the original distribution of free disclinations. In this view, then, the grain-boundary-percolation transition seen in our IS is the direct analog of the disclination-unbinding transition in equilibrium. In other words, if the disclinations unbind, the grain boundaries in the IS must percolate; and, on the other hand, bound disclinicity appears in the IS as localized (or no) grain boundaries. Our net-disclination algorithm then simply erases the grain-boundary network, revealing bound disclinations in our hexatic IS \[Fig. 6(a)\] and unbound disclinations in the liquid IS \[Fig. 6(b)\].

Thus we believe that we have seen both defect-unbinding transitions reflected in our IS. However this conclusion requires a chain of reasoning which is not airtight. Hence, as a further test of the hypothesis that disclination unbinding mediates the hexatic $``$ liquid transition, we have analyzed the disclination distribution of the two equilibrium phases, using the disclination “charge-charge” correlation function (CCF):

$`g_q(r)={\displaystyle \underset{ij}{}}\delta (rr_{ij})q_iq_j`$ (18)   
$`r_{ij}=|𝐫_i𝐫_j|`$ (19)   
$`q_i=6z_i,`$ (20)   
where the sum is over all pairs of atoms. According to Halperin , the absolute value of this function should exhibit (asymptotically) a power-law decay when the disclinations are bound (but interacting with a logarithmic potential, as in the hexatic phase), and an exponential decay for the case of unbound disclinations. As suggested by Fig. 6, typical equilibrium configurations contain relatively few candidates for free disclinations. Additionally, at equilibrium there are very many non-sixfold-coordinated atoms that do not comprise unbound disclinations—they are components of dislocations, both stable and “virtual”—so that the disclination CCF must be extracted from the considerable noise resulting from these other defects. For these reasons, the elucidation of the asymptotic behavior of the CCF requires long-time averaging over very many equilibrium configurations. This is an extremely time consuming process. We have obtained results which, we believe, offer some evidence for the unbinding of disclinations in the liquid phase. However, even after much time averaging, our CCFs show meaningful behavior only over a single decade of distance; hence we cannot claim that this evidence itself is conclusive.

Our CCF results for 36K particles are presented in Figure 7 (hexatic phase) and Figure 8 (liquid). The severe noise in these data is clearly evident. However we also believe that there is clear visual support in these data for the hypothesis that the hexatic CCF is better fit by a power-law or subexponential form, while the liquid CCF is better fit by an exponential. Obviously one can only draw very tentative conclusions from these data. However, from our experience we believe that further time averaging will not significantly reduce the noise in these data; instead, larger systems (with larger numbers of defects) need to be examined.

## IV Discussion

Our principal conclusion from these studies is that the basic premise of inherent-structures theory is well confirmed for 2D simple fluids. That is, inherent structures obtained from differing equilibrium phases differ from one another in qualitative and reproducible ways. Our own results extend those from previous IS studies of 2D fluids  in two ways. We have studied system sizes 1–2 orders of magnitude beyond those of previous studies. This increase in size has enabled us to obtain, for the first time, convincing differences among three distinct phases in 2D.

In particular, we find unambiguous evidence for the existence of hexatic inherent structures, which are also hexatic (structural) glasses. Such structures do not appear at the level $`N=4000`$ , but are found at $`N=36000`$—at which size the equilibrium hexatic phase is thermodynamically metastable . It is certain that the hexatic IS will persist for all larger $`N`$. The question is then (in the language of IST), can these IS reach a balance of entropy and energy (with increasing $`N`$) such that the hexatic phase becomes a true minimum of the free energy, in some part of the phase diagram? Beyond the minimal criterion of showing the existence of such IS, our present results do not answer this question. However, more detailed studies of the IS reported here may allow further progress.

Our results are also relevant to the KTHNY theory of two-stage melting. Besides the fact that we have shown three “phases” of inherent structures—with the same types of bond-orientational order as the corresponding equilibrium phases—we have also obtained graphic evidence for both defect-unbinding transitions. The unbinding of dislocations, as one moves from solid-derived to hexatic-derived IS, is apparent. In contrast, free disclinations are (apparently) not mechanically stable \[Fig. 4\], and so hidden in the liquid-derived IS. However, we have used a simple, deterministic algorithm which identifies and reveals net disclinicity (i.e., that not canceled by neighboring atoms) in any IS. This rule, applied to our hexatic- and liquid-derived IS, yields configurations of defects showing a clear disclination-unbinding transition. This latter transition shows up as a percolation transition for grain boundaries in the untransformed IS; the net disclinations are then (typically) found at the nodes of the grain-boundary network. We have also obtained some further evidence for disclination unbinding by studying the disclination charge-charge correlation function for the equilibrium fluid.

Finally, we reiterate that inherent structures are, in principle, structural glasses: nonequilibrium (and disordered) configurations trapped away from the equilibrium state by potential-energy barriers. Hence we find that there are two distinct types of structural glasses for one-component, Lennard-Jones fluids in 2D. These two types might be experimentally observable—although it seems likely that the required quench rates are too high. The idea however has considerable interest. Just as the novel phases found in 2D equilibrium phase diagrams have enhanced our knowledge of the phases of fluid matter in general, so may unusual structure in a 2D “glassy phase diagram” be expected to offer new insights into our general understanding of glassy matter.

Acknowledgments.– We thank David Nelson for helpful comments regarding free disclinations, and Ken Stephenson of the University of Tennessee’s Mathematics Department for generating starting structures such as that shown in Fig. 4(a). This work was supported in part by the Applied Mathematical Science Research Program, Office of Energy Research, and the Division of Materials Science, U.S. Department of Energy, through Contract No. DE-AC05-84OR21400 with Martin Marietta Energy Systems Inc. FLS and GSC were supported by the NSF under Grant # DMR-9413057.