KTHNY theory

In statistical mechanics, the KTHNY-theory describes the process of melting of crystals in two dimensions (2D). The name is derived from the initials of the surnames of John Michael Kosterlitz, David J. Thouless, Bertrand Halperin, David R. Nelson, and A. Peter Young, who developed the theory in the 1970s. It is, beside the Ising model in 2D and the XY model in 2D, one of the few theories, which can be solved analytically and which predicts a phase transition at a temperature $T>0$.

Main idea

Melting of 2D crystals is mediated by the dissociation of topological defects, which destroy the order of the crystal. In 2016, Michael Kosterlitz and David Thouless were awarded with the Nobel prize in physics for their idea, how thermally excited pairs of virtual dislocations induce a softening (described by renormalization group theory) of the crystal during heating. The shear elasticity disappears simultaneously with the dissociation of the dislocations, indicating a fluid phase. Based on this work, David Nelson and Bertrand Halperin showed, that the resulting hexatic phase is not yet an isotropic fluid. Starting from a hexagonal crystal (which is the densest packed structure in 2D), the hexatic phase has a six-folded director field, similar to liquid crystals. Orientational order only disappears due to the dissociations of a second class of topological defects, named disclinations. Peter Young calculated the critical exponent of the diverging correlations length at the transition between crystalline and hexatic. KTHNY theory predicts two continuous phase transitions, thus latent heat and phase coexistence is ruled out. The thermodynamic phases can be distinguished based on discrete versus continuous translational and orientational order. One of the transitions separates a solid phase with quasi-long range translational order and perfect long ranged orientational order from the hexatic phase. The hexatic phase shows short ranged translational order and quasi-long ranged orientational order. The second phase transition separates the hexatic phase from the isotropic fluid, where both, translational and orientational order is short ranged. The system is dominated by critical fluctuations, since for continuous transitions, the difference of energy between the thermodynamic phases disappears in the vicinity of the transition. This implies, that ordered and disordered regions fluctuate strongly in space and time. The size of those regions grows strongly near the transitions and diverges at the transition itself. At this point, the pattern of symmetry broken versus symmetric domains is fractal. Fractals are characterized by a scaling invariance – they appear similar on an arbitrary scale or by arbitrarily zooming in (this is true on any scale larger than the atomic distance). The scale invariance is the basis to use the renormalization group theory to describe the phase transitions. Both transitions are accompanied by spontaneous symmetry breaking. Unlike for melting in three dimensions, translational and orientational symmetry breaking does not need to appear simultaneously in 2D, since two different types of topological defects destroy the different types of order.

Background

Michael Kosterlitz and David Thouless tried to resolve a contradiction about 2D crystals: on one hand side, the Mermin-Wagner theorem claims that symmetry breaking of a continuous order-parameter cannot exist in two dimensions. This implies, that perfect long range positional order is ruled out in 2D crystals. On the other side, very early computer simulations of Berni Alder and Thomas E. Wainwright indicated crystallization in 2D. The KTHNY theory shows implicitly that periodicity is not a necessary criterion for a solid (this is already indicated by the existence of amorphous solids like glasses). Following M. Kosterlitz, a finite shear elasticity defines a 2D solid, including quasicrystals in this description.

Structure factor in 2D

Figure 1: Structure factor of a) an isotropic fluid, b) the hexatic phase, c) a crystal in two dimensions.

All three thermodynamic phases and their corresponding symmetries can be visualized using the structure factor:$S({\vec {q}})={\frac {1}{N}}\langle \sum _{ij}e^{-i{\vec {q}}({\vec {r}}_{i}-{\vec {r}}_{j})}\rangle$. The double sum runs over all positions of particle pairs i and j and the brackets denote an average about various configurations. The isotropic phase is characterized by concentric rings at $q=2\pi /a$, if $a=1/{\sqrt {\rho }}$ is the average particle distance calculated by the 2D particle density $\rho$. The (closed packed) crystalline phase is characterized by six-fold symmetry based on the orientational order. Unlike in 3D, where the peaks are arbitrarily sharp ($\delta$-peaks), the 2D peaks have a finite width described with a Lorenz-curve. This is due to the fact, that the translational order is only quasi-long ranged as predicted by the Mermin-Wagner theorem. The hexatic phase is characterized by six segments, which reflect the quasi-long ranged orientational order. The structure factor of Figure 1 is calculated from the positions of a colloidal monolayer (crosses at high intensity are artefacts from the Fourier transformation due to the finite (rectangular) field of view of the ensemble).

Interaction between dislocations

Figure 2: If Youngs modulus becomes $16\pi$, elasticity disappears discontinuously and the crystal melts.

To analyse melting due to the dissociation of dislocations, one starts with the energy $H_{loc}$ as function of distance between two dislocations. An isolated dislocation in 2D is a local distortions of the six-folded lattice, where neighbouring particles have five- and seven nearest neighbours, instead of six. It is important to note, that dislocations can only be created in pairs, due to topological reasons. A bound pair of dislocations is a local configuration with 5-7-7-5 neighbourhood.

$H_{loc}=-{\frac {a^{2}Y}{8\pi }}\sum _{k\neq l}{\Big [}{\vec {b}}({\vec {r}}_{k})\cdot {\vec {b}}({\vec {r}}_{l})\ln {\frac {\Delta {\vec {r}}_{k,l}}{a}}-{\frac {[{\vec {b}}({\vec {r}}_{k})\cdot \Delta {\vec {r}}_{k,l}][{\vec {b}}({\vec {r}}_{l})\cdot \Delta {\vec {r}}_{k,l}]}{\Delta r_{i,j}^{2}}}{\Big ]}+E_{c}\cdot N_{loc}.$

The double sum runs over all positions of defect pairs $k$ and $l$, $\Delta {\vec {r}}_{k,l}={\vec {r}}_{k}-{\vec {r}}_{l}$ measures the distance between the dislocations. ${\vec {b}}$ is the Burgers vector and denotes the orientation of the dislocation at position Orte ${\vec {r}}_{k}$. The second term in the brackets brings dislocations to arrange preferentially antiparallel due to energetic reasons. Its contribution is small and can be neglected for large distance between defects. The main contribution stems from the logarithmic term (the first one in the brackets) which describes, how the energy of a dislocation pair diverges with increasing distance. Since the shortest distance between two dislocations is given approximatively by the average particle distance $a$, the scaling of distances with $a$ prevents the logarithm $\ln {\frac {\Delta {\vec {r}}_{k,l}}{a}}$ to become negative. The strength of the interaction is proportional to Young's modulus $Y$ given by the stiffness of the crystal lattice. To create a dislocation from an undisturbed lattice, a small displacement on a scale smaller than the average particle distance $a$ is needed. The discrete energy associated with this displacement is usually called core energy Energie $E_{c}$ and has to be counted for each of the $N_{loc}$ dislocations individually (last term).

An easy argument for the dominating logarithmic term is, that the magnitude of the strain induced by an isolated dislocation decays according to $\propto {\frac {1}{r}}$ with distance. Assuming Hooke's approximation, the associated stress is linear with the strain. Integrating the strain ~1/r gives the energy proportional to the logarithm. The logarithmic distance dependence of the energy is the reason, why KTHNY-theory is one of the few theories of phase transitions which can be solved analytically: in statistical physics one has to calculate partition functions, e.g. the probability distribution for all possible configurations of dislocation pairs given by the Boltzmann distribution $e^{\frac {H_{loc}}{k_{B}T}}$. Here, $k_{B}T$ is the thermal energy with Boltzmann constant $k_{B}$. For the majority of problems in statistical physics one can hardly solve the partition function due to the enormous amount of particles and degrees of freedoms. This is different in KTHNY theory due to the logarithmic energy functions of dislocations $H_{loc}$ and the e-function from the Boltzmann factor as inverse which can be solved easily.

Example

We want to calculate the mean squared distance between two dislocations considering only the dominant logarithmic term for simplicity:

$\langle r^{2}\rangle ={\frac {\int r^{2}\cdot e^{-{\frac {Ya\ln(r/a)}{4\pi k_{B}T}}}d^{2}r}{\int e^{-{\frac {Ya\ln(r/a)}{4\pi k_{B}T}}}d^{2}r}}\sim {\frac {2-{\frac {Y\cdot a}{4\pi k_{B}T}}}{4-{\frac {Y\cdot a}{4\pi k_{B}T}}}}.$

This mean distance $\langle r^{2}\rangle \to 0$ tends to zero for low temperatures – dislocations will annihilate and the crystal is free of defects. The expression diverges $\langle r^{2}\rangle \to \infty$, if the denominator tends to zero. This happens, when ${\frac {Y\cdot a}{4\pi k_{B}T}}=4$. A diverging distance of dislocations implies, that they are dissociated and do not form a bound pair. The crystal is molten, if several isolated dislocations are thermally excited and the melting temperature $T_{m}$ is given by Young's modulus:

${\frac {Y\cdot a}{k_{B}T_{m}}}=16\pi .$

The dimensionless quantity $16\pi$ is a universal constant for melting in 2D and is independent of details of the system under investigation. This example investigated only an isolated pair of dislocations. In general, a multiplicity of dislocations will appear during melting. The strain field of an isolated dislocation will be shielded and the crystal will get softer in the vicinity of the phase transition; Young's modulus will decrease due to dislocations. In KTHNY theory, this feedback of dislocations on elasticity, and especially on Young's modulus acting as coupling constant in the energy function, is described within the framework of renormalization group theory.

Renormalization of elasticity

If a 2D crystal is heated, virtual dislocation pairs will be excited due to thermal fluctuations in the vicinity of the phase transition. Virtual means, that the average thermal energy is not large enough to overcome (two times) the core-energy and to dissociate (unbind) dislocation pairs. Nonetheless, dislocation pairs can appear locally on very short time scales due to thermal fluctuations, before they annihilate again. Although they annihilate, they have a detectable impact on elasticity: they soften the crystal. The principle is completely analogue to calculating the bare charge of the electron in quantum electrodynamics (QED). In QED, the charge of the electron is shielded due to virtual electron-positron pairs due to quantum fluctuations of the vacuum. Roughly spoken one can summarize: If the crystal is softened due to the presence of virtual pairs of dislocation, the probability (fugacity) $y$ for creating additional virtual dislocations is enhanced, proportional to the Boltzmann factor of the core-energy of a dislocation $y=e^{\frac {E_{C}}{k_{B}T}}$. If additional (virtual) dislocations are present, the crystal will get additionally softer. If the crystal is additionally softer, the fugacity will increase further... and so on and so forth. David Nelson, Bertrand Halperin and independently Peter Young formulated this in a mathematically precise way, using renormalization group theory for the fugacity and the elasticity: In the vicinity of the continuous phase transition, the system becomes critical – this means that it becomes self-similar on all length scales $\gg a$. Executing a transformation of all length scales by a factor of $l$, the energy $E\to E(l)$ and fugacity $y\to y(l)$ will depend on this factor, but the system has to appear identically, simultaneously due to the self similarity. Especially the energy function (Hamiltonian) of the dislocations have to be invariant in structure. The softening of the system after a length scale transformation (zooming out to visualize a larger area implies to count more dislocations) is now covered in a renormalized (reduced) elasticity. The recursion relation for elasticity and fugacity are:

${\frac {dY^{-1}(l)}{dl}}={\frac {3}{2}}\pi y^{2}e^{Y(l)/8\pi }I_{0}{\Big (}Y(l)/8\pi {\Big )}-{\frac {3}{4}}\pi y^{2}e^{Y(l)/8\pi }I_{1}{\Big (}Y(l)/8\pi {\Big )},$

${\frac {dy(l)}{dl}}={\Big (}2-{\frac {Y(l)}{8\pi }}{\Big )}y(l)+2\pi y^{2}e^{Y(l)/16\pi }I_{0}{\Big (}Y(l)/8\pi {\Big )}.$

Similar recursion relations can be derived for the shear modulus and the bulk modulus. $I_{0}$ and $I_{1}$ are Bessel functions, respectively. Depending on the starting point, the recursion relation can run into two directions. $y\to 0$ implies no defects, the ensemble is crystalline. $y\to \infty$, implies arbitrary many defects, the ensemble is fluid. The recursion relation have a fix-point at $y=0$ with $E_{R}/k_{B}T=16\pi$. Now, $E_{R}$ is the renormalized value instead of the bare one. Figure 2 shows Youngs’modulus as function of the dimensionless control parameter $\Gamma$. It measures the ratio of the repelling energy between two particles and the thermal energy (which was constant in this experiment). It can be interpreted as pressure or inverse temperature. The black curve is a thermodynamic calculation of a perfect hexagonal crystal at $T=0$. The blue curve is from computer simulations and shows a reduced elasticity due to lattice vibrations at $T>0$. The red curve is the renormalization following the recursion relations, Young's modulus disappears discontinuously to zero at $16\pi$. Turquoise symbols are from measurements of elasticity in a colloidal monolayer, and confirm the melting point at $Y_{R}=16\pi$.

Interaction between disclinations

Figure 3: Frank's constant in the hexatic phase: it falls below at melting to the isotropic fluid $72/\pi$, and diverges at the transition to the crystal.

The system enters the hexatic phase after the dissociation of dislocations. To reach the isotropic fluid, dislocations (5-7-pairs) have to dissociate into disclinations, consisting of isolated 5-folded and isolated 7-folded particles. Similar arguments for the interaction of disclinations compared to dislocations can be used. Again, disclinations can only be created as pairs due to topological reasons. Starting with the energy $H_{cli}$ as function of distance between two disclinations one finds:

$H_{cli}=-{\frac {F_{A}\cdot \pi }{36}}\sum _{k\neq l}s({\vec {r}}_{k})\cdot s({\vec {r}}_{l})\ln {\frac {\Delta {\vec {r}}_{k,l}}{a}}+E_{s}\cdot N_{cli}.$

The logarithmic term is again dominating. The sign of the interaction gives attraction or repulsion for the winding numbers $+\pi /3$ and $-\pi /3$ of the five- and seven-folded disclinations in a way that charges with opposite sign have attraction. The overall strength is given by the stiffness against twist. The coupling constant $F_{A}$ is called Frank's constant, following the theory of liquid crystals. $E_{s}$ is the discrete energy of a dislocation to dissociate into two disclinations. The squared distance of two disclinations can be calculated the same way, as for dislocations, only the prefactor, denoting the coupling constant, has to be changed accordingly. It diverges for ${\frac {F_{A}\cdot \pi }{36}}=4$. The system is molten from the hexatic phase into the isotropic liquid, if unbound disclinations are present. This transition temperature $T_{i}$ is given by Frank's constant:

${\frac {F_{A}}{k_{B}T_{i}}}=72/\pi .$

$72/\pi$ is again a universal constant. Figure 3 shows measurements of the orientational stiffness of a colloidal monolayer; Frank's constant drops below this universal constant at $T_{i}$.

Critical exponents

Typically, Kosterlitz–Thouless transitions have a continuum of critical points which can be characterised by self-similar grains of disordered and ordered regions. In second order phase transitions, the correlation length measuring the size of those regions diverges algebraically:

$\xi =\xi _{0}{\Big (}{\frac {T-T_{c}}{T_{c}}}{\Big )}^{-\nu }$.

Here, $T_{c}$ is the transition temperature and $\nu$ is a critical exponent. Another special feature of Kosterlitz–Thouless transitions is, that translational and orientational correlation length in 2D diverge exponentially (see also hexatic phase for the definition of those correlation functions):

$\xi =\xi _{0}\cdot e^{{\Big (}{\frac {T-T_{c}}{T_{c}}}{\Big )}^{-\nu }}.$

The critical exponent becomes ${\bar {\nu }}=0{,}36963\dots$ for the diverging translational correlation length at the hexatic – crystalline transition. D. Nelson and B. Halperin predicted, that Frank's constant diverges exponentially with ${\bar {\nu }}$ at $T_{m}$, too. The red curve shows a fit of experimental data covering the critical behaviour; the critical exponent is measured to be ${\bar {\nu }}=0{,}35\pm 0{,}02$. This value is compatible with the prediction of KTHNY theory within the error bars. The orientational correlation length at the hexatic – isotropic transition is predicted to diverge with an exponent $\nu =0{,}5$. This rational value is compatible with mean-field-theories and implies that a renormalization of Frank's constant is not necessary. The increasing shielding of orientational stiffness due to disclinations has not to be taken into account – this is already done by dislocations which are frequently present at $T_{i}$. Experiments measured a critical exponent of $\nu =0{,}5\pm 0{,}03$. KTHNY-theory has been tested in experiment and in computer simulations. For short range particle interaction (hard discs), simulations found a weakly first order transition for the hexatic – isotropic transition, slightly beyond KTHNY-theory.