Kosterlitz–Thouless physics: a review of key issues

In 1969, Thouless was first exposed to the essential ideas vital to theories of transitions driven by topological defects during a visit to Bell Labs, one of the most important centers of theoretical physics in the twentieth century. There he learned of the work by Anderson and colleagues on the transformation of the Kondo problem into a 1D Ising model with 1/r² interactions [4–8] and was asked if he knew anything about this intermediate case. It was known that the model with interactions falling off more slowly than this have long range order [3], or finite magnetization, at low temperatures and that systems with shorter ranged interactions are disordered [2] at all T  >  0.

To answer this challenge, Thouless constructed an argument [1] based on that of Peierls [12, 13] and Landau and Lifshitz [14]. In it, one considers a system of length L with a single domain wall between spins with s  =  +1 and spins with s  =  −1. The entropy is just

$$\begin{eqnarray}&&S={{k}_{\text{B}}}\text{ln}(L/a)+\mathcal{O}(1)\end{eqnarray} \tag{ 1 }$$

and the energy of this isolated domain wall at x is

$$\begin{eqnarray}&&E=J\int_{0}^{x-a/2}\text{d}{{x}_{2}}\int_{x+a/2}^{L}\frac{\text{d}{{y}_{1}}\text{d}{{y}_{2}}}{{{\left(\,{{y}_{1}}-{{y}_{2}}\right)}^{2}}}\approx J\text{ln}(L/a)+\mathcal{O}(1)\end{eqnarray} \tag{ 2 }$$

for the domain wall at $0\ll x\ll L$ well away from the ends of the chain. Thus, the free energy $\Delta F$ of an isolated domain wall is

$$\begin{eqnarray}&&\Delta F(T,L)=\left(J-{{k}_{\text{B}}}T\right)\text{ln}(L/a)+\mathcal{O}(1).\end{eqnarray} \tag{ 3 }$$

The probability of a single domain wall in a system of length $L\gg a$ is

$$\begin{eqnarray}P(t,L)&\sim& \frac{{{\text{e}}^{-\Delta F(L,T)/\left({{k}_{\text{B}}}T\right)}}}{1+{{\text{e}}^{-\Delta F(L,T)/\left({{k}_{\text{B}}}T\right)}}}\\ &\sim& \frac{1}{1+{{(L/a)}^{J/\left({{k}_{\text{B}}}T\right)-1}}}\to \left\{\begin{array}{c} 0,\,\,\,\,\, & {{k}_{\text{B}}}T<J, \\ 1,\,\, & {{k}_{\text{B}}}T>J. \end{array}\right. \end{eqnarray} \tag{ 4 }$$

This implies that one can identify the critical temperature as ${{k}_{\text{B}}}{{T}_{\text{c}}}=J$ since there will be no isolated domain walls in the equilibrium system when $T<{{T}_{\text{c}}}$ so that the system will be ordered. There will be finite domains of reversed spins below ${{T}_{\text{c}}}$ but still a finite magnetization. On the other hand, for $T>{{T}_{\text{c}}}$, there is a finite probability that there will be some isolated free domain walls implying disorder and vanishing magnetization. The original argument [1] went further to show that the magnetization has a discontinuity at ${{T}_{\text{c}}}$ although the transition is continuous. This was later confirmed by more rigorous work [15–17].

The free energy $\Delta F$ of equation (3) for an isolated domain wall in the 1/r² ferromagnetic Ising chain has important implications which were to become relevant to later developments in the defect theory of phase transitions. Domain walls can be regarded as topological defects in an Ising ferromagnet and it is obvious that a finite concentration of these in thermal equilibrium tend to disorder the system. Equation (4) shows that there will be no isolated domain walls in the system for ${{k}_{\text{B}}}T<J$ and there will be long range ferromagnetic order. On the other hand, when ${{k}_{\text{B}}}T>J$, a finite density of these topological defects is present in thermal equilibrium thus disordering the system. This argument is made more exact by Anderson, Yuval and Hamann [6, 8] who treated the 1/r² Ising chain by a renormalization group method. Because this very early version of a renormalization group treatment of a defect driven phase transition was so important to later developments by Thouless and myself [18–20], I review it here.

The first step is to rewrite the partition function Z(K, y, h) of the 1D Ising model with Hamiltonian [6, 8], assuming periodic boundary conditions

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{1}{2}{{K}_{0}}(T)\underset{i>j}{\mathop \sum}\,\frac{{{\left({{s}_{i}}-{{s}_{j}}\right)}^{2}}}{{{(i-j)}^{2}}}+\frac{1}{2}A\underset{i}{\mathop \sum}\,{{\left({{s}_{i}}-{{s}_{i+1}}\right)}^{2}}-h\underset{i}{\mathop \sum}\,{{s}_{i}}\end{eqnarray} \tag{ 5 }$$

in terms of a set of 2n alternating domain walls at ${{r}_{1}},\cdots {{r}_{2n}}$ where $\Delta s=\pm 2$ across a wall. The partition function becomes [6, 8]

$$\begin{eqnarray}&&Z=\underset{n=0}{\overset{\infty}{\mathop \sum}}\,{{y}^{2n}}\underset{i=1}{\overset{2n}{\mathop \prod}}\,\int_{{{r}_{i-1}}+a}^{{{r}_{i+1}}-a}\frac{\text{d}{{r}_{i}}}{a} \\ &&\times \text{exp}\left[{{K}_{0}}\underset{i>j}{\mathop \sum}\,{{q}_{i}}{{q}_{j}}\text{ln}\left(\frac{|{{r}_{i}}-{{r}_{j}}|}{a}\right)-h\underset{i}{\mathop \sum}\,{{q}_{i}}\frac{{{r}_{i}}}{a}\right],\end{eqnarray} \tag{ 6 }$$

where ${{q}_{i}}={{(-)}^{i}}$, ${{K}_{0}}(T)=\frac{J}{{{k}_{\text{B}}}T}$ and ${{y}_{0}}=\text{exp}\left(-\frac{A}{{{k}_{\text{B}}}T}-{{K}_{0}}(T)(1+\gamma )\right)$, with $\gamma =$ Euler's constant, is the domain wall fugacity, which is a measure of the defect density. As will become apparent, the value of the fugacity y₀ affects only the value of the transition temperature ${{T}_{\text{c}}}$ but not its nature. Note that A, the nearest neighbor coupling constant, affects only the the domain wall fugacity y and not the nature of the transition, provided A  >  0. The idea is to proceed perturbatively [4–6, 8] and an expansion parameter is needed. A system of very dilute domain walls is amenable to an expansion in powers of the fugacity y₀ which can be made small by choosing $A\gg 1$ and this expansion is the basis of the renormalization group procedure. This is a somewhat unusual expansion as ${{y}_{0}}\sim {{\text{e}}^{-\mu /{{k}_{\text{B}}}T}}$ is very singular at T  =  0 so, although this expansion is most reasonable at low T, it is not a standard low temperature expansion.

The partition function $Z\left({{K}_{0}},{{y}_{0}},h\right)$ defined by (6) is not analytically accessible but information can be obtained by a renormalization group procedure [5, 6, 8]. This consists of doing a partial trace by integrating over the domain wall positions in the interval $a<{{r}_{i+1}}-{{r}_{i}}<a{{\text{e}}^{\delta l}}$ and rescaling the short distance cutoff or lattice spacing a to obtain a new partition function describing fluctuations on length scales larger than $a(l)=a{{\text{e}}^{l}}$. The rescaled Hamiltonian has the same form as the original Hamiltonian with cutoff a(l) provided the interaction parameters ${{K}_{0}},{{y}_{0}},h$ are replaced by renormalized parameters K(l), y(l), h(l) which are solutions of the renormalization group (RG) flow equations for K(l), y(l), h(l) and the free energy per unit length f(l), [6, 8, 21]

$$\begin{eqnarray}&&\frac{\text{d}K}{\text{d}l}=-4K{{y}^{2}}+\mathcal{O}\left(\,{{y}^{4}}\right),\\ &&\frac{\text{d}y}{\text{d}l}=-(K-1)y+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{\text{d}h}{\text{d}l}=\left(1-2{{y}^{2}}\right)h+\mathcal{O}\left(\,{{y}^{4}}\right),\\ &&\frac{\text{d}f}{\text{d}l}=f+{{y}^{2}}+\mathcal{O}\left(\,{{y}^{4}},{{h}^{2}}\right).\end{eqnarray} \tag{ 7 }$$

It was a major revelation in 1971 when I understood that, by integrating these flow equations and by using a bit of physical intuition, it was possible to extract the behavior of some physically relevant quantities such as the magnetization, susceptibility, specific heat and the correlation length. This was fortunate as calculating the partition function $Z\left({{K}_{0}},{{y}_{0}},h\right)$ was well beyond my ability, but Anderson's scaling ideas provided a more feasible method to obtain information. The flow equations of (7) can be derived with some difficulty but within a finite time and integrating these gives

$$\begin{eqnarray}&&K(l)-\text{ln}\left[K(l)\right]-2{{y}^{2}}(l)={{K}_{0}}-\text{ln}\left[{{K}_{0}}\right]-2y_{0}^{2}=1+C/2\end{eqnarray} \tag{ 8 }$$

where the integration constant C is also a renormalization group invariant.

The transition temperature ${{T}_{\text{c}}}$ is determined by equation (8) with C  =  0 which gives

$$\begin{eqnarray}&&\frac{{{k}_{\text{B}}}{{T}_{\text{c}}}}{J}=1-2{{y}_{0}}+\frac{8}{3}y_{0}^{2}+\mathcal{O}\left(\,y_{0}^{3}\right)\end{eqnarray} \tag{ 9 }$$

which is the dashed line in figure 1. The solid line which is essentially indistinguishable from the dashed line is ${{T}_{\text{c}}}$ from (8) with C  =  0.

Zoom In Zoom Out Reset image size

The flows are plotted in figure 2 and the arrows denote the direction of increasing l. $C\left({{K}_{0}},{{y}_{0}}\right)$ is a constant of integration which depends on ${{K}_{0}}=J/\left({{k}_{\text{B}}}T\right)$ and the fugacity y₀, defined below equation (6). Since the flow line separating regions I and II flows into the point $K(\infty )=1,\,y(\infty )=0$, all systems with initial parameters ${{K}_{0}},{{y}_{0}}$ lying on this separatrix can be interpreted as being at the critical temperature ${{k}_{\text{B}}}{{T}_{\text{c}}}\left(\,{{y}_{0}}\right)/J=K_{\text{c}}^{-1}\left(\,{{y}_{0}}\right)$. For $T={{T}_{\text{c}}}(1+t)$, C  =  −t(4y₀  −  t) to lowest order in t and y₀. Systems in region I flow to some point on the fixed line at $K(\infty )\geqslant 1,\,y(\infty )=0$ which corresponds to a completely ordered system, while a system in region II or III flows to one where y(l) increases to $\mathcal{O}(1)$. This equivalent system has a large density of domain walls which is interpreted as complete disorder or $T>{{T}_{\text{c}}}$. The RG flows are sketched in figure 2 and the transition temperature as a function of y₀ in figure 1 from (8) with C  =  0.

Zoom In Zoom Out Reset image size

One of the big surprises was that the zero field magnetization m has a finite value as $T\to T_{\text{c}}^{-}$ so that there is a discontinuity in $m(T)=m\left({{T}_{\text{c}}}\right)\left(1-b\sqrt{{{T}_{\text{c}}}-T}\right)$ where b is a positive constant [8]. It is interesting to note that this discontinuity in the spontaneous magnetization was first predicted by Thouless [1] and derived from the RG flow equations (7) [8]. It was finally proved rigorously [16] that the model has a finite spontaneous magnetization at low T and later that the magnetization m was discontinuous at the transition [17].

Since the form of the RG flow equations (7) turned out to be very similar to those of the 2D planar rotor model [20], I rederive some of the important results of Anderson and Yuval [8] in more modern language. One of the most remarkable features of the Ising ferromagnet of (5) is that the magnetization jumps discontinuously to zero at ${{T}_{\text{c}}}$. This follows from the solution to the RG equations. Linearize (7) about K  =  1 by writing K  =  1  +  x so that to $\mathcal{O}(x)$,

$$\begin{eqnarray}&&\frac{\text{d}x}{\text{d}l}=-4{{y}^{2}}+\mathcal{O}\left(\,{{y}^{4}},x{{y}^{2}}\right),\\ &&\frac{\text{d}y}{\text{d}l}=-xy+\mathcal{O}\left(\,{{y}^{3}},{{x}^{2}}y\right),\\ &&\Rightarrow \,\,C(t)={{x}^{2}}(l)-4{{y}^{2}}(l)\,\,\,\text{where}\,\,\,t=\frac{T-{{T}_{\text{c}}}}{{{T}_{\text{c}}}}.\end{eqnarray} \tag{ 10 }$$

Here $C(t)=x_{0}^{2}-4y_{0}^{2}=\left({{x}_{0}}-2{{y}_{0}}\right)\left({{x}_{0}}+2{{y}_{0}}\right)$, the small $\left({{x}_{0}},{{y}_{0}}\right)$ expansion of C of equation (8), is a constant of integration which depends on the deviation t from the critical temperature ${{T}_{\text{c}}}$. The RG flows are a set of hyperbolae in the x, y plane in figure 3. The transition is at x  =  y  =  0 and the flows are towards the line $x(\infty )\geqslant 0,y(\infty )=0$ when $C(t)\geqslant 0$ and ${{x}_{0}}\geqslant 2{{y}_{0}}$, which defines region I. When C(t)  <  0, y(l) becomes relevant and flows to $y(l)=\mathcal{O}(1)$, defining region II. Finally, in region III where $C(t)\geqslant 0$ and ${{x}_{0}}\leqslant -2{{y}_{0}}$, the fixed line y  =  0 is unstable and y(l) is relevant, increasing with l. In particular, one identifies the transition as C  =  0 when x(l)  =  2y(l) and all points x(l)  >  2y(l) as $T<{{T}_{\text{c}}}$ since the domain wall fugacity flows to the y  =  0 line which corresponds to a completely ordered system with no spin flips. The high temperature phase is identified with the region of the x, y plane which flows to the region where $y(l)=\mathcal{O}(1)$. This is interpreted as a system with a high concentration of domain walls corresponding to $T>{{T}_{\text{c}}}$ with spin disorder. Thus, one can identify ${{x}_{0}}-2{{y}_{0}}=-t$ so that

$$\begin{eqnarray}&&C(t)=-t\left(4{{y}_{0}}-t\right).\end{eqnarray} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

The unstable separatrix $C=0={{x}_{0}}+2{{y}_{0}}$ corresponding to t  =  4y₀ or $T\gg {{T}_{\text{c}}}$ is just part of the disordered phase, see figure 3.

The solutions to the flow equations (10) fall into three distinct regions of the x(l), y(l) plane: region I, $C(t)\geqslant 0,\,\,\,\,x(l)\geqslant 2y(l)$, region II with $C(t)<0,\,\,\,\,2y(l)>|x(l)|$ and region III where $C>0\,\,\,\,x(l)<0$ with $|x(l)|\geqslant 2y(l)$. As will become apparent, region I corresponds to $T\leqslant {{T}_{\text{c}}}(\,y)$, regions II and III correspond to $T>{{T}_{\text{c}}}(\,y)$. In region I when $C(t)\geqslant 0$ and ${{x}_{0}}>2{{y}_{0}}$, the solutions to (10) are

$$\begin{eqnarray}&&x(l)=\sqrt{C}\frac{1+a{{\text{e}}^{-2\sqrt{C}l}}}{1-a{{\text{e}}^{-2\sqrt{C}l}}},\\ &&y(l)=\sqrt{C}\frac{\sqrt{a}{{\text{e}}^{-\sqrt{C}l}}}{1-a{{\text{e}}^{-2\sqrt{C}l}}},\,\,\,\,\,\text{where}\,\,\,\,a=\frac{{{x}_{0}}-\sqrt{C}}{{{x}_{0}}+\sqrt{C}}.\end{eqnarray} \tag{ 12 }$$

From this, $y\left(l=\infty \right)=0$ so that the scaled system at $l=\infty$ has no defects which implies that, when $T\leqslant {{T}_{\text{c}}}$, one can define a correlation length ${{\xi}_{-}}(T)={{\text{e}}^{{{l}^{\ast}}}}$ where ${{l}^{\ast}}=1/\left(2\sqrt{C}\right)$, so that

$$\begin{eqnarray}{{\xi}_{-}}(T)&=&{{\text{e}}^{1/\left(2\sqrt{C}\right)}}=\text{exp}\left(\frac{1}{2\sqrt{|t|\left(4{{y}_{0}}+|t|\right)}}\right)\\ &=&\left\{\begin{array}{c} \text{exp}\left(\frac{1}{4\sqrt{{{y}_{0}}|t|}}\right),\,\, & 4{{y}_{0}}\gg |t|\geqslant 0, \\ \text{exp}\left(\frac{1}{2|t|}\right),\,\, & 1\gg |t|\gg 4{{y}_{0}}. \end{array}\right. \end{eqnarray} \tag{ 13 }$$

Thus, the true critical region is $t<{{t}^{\ast}}=4{{y}_{0}}$ where one may expect to see ${{\xi}_{-}}(T)$ diverging as ${{\text{e}}^{1/\left(4\sqrt{{{y}_{0}}|t|}\right)}}$ but this region is extremely small since the crossover temperature is proportional to the defect fugacity ${{y}_{0}}\sim {{\text{e}}^{-\mu /\left({{k}_{\text{B}}}T\right)}}$ which is exponentially small. Above this temperature but still in the region $|t|\ll 1$, a gradual crossover will occur to the ${{\text{e}}^{1/\left(2|t|\right)}}$ behavior. This same phenomenon will be found in the 2D planar rotor model and related systems. The correlation length ${{\xi}_{-}}(T)$ can be interpreted as the mean separation between spin flips of opposite sign, or as the size of domains of reversed spins which diverges as $T\to T_{\text{c}}^{-}$.

In this Ising ferromagnet with 1/r² interactions, one can use the RG to show that the magnetization m(T) is discontinuous at ${{T}_{\text{c}}}$ [5]. From equation (7),

$$\begin{eqnarray}&&\frac{\text{d}m}{\text{d}l}=2{{y}^{2}}m+\mathcal{O}\left(\,{{y}^{4}}\right)\\ &&\Rightarrow \,\,\frac{\text{d}}{\text{d}l}\left(K{{m}^{2}}\right)=0\,\,\,\Rightarrow \,\,\,K(0){{m}^{2}}(0)=K(l){{m}^{2}}(l).\end{eqnarray} \tag{ 14 }$$

This shows that K(l)m²(l) is a renormalization group invariant, at least to $\mathcal{O}\left(\,{{y}^{2}}\right)$, so that one can obtain the value of the physical quantity K(0)m²(0) by a judicious choice of the length scale l. When $T\leqslant {{T}_{\text{c}}}$, $y(\infty )=0$ so that the effective system at this scale contains no defects which implies that $m(\infty )={{m}_{\text{sat}}}=1$ and $K(\infty )=1+\sqrt{C}=1+\sqrt{|t|\left(4{{y}_{0}}\,+\,|t|\right)}$ so that

$$\begin{eqnarray}{{K}_{0}}(T){{m}^{2}}(T)&=&K(\infty ){{m}^{2}}(\infty )=1+\sqrt{|t|\left(4{{y}_{0}}\,+\,|t|\right)},\,\,\,\,|t|\ll 1,\\ m(T)&=&\left\{\begin{array}{c} K_{0}^{-1/2}\left({{T}_{\text{c}}}\right)\left(1+\sqrt{|t|\left(4{{y}_{0}}\,+\,|t|\right)}\right)/2,\,\, & -1\ll t\leqslant 0, \\ 0 & t>0. \end{array}\right. \end{eqnarray} \tag{ 15 }$$

Thus, the zero field magnetization m(T) has a discontinuous jump to $m\left(T>{{T}_{\text{c}}}\right)=0$ of size $K_{0}^{-1/2}\left({{T}_{\text{c}}}\right)=\sqrt{\frac{{{k}_{\text{B}}}{{T}_{\text{c}}}}{J}}$, as first argued by Thouless [1], derived by the RG method [6, 8] and proved rigorously by Aizenman et al [17]. As mentioned above, the combination ${{K}_{0}}{{m}^{2}}$ plays the role played by the renormalized stiffness constant in the planar rotor model or in ⁴He in 2D as in figure 4.

Zoom In Zoom Out Reset image size

To this point, the discussion has been limited to region I where $y(l)\to 0$ as $l\to \infty$ so that $T\leqslant {{T}_{\text{c}}}$. When $T>{{T}_{\text{c}}}$ but $t\leqslant 4{{y}_{0}}$, the system is in region II. As the length scale l increases, the fugacity decreases but, at larger l, y(l) begins to increase and continues to increase as l does. This is readily seen from the linearized flow equations of (10) when C  <  0 and $4{{y}^{2}}-{{x}^{2}}=\,|C|>0$, so that

$$\begin{eqnarray}&&\frac{\text{d}y}{\text{d}l}=-xy+\mathcal{O}\left(\,{{y}^{4}}\right),\\ &&\frac{\text{d}x}{\text{d}l}=-4{{y}^{2}}+\mathcal{O}\left(\,{{y}^{4}}\right)=-\left({{x}^{2}}\,+\,|C|\right),\\ &&l=-\int_{{{x}_{0}}}^{x(l)}\frac{\text{d}x}{{{x}^{2}}+|C|}=\frac{1}{\sqrt{|C|}}\left(\text{ta}{{\text{n}}^{-1}}\left(\frac{{{x}_{0}}}{\sqrt{|C|}}\right)-\text{ta}{{\text{n}}^{-1}}\left(\frac{x(l)}{\sqrt{|C|}}\right)\right),\\ &&x(l)=\sqrt{|C|}\text{cot}\left(\sqrt{|C|}l+\text{co}{{\text{t}}^{-1}}\frac{{{x}_{0}}}{\sqrt{|C|}}\right).\end{eqnarray} \tag{ 16 }$$

This is valid for any l, provided only that the RG equations (10) remain valid which requires that the renormalized fugacity y(l)  <  1. Thus, a suitable choice of l is l^* where $x\left({{l}^{\ast}}\right)=-{{x}_{0}}$, $y\left({{l}^{\ast}}\right)={{y}_{0}}$ with ${{x}_{0}}\gg \sqrt{|C|}$ and ${{y}_{0}}\gg \sqrt{|C|}$, so that

$$\begin{eqnarray}&&{{l}^{\ast}}=\frac{2}{\sqrt{|C|}}\text{ta}{{\text{n}}^{-1}}\frac{{{x}_{0}}}{\sqrt{|C|}}=\frac{\pi}{\sqrt{|C|}}=\frac{\pi}{\sqrt{t\left(4{{y}_{0}}-t\right)}}.\\ &&{{\xi}_{+}}(t)={{\text{e}}^{{{l}^{\ast}}}}={{\text{e}}^{\pi /\sqrt{4{{y}_{0}}t}}},\,\,\,\,0<t\ll 4{{y}_{0}}.\end{eqnarray} \tag{ 17 }$$

When t  >  4y₀, the RG equations (10) must be used with $C\left(t,{{y}_{0}}\right)=t\left(t-4{{y}_{0}}\right)>0$ in region III. Most of this region correponds to systems with $T\gg {{T}_{\text{c}}}$ which is expected to display no critical behavior.

These results for the correlation lengths ${{\xi}_{\pm }}(t)$ have a physical interpretation for this 1D Ising model with 1/r² interactions. In the low temperature region I, ${{\xi}_{-}}(T)$ may be interpreted as the mean separation of a pair of opposite domain walls or the size of a domain of reversed spins. In region II, the fugacity y(l) decreases from its initial value y₀ to some minimum but increases for larger l to some value of $\mathcal{O}(1)$ which is outside the domain of validity of the renormalization group equations (7) and (10). Thus, choose l^* in (17) so that $x\left({{l}^{\ast}}\right)=-{{x}_{0}}$ and $y\left({{l}^{\ast}}\right)={{y}_{0}}={{x}_{0}}+\mathcal{O}(C)$ so that $\text{ta}{{\text{n}}^{-1}}\left({{x}_{0}}/\sqrt{|C|}\right)=-\text{ta}{{\text{n}}^{-1}}\left(x\left({{l}^{\ast}}\right)/\sqrt{|C|}\right)=\pi /2$ and ${{\xi}_{+}}(T)$ is interpreted as the typical size of an ordered domain of spins when $T>{{T}_{\text{c}}}$.

The singular part of the free energy f_s is obtained from the flow equations

$$\begin{eqnarray}&&f\left({{l}^{\ast}}\right){{\text{e}}^{-{{l}^{\ast}}}}-f(0)=\int_{0}^{{{l}^{\ast}}}\text{d}l{{\text{e}}^{-l}}{{y}^{2}}(l).\end{eqnarray} \tag{ 18 }$$

The trick is to choose a value of l^* such that the integral and f(l^*) are easily obtained and a suitable choice is ${{l}^{\ast}}=\text{ln}{{\xi}_{\pm }}(T)$. Thus the the singular part of the physical free energy f_s  =  f(l  =  0) is given by

$$\begin{eqnarray}&&{{f}_{s}}={{\text{e}}^{-{{l}^{\ast}}}}f\left({{l}^{\ast}}\right)-\int_{0}^{{{l}^{\ast}}}\text{d}l{{\text{e}}^{-l}}{{y}^{2}}(l)\sim {{\text{e}}^{-{{l}^{\ast}}}}\sim \xi _{\pm }^{-1}(T).\end{eqnarray} \tag{ 19 }$$

This result was first obtained by an RG procedure [8]. This means that the specific heat c_h(t) has the same essential singularity at ${{T}_{\text{c}}}$ as ${{\xi}_{\pm }}(T)$ which is undetectable since all derivatives of the specific heat are zero at ${{T}_{\text{c}}}$ [8].

For both an equilibrium film of ⁴He and a 2D planar rotor magnet, the Hamiltonian is

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{{{K}_{0}}(T)}{2}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{a_{0}^{2}}{{\left(\nabla \theta \left(\mathbf{r}\right)\right)}^{2}}\end{eqnarray} \tag{ 20 }$$

where a₀ is the lattice spacing and

$$\begin{eqnarray}{{K}_{0}}(T)=\left\{\begin{array}{c} \frac{J}{{{k}_{\text{B}}}T} & \,\,\,\,\text{planar}\,\,\,\,\text{rotor}, \\ \frac{{{\hslash}^{2}}\rho _{s}^{0}(T)}{{{m}^{2}}{{k}_{\text{B}}}T} & \,\,\,\,\text{superfluid}\,\,\,\,\text{film}. \end{array}\right. \end{eqnarray}$$

Here, J is the exchange interaction between nearest neighbor unit length spins, $\mathbf{s}\left(\mathbf{r}\right)=\left(\text{cos}\theta \left(\mathbf{r}\right),\text{sin}\theta \left(\mathbf{r}\right)\right)$, $H\left[\left\{\mathbf{s}\right\}\right]=$$\frac{J}{2}\sum_{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{{\left[\mathbf{s}\left(\mathbf{r}\right)-\mathbf{s}\left({{\mathbf{r}}^{\prime}}\right)\right]}^{2}}=J\sum_{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}\left[1-\text{cos}\left(\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)\right)\right]$. At low T, adjacent spins will be almost parallel so that $\text{cos}(\delta \theta )\approx$$1-(\delta \theta )/2$ giving (20). For a ⁴He film,

$$\begin{eqnarray}&&H=\frac{1}{2}\rho _{s}^{0}(T)\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{a_{0}^{2}}{{\mathbf{v}}_{s}}{{\left(\mathbf{r}\right)}^{2}},\end{eqnarray} \tag{ 21 }$$

where ${{\mathbf{v}}_{s}}\left(\mathbf{r}\right)$ is the local superfluid velocity relative to the substrate. Since ${{\mathbf{v}}_{s}}=\frac{\hslash}{m}\nabla \theta$ where θ is the phase of the condensate wave function $\psi \left(\mathbf{r}\right)=|\psi \left(\mathbf{r}\right)|{{\text{e}}^{\text{i}\theta \left(\mathbf{r}\right)}}$. Thus, one argues that (20) is a good description of a planar rotor ferromagnet and of a thin 2D ⁴He film and this form is extensively used.

However, it is important to remember that this Hamiltonian is periodic under $\theta \left(\mathbf{r}\right)\to \theta \left(\mathbf{r}\right)+2\pi$, which allows vortex excitations. Although these have large energy and are therefore improbable, vortices are the only excitations capable of disordering the system or, equivalently, of destroying a uniform superfluid flow. The drift of vortices across the flow to the edges of the system unwinds the phase difference of the superfluid order parameter $\psi ={{\text{e}}^{\text{i}\theta}}$ between the ends of the system. This simple observation implies that these topological excitations must be taken seriously and cannot be ignored because the are very improbable. Smooth fluctuations in the phase are the most probable excitations since these have arbitrarily small energy and are third sound modes in ⁴He films and spin waves in magnets but have no effect on the rigidity of the 2D system. These cannot have anything to do with the phase transition. Very similar but technically more complex arguments hold for a 2D crystal in which phonons are the low energy excitations which do not affect the transition to a fluid while the high energy topological excitations, dislocations, are responsible for the melting transition in 2D.

Note that the Hamiltonian of (20) can be obtained from a Ginsburg-Landau free energy density functional $\mathcal{F}\,[\psi ]$ [45]

$$\begin{eqnarray}&&\mathcal{F}\,\left[\psi \left(\mathbf{r}\right)\right]=\frac{1}{2}|\nabla \psi \left(\mathbf{r}\right){{|}^{2}}+\frac{r(T)}{2}|\psi {{|}^{2}}+\frac{u}{4}|\psi {{|}^{4}}+\cdots \end{eqnarray} \tag{ 22 }$$

Since the system is assumed to be well below the mean field transition temperature, one take $r(T)\ll 0$ and $u\gg 0$. The fluctuations in the amplitude $|\psi \left(\mathbf{r}\right)|$ are negligible and can be ignored, so that $|\psi {{|}^{2}}\,=\,|r(T)|/u$ which is a constant $\mathcal{O}(1)$. Thus, in this phase only approximation, the Hamiltonian (20) results.

Ignoring effects of the $2\pi$ periodicity, at low T the obvious approximation is to regard the system as Gaussian with $-\infty <\theta \left(\mathbf{r}\right)<+\infty$ so that, with ${{s}_{x}}+\text{i}{{s}_{y}}={{\text{e}}^{\text{i}\theta}}=\psi$

$$\begin{eqnarray}&&\langle \psi \left(\mathbf{r}\right){{\psi}^{\ast}}\left({{\mathbf{r}}^{\prime}}\right)\rangle =\langle \mathbf{s}\left(\mathbf{r}\right)\centerdot \mathbf{s}\left({{\mathbf{r}}^{\prime}}\right)\rangle =\langle {{\text{e}}^{\text{i}\left(\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)\right)}}\rangle \sim {{\left(\frac{|\mathbf{r}-{{\mathbf{r}}^{\prime}}|}{{{a}_{0}}}\right)}^{-1/2\pi {{K}_{0}}(T)}}.\end{eqnarray} \tag{ 23 }$$

This result means that, as expected, there is no true long range order at any finite T  >  0 in the planar rotor magnet [23, 46, 47], no Bose Einstein condensation with true long range phase coherence in 2D ⁴He films [48–50] and no long range order in 2D crystals [44, 51]. Also, the power law decay indicates that D  =  2 is the lower critical dimension.

The question of the lower critical dimension in models with a continuous symmetry with a Hamiltonian of the form

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=-K\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,\mathbf{s}\left(\mathbf{r}\right)\centerdot \mathbf{s}\left({{\mathbf{r}}^{\prime}}\right),\,\,\,\,\mathbf{s}=\left({{s}_{1}},{{s}_{2}},\cdots,{{s}_{n}}\right),\,\,\,\,|\mathbf{s}|\,=1,\end{eqnarray} \tag{ 24 }$$

was finally settled by Polyakov [52] by a renormalization group method. The RG recursion relation in dimension d for the temperature variable K⁻¹ is

$$\begin{eqnarray}&&\frac{\text{d}{{K}^{-1}}}{\text{d}l}=-(d-2){{K}^{-1}}+\frac{n-2}{2\pi}{{K}^{-2}}+\mathcal{O}\left({{(d-2)}^{2}},{{K}^{-3}}\right).\end{eqnarray} \tag{ 25 }$$

This yields a critical fixed point at $1/{{K}^{\ast}}=2\pi (d-2)/(n-2)$ which vanishes as $d\to 2$. In exactly d  =  2, the effective temperature K⁻¹(l) increases with l implying that the system is in a disordered phase at any finite K⁻¹  >  0. On the other hand, in exactly 2D, the n  =  2 model is marginal with $\text{d}{{K}^{-1}}/\text{d}l=0$ to all orders in K⁻¹ [53]. The correlation function of (23) agrees with this result and is exact within the gaussian or spin wave approximation.

The same considerations apply to a 2D crystal whose ground state is generally a triangular lattice with sites denoted by $\mathbf{r}$ and deviations from these by $\mathbf{u}\left(\mathbf{r}\right)$ so that the elastic energy can be written in continuum notation as [54]

$$\begin{eqnarray}&&H=\frac{1}{2}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\left[2\mu u_{ij}^{2}\left(\mathbf{r}\right)+\lambda u_{kk}^{2}\left(\mathbf{r}\right)\right].\end{eqnarray} \tag{ 26 }$$

Here, μ and λ are Lamé coefficients and the strain tensor ${{u}_{ij}}=\frac{1}{2}\left(\frac{\partial {{u}_{i}}}{\partial {{r}_{j}}}+\frac{\partial {{u}_{j}}}{\partial {{r}_{i}}}\right)$. Ignoring the periodicity under $\mathbf{u}\to \mathbf{u}+\mathbf{a}$, where $\mathbf{a}$ is a primitive lattice vector, one obtains the Debye Waller factor which is a measure of translational order. Defining the order parameter ${{\rho}_{\mathbf{G}}}\left(\mathbf{r}\right)={{\text{e}}^{\text{i}\mathbf{G}\centerdot \mathbf{r}}}={{\text{e}}^{\text{i}\mathbf{G}\centerdot \mathbf{u}\left(\mathbf{r}\right)}}$ with $\mathbf{G}$ a reciprocal lattice vector, one has

$$\begin{eqnarray}&&{{C}_{\text{G}}}(\mathbf{r})=\langle {{\rho}_{\mathbf{G}}}(\mathbf{r})\rho _{\mathbf{G}}^{\ast}(0)\rangle \sim {{r}^{-{{\eta}_{\mathbf{G}}}(T)}} \\ &&{{\eta}_{\mathbf{G}}}(T)=\frac{{{k}_{\text{B}}}T{{G}^{2}}(3\mu +\lambda )}{4\pi \mu (2\mu +\lambda )}\end{eqnarray} \tag{ 27 }$$

This means that the long range translational order in the ground state is destroyed by the low energy phonon excitations at any T  >  0 in a 2D crystal, in complete analogy to spin waves in the 2D planar rotor ferromagnet.

The planar rotor and related models are special cases in which topological defects play a vital role [18, 19]. Minimizing H of (20), one obtains

$$\begin{eqnarray}&&{{\nabla}^{2}}\theta \left(\mathbf{r}\right)=0.\end{eqnarray} \tag{ 28 }$$

Because of the periodicity under $\theta \left(\mathbf{r}\right)\to \theta \left(\mathbf{r}\right)+2\pi n$, the solution to (28) obeys

$$\begin{eqnarray}&&{{\mathop{\oint}^{}}_{\text{C}}}\text{d}\theta =2\pi n,\,\,\,\,n=0,\pm 1,\pm 2,\cdots \end{eqnarray} \tag{ 29 }$$

where C encloses the vortex core of circulation n. It is easily seen that

$$\begin{eqnarray}&&\,\,\,\,\,\,\,\,\theta \left(\mathbf{r}\right)=\phi \left(\mathbf{r}\right)+\Theta\left(\mathbf{r}\right)=\phi \left(\mathbf{r}\right)+\underset{\mathbf{R}}{\mathop \sum}\,n\left(\mathbf{R}\right)\Theta\left(\mathbf{r},\mathbf{R}\right),\\ &&\Theta\left(\mathbf{r},\mathbf{R}\right)=\text{ta}{{\text{n}}^{-1}}\frac{Y-y}{X-x},\,\,\,n\left(\mathbf{R}\right)=\underset{i}{\mathop \sum}\,{{n}_{i}}\delta \left(\mathbf{R}-{{\mathbf{R}}_{\mathbf{i}}}\right).\end{eqnarray} \tag{ 30 }$$

Here, $\Theta\left(\mathbf{r},\mathbf{R}\right)$ is the phase at $\mathbf{r}=(x,y)$ due to a unit strength vortex at $\mathbf{R}=(X,Y)$ on the dual lattice [19, 55, 56] corresponding to a local energy minimum, n_i is the quantized circulation of the vortex at ${{\mathbf{R}}_{i}}$ and ϕ is the smooth deviation from the local minimum configuration, corresponding to spin waves in a magnet and to third sound in a ⁴He film. After some tedious but straghtforward algebra, one obtains an expression for the energy of a configuration of phase angles in the presence of an arbitrary set of vortices $n\left(\mathbf{R}\right)$ [19],

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{{{H}_{\text{sw}}}}{{{k}_{\text{B}}}T}+\frac{{{H}_{v}}}{{{k}_{\text{B}}}T}+\frac{{{H}_{\text{int}}}}{{{k}_{\text{B}}}T}.\\ &&\frac{{{H}_{\text{sw}}}}{{{k}_{\text{B}}}T}=\frac{1}{2}{{K}_{0}}(T)\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla \phi \left(\mathbf{r}\right)\right)}^{2}},\\ &&\frac{{{H}_{v}}}{{{k}_{\text{B}}}T}=-\pi {{K}_{0}}(T){{\mathop{\int}^{}}_{|\mathbf{R}-{{\mathbf{R}}^{\prime}}|>a}}{{\text{d}}^{2}}\mathbf{R}{{\text{d}}^{2}}{{\mathbf{R}}^{\prime}}n\left(\mathbf{R}\right)n\left({{\mathbf{R}}^{\prime}}\right)\text{ln}\frac{|\mathbf{R}-{{\mathbf{R}}^{\prime}}|}{a}\\ &&-\text{ln}\,{{y}_{0}}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{R}{{n}^{2}}\left(\mathbf{R}\right)+\pi {{K}_{0}}(T){{\left(\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{R}n\left(\mathbf{R}\right)\right)}^{2}}\text{ln}\frac{L}{a},\\ &&{{H}_{\text{int}}}=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\nabla \phi \left(\mathbf{r}\right)\centerdot \nabla \Theta\left(\mathbf{r}\right)=-\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\phi \left(\mathbf{r}\right){{\nabla}^{2}}\Theta\left(\mathbf{r}\right)=0.\end{eqnarray} \tag{ 31 }$$

Note that the spin wave and vortex contributions to H decouple because ${{\nabla}^{2}}\Theta\left(\mathbf{r}\right)=0$.

This decoupling of the spin wave and vortex degrees of freedom is realized explicitly in the Villain model on a discrete lattice [55, 56] where the decoupling of the spin wave and the vortex parts is made explicit. The partition function is

$$\begin{eqnarray}&&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,Z=\underset{\mathbf{r}}{\mathop \prod}\,\int_{0}^{2\pi}\frac{\text{d}\theta \left(\mathbf{r}\right)}{2\pi}\text{exp}\left(\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,V\left(\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)\right)\right),\\ &&{{\text{e}}^{V(\theta )}}=\underset{m=-\infty}{\overset{+\infty}{\mathop \sum}}\,{{\text{e}}^{-\frac{1}{2}K{{(\theta -2\pi m)}^{2}}}},\\ &&m\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)=\tilde{\mathop{m}}\,\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)+m\left(\mathbf{r}\right)-m\left({{\mathbf{r}}^{\prime}}\right);\,\,\,\,\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>\in \mathbf{R}}{\mathop \sum}\,\tilde{\mathop{m}}\,\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)=n\left(\mathbf{R}\right),\\ &&Z\left(K,{{y}_{0}}\right)=\underset{\mathbf{r}}{\mathop \prod}\,\int_{0}^{2\pi}\frac{\text{d}\theta \left(\mathbf{r}\right)}{2\pi}\underset{m\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)}{\mathop \sum}\,{{\text{e}}^{\mathcal{A}[\theta,m]}},\\ &&\,\mathcal{A}[\theta,m]=-\frac{K}{2}\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,{{\left[\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)-2\pi m\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)\right]}^{2}}+\text{ln}{{y}_{0}}\underset{\text{R}}{\mathop \sum}\,{{n}^{2}}\left(\mathbf{R}\right)\end{eqnarray} \tag{ 32 }$$

where the integer variable $m\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)$ lies on the bond between nearest neighbor sites $\mathbf{r}$ and ${{\mathbf{r}}^{\prime}}$. One can write $m\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)=\tilde{m}\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)+m\left(\mathbf{r}\right)-m\left({{\mathbf{r}}^{\prime}}\right)$, by making a discrete gauge transformation, and then define the site variable as $\phi \left(\mathbf{r}\right)=\theta \left(\mathbf{r}\right)+2\pi m\left(\mathbf{r}\right)$. This form of the interaction energy $V(\theta )$ is periodic under $\theta \left(\mathbf{r}\right)\to \theta \left(\mathbf{r}\right)+2\pi n\left(\mathbf{r}\right)$ and approximates ${{K}_{0}}\left(1-\text{cos}\theta \right)$ [55]. The difference between the Villain form and (31) is quantitative but not qualitative and, invoking the concept of universality [57], one expects all models with $2\pi$ periodicity and short range interactions to have identical behavior in the critical region.

One is now faced with computing the grand canonical partition function $Z\left({{K}_{0}},{{y}_{0}}\right)=\sum_{\left\{n\left(\mathbf{R}\right)\right\}}\text{exp}\left(-H/{{k}_{\text{B}}}T\right)$. Here, L is the linear size of the system and a is a short distance cutoff which is the lattice spacing in a magnet and the vortex core size in a superfluid so that $L/a\gg 1$. The last term in H_v of equation (31) requires that $\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{R}n\left(\mathbf{R}\right)=0$ to ensure that the energy is finite in the thermodynamic limit $L\to \infty$. The resulting expression is the energy of a neutral Coulomb plasma [18–20] since the interaction between charged particles in 2D is ${{q}_{i}}{{q}_{j}}\text{ln}\left({{r}_{ij}}/a\right)$ and the problem is reduced to the statistical mechanics of the Coulomb plasma.

The parameter y₀ in (31) is the fugacity of a vortex since the probability of a vortex of circulation $2\pi n$ is proportional to $y_{0}^{{{n}^{2}}}$. To proceed, assume the vortex fugacity ${{y}_{0}}\ll 1$ so that the system has a very small concentration of vortices. This implies that one might try a perturbation expansion to $\mathcal{O}\left(\,y_{0}^{2}\right)$ where ${{y}_{0}}={{\text{e}}^{-{{E}_{\text{c}}}/{{k}_{\text{B}}}T}}$. The vortex core energy ${{E}_{\text{c}}}$ is the interaction energy in the region near the vortex center where the quadratic approximation ${{(\nabla \theta )}^{2}}$ is inadequate because $\nabla \theta$ is large. This also defines the vortex core size a₀ which is a few lattice spacings in a magnet and a few interparticle spacings in a superfluid. Near the center of a vortex, ${{v}_{s}}\sim 1/r$ which diverges as $r\to 0$, which implies that the fluid is normal at the vortex core. Whatever the microscopic structure of a vortex core, there is a large energy density in the core region which is not described by ${{(\nabla \theta )}^{2}}$ and one assumes the contribution to the total energy from this region can be simply described by a core energy ${{E}_{\text{c}}}\gg {{k}_{\text{B}}}T$. Furthermore, assume that the fugacity ${{y}_{0}}\ll 1$ and is a temperature independent constant. The description of the system by equation (31) is a very simplistic one containing only two sets of excitations which are the low energy spin waves and the disordering high energy vortices [18–20]. All others are neglected and it is assumed that this description is adequate at low T and near ${{T}_{\text{c}}}$ since the transition is driven by the vortex excitations.

The first attempt [18] to solve the transition problem in the 2D planar rotor model was based on the free energy of a single isolated vortex of unit circulation. From (31), $\Delta E(L)/\left({{k}_{\text{B}}}T\right)=\pi {{K}_{0}}(T)\text{ln}\left(L/{{a}_{0}}\right)+\mathcal{O}(1)$. The entropy is $\Delta S(L)/{{k}_{\text{B}}}=\text{ln}\left({{L}^{2}}/a_{0}^{2}\right)+\mathcal{O}(1)$ so that the free energy is $\Delta F(L)/$$\left({{k}_{\text{B}}}T\right)=\left(\Delta E(L)-\Delta S(L)\right)/\left({{k}_{\text{B}}}T\right)=\left(\pi {{K}_{0}}(T)-2\right)\text{ln}\left(L/{{a}_{0}}\right)+\mathcal{O}(1)$. Thus, in the thermodynamic limit $L/{{a}_{0}}\to \infty$, the vortex free energy $\Delta F(L)\to -\infty$ when $\pi {{K}_{0}}(T)>2$ and $\Delta F(L)\to +\infty$ so that the probability of an isolated vortex being present in equilibrium is zero. On the other hand, when $\pi {{K}_{0}}(T)<2$, $\Delta F(L)\to -\infty$ so that the probability of an isolated vortex in equilibrium is unity. Of course, there will be a non zero density of finite clusters of vortices of zero total vorticity at all temperatures since the free energy of a finite size zero vorticity cluster is finite. However, the transition temperature ${{T}_{\text{c}}}$ is determined by the renormalized stiffness ${{K}_{\text{R}}}(T)$ such that $\pi {{K}_{\text{R}}}\left({{T}_{\text{c}}}\right)=2$, as discussed in the following section.

David and I congratulated ourselves on finding important new physics but our euphoria soon dissipated. We were informed that Berezinskii [47] had discussed the the vortex driven transition in a superfluid film a year earlier than our paper [18, 19]. Since neither of us knew any Russian we were blissfully unaware of this work while we were developing the basic physics of the vortex driven transition. For some unknown reason, our work seems to have had much greater impact than that of Berezinskii.

The careful reader will have realized that the discussion is rather sloppy and the relation between the integer field $n\left(\mathbf{R}\right)$ and the vorticity is not obvious. Also, the relation between the topological charges $q\left(\mathbf{R}\right)=1/(2\pi )\mathop{\oint}^{}\text{d}\theta =0,\pm 1$ and the effect of boundary conditions are unclear. These issues have been addressed in [58] where the relation between the topological charge $q\left(\mathbf{R}\right)$ and the integer variables $n\left(\mathbf{R}\right)$ is clarified. This issue was first raised by Savit [59]. The effects of boundary conditions in the Coulomb gas language is also discussed in [58] which is important for the Coulomb gas at very low temperature and for finding ground state configurations [60–62].

In this section, I discuss the planar rotor magnet described by equation (31) at low T and fugacity ${{y}_{0}}\ll 1$. The first attempt with David Thouless considered a neutral set of charges $n\left(\mathbf{R}\right)=\pm 1$ interacting by a 2D Coulomb potential (31). Since the fugacity ${{y}_{0}}\ll 1$, charges with larger n are ignored because they are less probable. A neutral pair of unit charges separated by r interacting by a lnr Coulomb interaction will be screened by the polarization of smaller pairs separated by less than r which are screened by even smaller pairs. This leads to a scale dependent dielectric constant $\epsilon (r)={{K}_{0}}/K(r)$ where the force between the charges is $2\pi K(r)/r$. The energy of this pair is

$$\begin{eqnarray}&&E(r)=2\pi \int_{a}^{r}\text{d}{{r}^{\prime}}\frac{K\left({{r}^{\prime}}\right)}{{{r}^{\prime}}}\equiv 2\pi f(r)\text{ln}\left(\frac{r}{a}\right),\\ &&\Rightarrow {{K}^{-1}}(r)=K_{0}^{-1}+4{{\pi}^{3}}y_{0}^{2}\int_{a}^{r}\frac{\text{d}{{r}^{\prime}}}{a}{{\left(\frac{{{r}^{\prime}}}{a}\right)}^{3-2\pi f\left({{r}^{\prime}}\right)}}.\end{eqnarray} \tag{ 33 }$$

We were unable to solve this integral equation so we made an unfortunate unnecessary approximation [19] by replacing $f\left({{r}^{\prime}}\right)$ by f(r) in (33). This was justified by $f\left({{r}^{\prime}}\right)-f(r)\ll 1$, but led to incorrect results.

In an important paper [63], Young demonstrated that our incorrect approximation was unnecessary and presented the correct procedure by defining a scale dependent fugacity

$$\begin{eqnarray}&&y(r)={{y}_{0}}{{\left(\frac{r}{a}\right)}^{2-\pi f(r)}}={{y}_{0}}\text{exp}\left(2\text{ln}(r/a)-\pi \int_{a}^{r}\text{d}{{r}^{\prime}}\frac{K\left({{r}^{\prime}}\right)}{{{r}^{\prime}}}\right),\\ &&{{K}^{-1}}(r)=K_{0}^{-1}+4{{\pi}^{3}}\int_{a}^{r}\frac{\text{d}{{r}^{\prime}}}{{{r}^{\prime}}}{{y}^{2}}\left({{r}^{\prime}}\right).\end{eqnarray} \tag{ 34 }$$

Differentiating these with respect to $l=\text{ln}r$ yields

$$\begin{eqnarray}&&\frac{\text{d}y(l)}{\text{d}l}=\left(2-\pi K(l)\right)y(l),\\ &&\frac{\text{d}{{K}^{-1}}(l)}{\text{d}l}=4{{\pi}^{3}}{{y}^{2}}(l).\end{eqnarray} \tag{ 35 }$$

These are the recursion relations derived earlier [20] by a renormalization group method. However, even if Thouless and I had not made our incorrect approximation in our 1973 paper, but obtained the correct results of (35), we would not have known what to do with them as we knew nothing about scale dependent coupling constants or renormalization groups at that time.

I was not happy with the self consistent results in [19], and realized that the scaling methods of Anderson and coworkers [6, 8] for the 1D Ising model might be applicable to the vortex problem of equation (31). The thinking was that in both the 1D Ising ferromagnet in the previous section and the 2D planar rotor model, the interactions between the important topological point defects are logarithmic and the transitions might be very similar. Of course, in the early 1970's, the importance of the spatial dimension and the range of the interactions [57] was not known to me so, motivated by Anderson's solution of the Kondo problem [6, 8] the renormalization group equations for the 2D Coulomb gas were derived [20] by working with the partition function $Z\left({{K}_{0}},{{y}_{0}}\right)$. This decouples into a product of a gaussian spin wave and a Coulomb gas partition functions,

$$\begin{eqnarray}&&\,\,\,\,\,Z={{Z}_{\text{sw}}}\left({{K}_{0}}\right){{Z}_{\text{c}}}\left({{K}_{0}},{{y}_{0}}\right),\\ &&\,\,\,{{Z}_{\text{c}}}=\underset{n=0}{\overset{\infty}{\mathop \sum}}\,\frac{y_{0}^{2n}}{n{{!}^{2}}}{{\mathop{\int}^{}}_{|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}|>a}}\underset{i=1}{\overset{2n}{\mathop \prod}}\,\frac{{{\text{d}}^{2}}{{\mathbf{R}}_{i}}}{{{a}^{2}}}\text{exp}\left[+\pi {{K}_{0}}\underset{i\ne j}{\mathop \sum}\,{{q}_{i}}{{q}_{j}}\text{ln}\left(\frac{|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}|}{a}\right)\right],\\ &&{{Z}_{\text{sw}}}=\mathop{\int}^{}\mathcal{D}\phi \text{exp}\left(-\frac{1}{2}{{K}_{0}}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla \phi \left(\mathbf{r}\right)\right)}^{2}}\right).\end{eqnarray} \tag{ 36 }$$

${{Z}_{\text{sw}}}\left({{K}_{0}}\right)$ is a trivial Gaussian partition function. In the Coulomb gas partition function ${{Z}_{\text{c}}}\left({{K}_{0}},{{y}_{0}}\right)$, there are n q  =  +1 charges and n q  =  −1 charges. These correspond to vortices of $\pm 1$ circulation which interact by a 2D Coulomb interaction. The integrations over the vortex positions ${{\mathbf{R}}_{i}}$ are restricted by the hard core repulsion $|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}|>a$ where the short distance cut-off length scale a can be taken to be a lattice spacing or the diameter of a vortex core. It turns out that this is quite arbitrary, but it must be kept finite until the end. The infrared cut-off is the system size L must also be kept finite during the calculations to avoid unphysical divergences. It is of some interest to note that both cut-offs, L and a, have clear physical meanings so there is no temptation to take the limit $L\to \infty$ or $a\to 0$ unless it is obvious that these limits make physical sense. Actually, these appear in the combination L/a and the real physical limit is $L/a\to \infty$.

The original derivation of the RG flow equations was done by an extremely involved method [20] based on the method pioneered by Anderson and coworkers [6]. This involves integrating out the short distance degress of freedom in the partition function ${{Z}_{\text{c}}}\left({{K}_{0}},{{y}_{0}}\right)$ of equation (36) and rescaling the cut-off $a\to a{{\text{e}}^{l}}$ resulting in the RG flow equations

$$\begin{eqnarray}&&\frac{\text{d}{{K}^{-1}}}{\text{d}l}=4{{\pi}^{3}}{{y}^{2}}+\mathcal{O}\left(\,{{y}^{4}}\right),\\ &&\frac{\text{d}y}{\text{d}l}=(2-\pi K)y+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{\text{d}f}{\text{d}l}=2f-2\pi {{y}^{2}}+\mathcal{O}\left(\,{{y}^{4}}\right),\end{eqnarray} \tag{ 37 }$$

where the lattice spacing at length scale l is $a(l)=a{{\text{e}}^{l}}$. The interaction is K(l) and the effective fugacity is y(l). The reduced free energy density $f\left({{K}_{0}},{{y}_{0}}\right)$ and specific heat $c\left({{K}_{0}},{{y}_{0}}\right)$ are

$$\begin{eqnarray}&&f\left({{K}_{0}},{{y}_{0}}\right)\equiv -\frac{\mathcal{F}}{{{k}_{\text{B}}}T}=\frac{\text{ln}Z\left({{K}_{0}},{{y}_{0}}\right)}{N},\\ &&c\left({{K}_{0}},{{y}_{0}}\right)={{k}_{\text{B}}}K_{0}^{2}\frac{{{\partial}^{2}}f}{\partial K_{0}^{2}}.\end{eqnarray} \tag{ 38 }$$

They can be integrated to give the renormalization group invariant

$$\begin{eqnarray}&&\text{ln}\left(\frac{\pi K(l)}{2}\right)+\frac{2}{\pi K(l)}-2{{\pi}^{2}}{{y}^{2}}(l)=1+\frac{C\left(t,{{y}_{0}}\right)}{8}\\ &&C\left(t,{{y}_{0}}\right)=-t\left(8\pi {{y}_{0}}-t\right)+\mathcal{O}\left({{t}^{3}},{{t}^{2}}{{y}_{0}},ty_{0}^{2}\right).\end{eqnarray} \tag{ 39 }$$

The RG flows are shown in figure 7 for various values of C. The critical line ${{K}_{\text{c}}}(\,y)$ is the separatrix between regions I and II in figure 7 which is given by equation (39) with C(t, y₀)  =  0. Since $y\ll 1$,

$$\begin{eqnarray}&&\pi {{K}_{\text{c}}}=2+4\pi y+\frac{16}{3}{{\pi}^{2}}{{y}^{2}}+\mathcal{O}\left(\,{{y}^{3}}\right)\end{eqnarray} \tag{ 40 }$$

and ${{K}_{\text{c}}}(\,y)$ from equations (39) and (40) is shown in figure 8

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These recursion relations were later derived by a technically much simpler but more sophisticated method [56, 64] which starts with the Fourier transformed superfluid momentum density ${{\mathbf{g}}^{s}}\left(\mathbf{q}\right)$ correlation function

$$\begin{eqnarray}&&{{C}_{\alpha,\beta}}\left(\mathbf{q},{{K}_{0}},{{y}_{0}}\right)=\langle g_{\alpha}^{s}\left(\mathbf{q}\right)g_{\beta}^{s}\left(-\mathbf{q}\right)\rangle =A(q) \\ &&\frac{{{q}_{\alpha}}{{q}_{\beta}}}{{{q}^{2}}}+B(q)\left({{\delta}_{\alpha \beta}}-\frac{{{q}_{\alpha}}{{q}_{\beta}}}{{{q}^{2}}}\right).\end{eqnarray} \tag{ 41 }$$

Consider the renormalized stiffness constant proportional to the renormalized or measured superfluid density $\rho _{s}^{\text{R}}(T)$, which can be written as a correlation function [65]

$$\begin{eqnarray}&&{{K}_{\text{R}}}(T)=\frac{{{\hslash}^{2}}\rho _{s}^{\text{R}}(T)}{{{m}^{2}}{{k}_{\text{B}}}T}=\underset{q\to 0}{\mathop{\lim}}\,\left(A(q)-B(q)\right)\\&&={{K}_{0}}-4{{\pi}^{2}}K_{0}^{2}\underset{q\to 0}{\mathop{\lim}}\,\frac{\langle n\left(\mathbf{q}\right)n\left(-\mathbf{q}\right)\rangle}{{{q}^{2}}}\\&&={{K}_{0}}+{{\pi}^{2}}K_{0}^{2}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{r}^{2}}\langle n\left(\mathbf{r}\right)n(0)\rangle.\end{eqnarray} \tag{ 42 }$$

Now, the vorticity–vorticity correlation function in (42) can be calculated as a power series in the fugacity y₀ to leading order to obtain

$$\begin{eqnarray}&&\langle n\left(\mathbf{r}\right)n(0)\rangle =-2y_{0}^{2}{{r}^{-2\pi {{K}_{0}}}}+\mathcal{O}\left(\,y_{0}^{4}\right)\\ &&K_{\text{R}}^{-1}=K_{0}^{-1}+4{{\pi}^{3}}y_{0}^{2}\int_{a}^{\infty}\frac{\text{d}r}{a}{{\left(\frac{r}{a}\right)}^{3-2\pi {{K}_{0}}}}+\mathcal{O}\left(\,y_{0}^{4}\right)\end{eqnarray} \tag{ 43 }$$

Since the renormalized stiffness constant ${{K}_{\text{R}}}$ is independent of the cut-off a, one can redefine the cut-off $a\to a{{\text{e}}^{l}}$ to obtain the renormalization group flows of equations (37) which is the same as the self consistent equation (33) which is correct to $\mathcal{O}\left(\,y_{0}^{2}\right)$ only. In principle, equation (43) can be calculated to any desired order. Since ${{K}_{\text{R}}}$ is independent of the cut-off, by redefining $a\to a{{\text{e}}^{l}}$ we can call K(l) and y(l) the effective coupling constat and fugacity at length scale l. This can be expressed as renormalization group equations (37) by standard methods.

This derivation has several advantages over the self consistent one as it makes clear that the RG equations (35) and (37) are valid to $\mathcal{O}\left(\,{{y}^{2}}\right)$ only. Also, this derivation also has the advantage that it uses the perturbative expansion of ${{K}_{\text{R}}}(K,y)$ which yields the exact relation [64]

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left({{K}_{0}},{{y}_{0}}\right)={{K}_{\text{R}}}\left(K(l),y(l)\right)\end{eqnarray} \tag{ 44 }$$

which is the Josephson relation [66] which, in d dimensions, is

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left({{K}_{0}},{{y}_{0}}\right)={{\text{e}}^{(2-d)l}}{{K}_{\text{R}}}\left(K(l),y(l)\right).\end{eqnarray} \tag{ 45 }$$

The first major physical prediction can be made from (37) and (44) at $T_{\text{c}}^{-}$. From figure 8, one sees that $K\left(T_{\text{c}}^{-}\right)$ is an irrelevant variable which flows to the stable fixed point ${{K}_{\text{c}}}\left(l=\infty \right)=2/\pi$ and $y(\infty )=0$. Thus, using (44), we have [64]

$$\begin{eqnarray}&&\frac{{{\hslash}^{2}}\rho _{s}^{\text{R}}\left(T_{\text{c}}^{-}\right)}{{{m}^{2}}{{k}_{\text{B}}}{{T}_{\text{c}}}}={{K}_{\text{R}}}\left({{K}_{0}},{{y}_{0}}\right)={{K}_{\text{R}}}\left(K(\infty ),y(\infty )=0\right)=K(\infty )=\frac{2}{\pi},\\ &&\frac{\rho _{s}^{\text{R}}\left(T_{\text{c}}^{-}\right)}{{{T}_{\text{c}}}}=\frac{2{{m}^{2}}{{k}_{\text{B}}}}{\pi {{\hslash}^{2}}}=3.491\times {{10}^{-8}}~\text{g} ~\text{c}{{\text{m}}^{-2}}~{{\text{K}}^{-1}}.\end{eqnarray} \tag{ 46 }$$

This prediction has been checked experimentally [67, 68] and the data from several different experiments [69–73] is presented in figure 9. It is of interest to notice that the experimental data was obtained and plotted before the authors were aware of the theoretical prediction of [64]. When the theoretical line was inserted, it fitted the plotted data to about 10% accuracy. This is actually a remarkable result because this theoretical prediction is a smoking gun prediction which is an inescapable consequence of the theory. As an aside, note that the stiffness constant ${{K}_{\text{R}}}(T)$ of (44) is sometimes called the helicity modulus $\Upsilon (T)$ [74, 75] which is related to ${{K}_{\text{R}}}(T)$ by

$$\begin{eqnarray}&&\Upsilon (T)={{(\hslash /m)}^{2}}\rho _{s}^{\text{R}}(T)={{k}_{\text{B}}}T{{K}_{\text{R}}}(T).\end{eqnarray} \tag{ 47 }$$

Zoom In Zoom Out Reset image size

The RG equations (37) have a fixed line y(l)  =  0 which is stable for $\pi K(l)\geqslant 2$ and unstable for $\pi K(l)<2$ as shown in figure 7. Thus, one interprets region I of figure 7 as the low temperature phase $T<{{T}_{\text{c}}}$, regions II and III as the disordered $T>{{T}_{\text{c}}}$ phase and the separatrix between regions I and II as the critical line $T={{T}_{\text{c}}}\left(\,{{y}_{0}}\right)$. All points in the ${{K}_{0}}(T),{{y}_{0}}$ plane on the separatrix and in region I flow to the stable portion of the fixed line y  =  0, while all points in regions II and III flow to a region where the effective vortex fugacity y(l) is $\mathcal{O}(1)$ implying the system is disordered.

To proceed further, it is convenient to concentrate on the region near $y=0,\pi K=2$ by defining $\pi K(l)=2+x(l)$ and writing the flow equations (37) to lowest non-trivial order as [20, 56, 64]

$$\begin{eqnarray}&&\frac{\text{d}x}{\text{d}l}=-16{{\pi}^{2}}{{y}^{2}}+\mathcal{O}\left(x{{y}^{2}},{{y}^{4}}\right)\\ &&\frac{\text{d}y}{\text{d}l}=-xy+\mathcal{O}\left({{x}^{2}}y,{{y}^{3}}\right)\end{eqnarray} \tag{ 48 }$$

Integrating these leads to the renormalization group invariant

$$\begin{eqnarray}&&{{x}^{2}}(l)-16{{\pi}^{2}}{{y}^{2}}(l)=C\left(\,{{y}_{0}},t\right)\\ &&=-t\left(8\pi {{y}_{0}}-t\right)+\mathcal{O}\left({{t}^{3}},{{t}^{2}}{{y}_{0}},ty_{0}^{2}\right).\end{eqnarray} \tag{ 49 }$$

The flows are plotted in figure 10 which is nothing but a linearized version of figure 7. Note that the renormalization group equations and consequences are very similar to those for the 1D long range Ising magnet [6, 8]. The line y  =  0 is identified as the spin wave phase since there are no free vortices when the fugacity y  =  0. This phase has finite rigidity and can be identified as the low T quasi ordered phase. On the other hand, when y  >  0 there will be small density of free, mobile vortices so that the rigidity vanishes at long length scales implying this is a high T disordered phase.

Zoom In Zoom Out Reset image size

The (x(l), y(l)) plane shown in figure 10 divides naturally into three separate regions; (I) $C\left(t,{{y}_{0}}\right)\geqslant 0,\,x(l)\geqslant 4\pi y(l)\geqslant 0$, (II) $C\left(t,{{y}_{0}}\right)<0,\,|x(l)|<4\pi y(l)$, and (III) $C\left(t,{{y}_{0}}\right)>0,\,x(l)<$$-4\pi y(l)$, which need separate solutions of equation (48). In region (I) $x(l)\geqslant 4\pi y(l)>0$ and $C\left(t,{{y}_{0}}\right)\geqslant 0$, so that

$$\begin{eqnarray}&&\frac{\text{d}x}{\text{d}l}=-16{{\pi}^{2}}{{y}^{2}},\\ &&\frac{\text{d}y}{\text{d}l}=-xy,\Rightarrow \,{{x}^{2}}(l)-16{{\pi}^{2}}{{y}^{2}}(l)=C>0\\ &&\frac{\text{d}x}{\text{d}l}=C-{{x}^{2}}\,\Rightarrow \,l=\frac{1}{2\sqrt{C}}\int_{{{x}_{0}}}^{x(l)}\text{d}x\left(\frac{1}{x+\sqrt{C}}-\frac{1}{x-\sqrt{C}}\right)\\ &&x(l)=\sqrt{C}\frac{{{x}_{0}}+\sqrt{C}+\left({{x}_{0}}-\sqrt{C}\right){{\text{e}}^{-2l\sqrt{C}}}}{{{x}_{0}}+\sqrt{C}-\left({{x}_{0}}-\sqrt{C}\right){{\text{e}}^{-2l\sqrt{C}}}}.\end{eqnarray} \tag{ 50 }$$

Now examine the behavior of x(l) in various limits by the deviation from ${{T}_{\text{c}}}$ small so that $C\left(t,{{y}_{0}}\right)\ll 1$ and the initial value ${{x}_{0}}\gg \sqrt{C}$. From (50),

$$\begin{eqnarray}x(l)=\left\{\begin{array}{c} \frac{{{x}_{0}}}{1+{{x}_{0}}l}\left(1+\mathcal{O}\left(l\sqrt{C}\right)\right) & \,\,\text{for}\,\,2l\sqrt{C}<1, \\ {} & {} \\ \sqrt{C}\left(1+\mathcal{O}\left({{\text{e}}^{-2l\sqrt{C}}}\right)\right) & \,\,\text{for}\,\,2l\sqrt{C}>1. \end{array}\right. \end{eqnarray} \tag{ 51 }$$

This allows one to define a correlation length for $T\leqslant {{T}_{\text{c}}}$ by ${{\xi}_{-}}(T)={{\text{e}}^{{{l}^{\ast}}}}$ where l^* is defined by $2{{l}^{\ast}}\sqrt{C}=1$, which gives

$$\begin{eqnarray}{{\xi}_{-}}(T)&=&{{\text{e}}^{1/2\sqrt{C}}}=\text{exp}\left(\frac{1}{2\sqrt{|t|\left(8\pi {{y}_{0}}+|t|\right)}}\right)\\ &=&\left\{\begin{array}{c} {{\text{e}}^{1/\left(2\sqrt{8\pi {{y}_{0}}|t|}\right)}} & \,|t|\ll 8\pi {{y}_{0}}, \\ {{\text{e}}^{1/\left(2|t|\right)}} & \,|t|\gg 8\pi {{y}_{0}}. \end{array}\right. \end{eqnarray} \tag{ 52 }$$

The renormalized, or measured, stiffness constant for $T\leqslant {{T}_{\text{c}}}$ is

$$\begin{eqnarray}&&{{K}_{\text{R}}}(T)={{K}_{\text{R}}}\left(K(l),y(l)\right)={{K}_{\text{R}}}\left(K(\infty ),y(\infty )=0\right)\\ &&=K\left(l=\infty \right)=\frac{1}{\pi}\left(2+\sqrt{C}\right)\\ &&=\frac{1}{\pi}\left(2+\sqrt{|t|\left(8\pi {{y}_{0}}\,+\,|t|\right)}\right).\end{eqnarray} \tag{ 53 }$$

Equation (53) seems to give an experimentally measurable prediction for the superfluid density $\rho _{s}^{\text{R}}(T)$ shown in figure 11

$$\begin{eqnarray}&&\frac{\rho _{s}^{\text{R}}(T)}{T}=\frac{{{m}^{2}}{{k}_{\text{B}}}}{\pi {{\hslash}^{2}}}\left(2+\sqrt{|t|\left(8\pi {{y}_{0}}\,+\,|t|\right)}\right).\end{eqnarray} \tag{ 54 }$$

Zoom In Zoom Out Reset image size

However, the superfluid density $\rho _{s}^{\text{R}}(T)$ is, unfortunately, not accessible experimentally as measurements on ⁴He films have to performed with either, or both, the frequency $\omega \ne 0$ or the superfluid velocity ${{\mathbf{v}}_{s}}\ne 0$. Of course, because the system is of finite size or because third sound waves are present $\mathbf{q}\ne 0$, but this effect is not as important. The theoretical prediction here is the $q=\omega ={{\mathbf{v}}_{s}}=0$ component of the dynamical response function $\rho _{s}^{\text{R}}\left(q,\omega,{{\mathbf{v}}_{s}}\right)$ which, unfortunately, is inaccessible to experiment. In principle, one can deduce the behavior of the stiffness constant ${{K}_{\text{R}}}(T)$ in the whole region $T\leqslant {{T}_{\text{c}}}$. For $T\ll {{T}_{\text{c}}}$, ${{K}_{\text{R}}}(T)$ should decrease linearly as T increases until about $10\%$ below ${{T}_{\text{c}}}$ but then vortex fluctuations will take over and in the final approach to ${{T}_{\text{c}}}$, the stiffness constant drops to the universal value $2/\pi$ as $\sqrt{|t|}$ [64]. The theoretical form of ${{\rho}_{s}}(T)$ is shown in figure 10 and an experimental measurement in figure 11. Note that the experimental ${{\rho}_{s}}(T)$ does not drop discontinuously to zero at ${{T}_{\text{c}}}$ but does drop rapidly. This is because the measurements [68] are performed at finite frequency while the theory [64] is for $\omega =0$. A dynamical extension of the static theory is needed [76] to make a more convincing interpretation of the experiments on thin Helium films since the measurements are, of necessity, performed with at least one of $\left(\omega,{{\mathbf{v}}_{s}}\right)$ non zero, while the static theory is for $\omega =0={{\mathbf{v}}_{s}}$ only. The theoretical behavior of ${{\rho}_{s}}(T)$ as a function of temperature T is shown in figure 11 for different coverages.

In region II, C(t, y₀)  <  0 in the small temperature region $0<t<8\pi {{y}_{0}}$ and x(l) can have either sign, so that from (50) [76]

$$\begin{eqnarray}&&16{{\pi}^{2}}{{y}^{2}}={{x}^{2}}\,+\,|C|,\\ &&-l=\int_{{{x}_{0}}}^{x(l)}\frac{\text{d}x}{{{x}^{2}}+|C|}=\frac{1}{\sqrt{|}C|}\left[\text{ta}{{\text{n}}^{-1}}\left(\frac{x(l)}{\sqrt{|C|}}\right)-\text{ta}{{\text{n}}^{-1}}\left(\frac{{{x}_{0}}}{\sqrt{|C|}}\right)\right],\\ &&\frac{x(l)}{\sqrt{|C|}}=\text{cot}\left[\sqrt{|C|}l+\text{co}{{\text{t}}^{-1}}\frac{{{x}_{0}}}{\sqrt{|C|}}\right],\\ &&\frac{4\pi y(l)}{\sqrt{|C|}}=\text{cosec}\left[\sqrt{|C|}l+\text{co}{{\text{t}}^{-1}}\frac{{{x}_{0}}}{\sqrt{|C|}}\right].\end{eqnarray} \tag{ 55 }$$

Here, the upper integration limit is chosen as $x\left({{l}_{+}}\right)=-{{x}_{0}}<0$ and, since the deviation from ${{T}_{\text{c}}}$ is very small when ${{x}_{0}}\gg \sqrt{|C\left(t,{{y}_{0}}\right)|}$, so that y(l₊)  =  y(0) when ${{l}_{+}}=\pi /\sqrt{|C|}$. This allows one to define another length scale when t  >  0

$$\begin{eqnarray}&&{{\xi}_{+}}(t)={{\text{e}}^{{{l}_{+}}}}={{\text{e}}^{\pi /\sqrt{|C|}}}={{\xi}_{-}}{{(t)}^{2\pi}}.\end{eqnarray} \tag{ 56 }$$

For length scales larger than ${{\xi}_{+}}$, the effective fugacity y(l  >  l₊) increases with l to large values outside the validity of our RG equations and the system is interpreted as completely disordered due to the large concentration of free unbound vortices.

The very reasonable assumption that ${{\rho}_{s}}(t>0)=0$ leads to the universal jump in ${{K}_{\text{R}}}(T)$ at ${{T}_{\text{c}}}$. One can show that $\rho _{s}^{\text{R}}=0$ for $T>{{T}_{\text{c}}}$ by defining the generalized stiffness constant ${{K}_{\text{R}}}\left(q,{{K}_{0}},{{y}_{0}}\right)$ at finite q by, following Nelson [77],

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left(q,{{K}_{0}},{{y}_{0}}\right)={{K}_{0}}-4{{\pi}^{2}}K_{0}^{2}\frac{\langle n(q)n(-q)\rangle}{{{q}^{2}}}.\end{eqnarray} \tag{ 57 }$$

Equation (44) gives

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left(q,{{K}_{0}},{{y}_{0}}\right)={{K}_{\text{R}}}\left(q{{\text{e}}^{l}},K(l),y(l)\right)\end{eqnarray} \tag{ 58 }$$

and we can evaluate this approximately at the scale ${{\text{e}}^{{{l}^{\ast}}}}={{\xi}_{+}}$ by a Debye–Hückel [78, 79] approximation

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left(q,{{K}_{0}},{{y}_{0}}\right)=K\left({{l}^{\ast}}\right)-\frac{4{{\pi}^{2}}{{K}^{2}}\left({{l}^{\ast}}\right)}{{{\left(q{{\text{e}}^{{{l}^{\ast}}}}\right)}^{2}}}\langle n\left(q{{\text{e}}^{{{l}^{\ast}}}}\right)n\left(-q{{\text{e}}^{{{l}^{\ast}}}}\right)\rangle \end{eqnarray} \tag{ 59 }$$

where the average $\langle n\left(q{{\xi}_{+}}\right)n\left(-q{{\xi}_{+}}\right)\rangle$ must be calculated in the Coulomb gas ensemble with Hamiltonian

$$\begin{eqnarray}&&\frac{{{H}_{\text{V}}}\left({{l}^{\ast}}\right)}{{{k}_{\text{B}}}T}=\frac{1}{2}\mathop{\int}^{}{{\text{d}}^{2}}k\left(\frac{4{{\pi}^{2}}K\left({{l}^{\ast}}\right)}{{{k}^{2}}}+B\left({{l}^{\ast}}\right)\right)n(k)n(-k).\end{eqnarray} \tag{ 60 }$$

At the scale ${{\text{e}}^{{{l}^{\ast}}}}={{\xi}_{+}}$ the system has large density of vortices with all integer values of ${{n}_{\text{r}}}$, when a reasonable approximation is to integrate instead of summing over the n(q) to obtain

$$\begin{eqnarray}&&{{\langle n\left(q{{\xi}_{+}}\right)n\left(-q{{\xi}_{+}}\right)\rangle}_{{{l}^{\ast}}}}=\frac{1}{4{{\pi}^{2}}K\left({{l}^{\ast}}\right)/{{\left(q{{\xi}_{+}}\right)}^{2}}+B\left({{l}^{\ast}}\right)}.\end{eqnarray} \tag{ 61 }$$

Here, $B\left({{l}^{\ast}}\right)\approx -\text{ln}\,y\left({{l}^{\ast}}\right)=\mathcal{O}(1)$, but its precise value is unimportant. Thus, one obtains

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left(q,T,{{y}_{0}}\right)=\frac{B\left({{l}^{\ast}}\right)K\left({{l}^{\ast}}\right)}{B\left({{l}^{\ast}}\right)+4{{\pi}^{2}}K\left({{l}^{\ast}}\right)/{{\left(q{{\xi}_{+}}\right)}^{2}}}\sim {{\left(q{{\xi}_{+}}\right)}^{2}},\end{eqnarray} \tag{ 62 }$$

as $q{{\xi}_{+}}\to 0$. This means that ${{K}_{\text{R}}}(T)={{K}_{\text{R}}}\left(q=0,T,{{y}_{0}}\right)=0$ when $T>{{T}_{\text{c}}}$ as expected, which completes the argument for the universal jump ${{\rho}_{s}}\left(T_{\text{c}}^{-}\right)/{{T}_{\text{c}}}=3.491\times {{10}^{-9}}~\text{g}~\text{c}{{\text{m}}^{-2}}~{{\text{K}}^{-1}}$ [64].

In the context of the renormalization group, the length scale ${{\xi}_{-}}(t)$ is a measure of when the RG trajectory, which initially flows parallel to the critical separatrix, deviates significantly from the straight line. A more physical interpretation is that ${{\xi}_{-}}(t)$ for $t\leqslant 0$ is the maximum size of neutral bound vortex pairs and there are no vortex pairs with separation $r>{{\xi}_{-}}$. Above ${{T}_{\text{c}}}$, t  >  0 and ${{\xi}_{+}}(t)$ is the largest separation of bound vortex pairs. All pairs with separation $r>{{\xi}_{+}}(t)$ must be regarded as unbound free vortices which can move freely to dissipate superfluid flow or to disorder the system. It is important to understand that ${{\xi}_{+}}(t)<\infty$ for t  >  0 and ${{\xi}_{+}}(t)=\infty$ for $t\leqslant 0$. On the other hand, the shorter length scale ${{\xi}_{-}}(t)<\infty$ both above and below ${{T}_{\text{c}}}$ and diverges as ${{\text{e}}^{\left(b/\sqrt{|t|}\right)}}$ as $|t|\to 0$. It is important to note that both ${{\xi}_{+}}(t)$ and ${{\xi}_{-}}\left(|t|\right)$ diverge exponentially for $8\pi {{y}_{0}}>|t|\geqslant 0$ as (56)

$$\begin{eqnarray}&&{{\xi}_{\pm }}(t)\sim \text{exp}\left({{b}_{\pm }}{{t}^{-1/2}}\right)\,\,\,\text{where}\,\,\,\frac{{{b}_{+}}}{{{b}_{-}}}=2\pi.\end{eqnarray} \tag{ 63 }$$

This behavior is commonly taken to define the critical region of the 2D planar rotor model and of a superfluid film but this region is limited to $|t|\ll 8\pi {{y}_{0}}$ and the vortex fugacity ${{y}_{0}}={{\text{e}}^{-{{E}_{\text{c}}}/{{k}_{\text{B}}}T}}$ is extremely small since the core energy ${{E}_{\text{c}}}\sim J$. The real critical region should be taken to be the range of t for which ${{\xi}_{\pm }}\gg 1$ which is the region $\sqrt{|C(t,{{y}_{0}})|}\ll 1$ where $|C\left(t,{{y}_{0}}\right)|\,=\,|t\left(8\pi {{y}_{0}}-t\right)|\ll 1$ in (63). There is a discussion of this crossover behavior of ${{\xi}_{\pm }}$ and the behavior outside the critical region in [80].

One can calculate other thermodynamic quantities from renormalization group considerations, particularly for $T\leqslant {{T}_{\text{c}}}$ when the RG maps the physical system into a simple Gaussian model when the vortex fugacity $y(l)\to 0$. The two point correlation function G(r, K, y) was calculated by Kosterlitz [20] who obtained

$$\begin{eqnarray}&&G(r,K,y)=\langle {{\text{e}}^{\text{i}\left(\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)\right)}}\rangle \sim {{r}^{-1/4}}{{\left(\text{ln}\,r\right)}^{1/8}},\end{eqnarray} \tag{ 64 }$$

when $T={{T}_{\text{c}}}$. I sketch the derivation [81] below by first noting the scaling equation for d  =  2,

$$\begin{eqnarray}&&\frac{\text{d}G(l)}{\text{d}l}=\eta (l)G(l)\,\,\,\text{where}\,\,\,\eta (l)=\frac{1}{2\pi K(l)}=\frac{1}{4}-\frac{x(l)}{8}+\mathcal{O}\left({{x}^{2}}\right),\\ &&G(0)=G(l)\text{exp}\left(-\int_{0}^{l}\text{d}{{l}^{\prime}}\eta \left({{l}^{\prime}}\right)\right),\\ &&G(l)=G\left(r{{\text{e}}^{-l}},K(l),y(l)\right).\end{eqnarray} \tag{ 65 }$$

To obtain the physical correlation function $G(0)=G\left(r,{{K}_{0}},{{y}_{0}}\right)$ in terms of the physical parameters K₀ and y₀, one needs G(l) which is to be computed from the scaled Hamiltonian. The RG flow equations (48) are valid provided the fugacity y(l) remains sufficiently small, so l in (65) is chosen so that $r{{\text{e}}^{-l}}=1$ and G(1, K(l), y(l)) is needed. This is the correlation function $\langle \text{exp}\left[\text{i}\left(\theta \left(\mathbf{r}+\mathbf{e}\right)-\theta \left(\mathbf{r}\right)\right)\right]\rangle \approx 1$ because, at low T, the correlation function of nearest neighbor spins of unit length is unity. Thus, the singular part of the correlation function is

$$\begin{eqnarray}&&G\left(r,{{K}_{0}},{{y}_{0}}\right)=\text{exp}\left(-\int_{0}^{\text{ln}r}\text{d}l\eta (l)\right).\end{eqnarray} \tag{ 66 }$$

To evaluate this, use $\eta (l)=1/4-x(l)/8+\mathcal{O}\left({{x}^{2}}(l)\right)$ with x(l) from equation (50) to obtain

$$\begin{eqnarray}&&\int_{0}^{l}\text{d}{{l}^{\prime}}x\left({{l}^{\prime}}\right)=\sqrt{C}l+\text{ln}\left[\frac{1-A{{\text{e}}^{-2\sqrt{C}l}}}{1-A}\right]\,\,\,\text{where}\,\,\,A=\frac{{{x}_{0}}-\sqrt{C}}{{{x}_{0}}+\sqrt{C}},\\ &&\approx \sqrt{C}l+\text{ln}\left(1+{{x}_{0}}l\right)\,\,\,\text{when}\,\,\,l\ll \frac{1}{2\sqrt{C}}=\text{ln}{{\xi}_{-}}(T),\\ &&\approx \sqrt{C}l+\text{ln}\frac{{{x}_{0}}}{2\sqrt{C}}\,\,\,\,\,\,\text{when}\,\,\,l\gg \frac{1}{2\sqrt{C}}=\text{ln}{{\xi}_{-}}(T).\end{eqnarray} \tag{ 67 }$$

Thus, the correlation function for $T\leqslant {{T}_{\text{c}}}$ is

$$\begin{eqnarray}&&G\left(\mathbf{r},T\right)\sim {{r}^{-\eta (T)}}{{\left(\text{ln}\,r\right)}^{1/8}}\,\,\,\text{for}\,\,\,1\ll r\ll {{\xi}_{-}}(T),\\ &&G\left(\mathbf{r},T\right)\sim {{r}^{-\eta (T)}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\text{for}\,\,\,1\ll {{\xi}_{-}}(T)\ll r,\end{eqnarray} \tag{ 68 }$$

where the temperature dependent exponent $\eta (T)$ is

$$\begin{eqnarray}\eta (T)=\frac{1}{2\pi {{K}_{\text{R}}}(T)}=\left\{\begin{array}{c} {{k}_{\text{B}}}T/(2\pi J), & \,\,T\ll {{T}_{\text{c}}}, \\ 1/4-(\sqrt{|t|(8\pi {{y}_{0}}\,+\,|t|})/8), & \,\,-1\ll t\leqslant 0. \end{array}.\right.\end{eqnarray} \tag{ 69 }$$

A plot of $\eta (T)$ is shown in figure 12 where one can see the $\sqrt{|t|}$ approach, for $|t|\ll 8\pi {{y}_{0}}$, to the asymptotic value $\eta \left({{T}_{\text{c}}}\right)=1/4$. Note that, for very large $r\gg {{\xi}_{-}}(T)$, the correlation function has a pure power law decay with a temperature dependent exponent, $\eta (T)$, which comes from an effective Gaussian model with coupling constant ${{K}_{\text{R}}}(T)$. For $r\ll {{\xi}_{-}}(T)$, $G\left(\mathbf{r},T\right)$ has the same power law decay but modified by a universal logarithmic correction in (68).

Zoom In Zoom Out Reset image size

The situation for $G\left(\mathbf{r},T\right)$ when $T>{{T}_{\text{c}}}$ is not so clear. Use equation (65) to obtain

$$\begin{eqnarray}&&G\left(r,{{K}_{0}},{{y}_{0}}\right)=G\left(r{{\text{e}}^{-l}},K(l),y(l)\right)\text{exp}\left(-\int_{0}^{l}\text{d}{{l}^{\prime}}\,\eta \left({{l}^{\prime}}\right)\right).\end{eqnarray} \tag{ 70 }$$

When $r\gg {{\xi}_{+}}(T)$, choose $l=\text{ln}{{\xi}_{+}}$ but one still needs $G\left(r/{{\xi}_{+}},K\left(\text{ln}{{\xi}_{+}}\right),y\left(\text{ln}{{\xi}_{+}}\right)\right)$ which is expected to be $\mathcal{O}\left(\text{exp}\left(-r/{{\xi}_{+}}\right)\right)$ since the correlation function should decay exponentially. Thus, when $r\gg {{\xi}_{+}}(T)$,

$$\begin{eqnarray}&&G(r,T)\sim \text{exp}\left(-\frac{r}{{{\xi}_{+}}(T)}\right)\text{exp}\left(-\int_{0}^{\text{ln}{{\xi}_{+}}}\text{d}l\eta (l)\right)\end{eqnarray} \tag{ 71 }$$

which was obtained explicitly by using a duality transformation [56]. On the other hand, when ${{\text{e}}^{l}}\ll {{\xi}_{+}}$, even when $T>{{T}_{\text{c}}}$, the correlation function $G\left(r/{{\xi}_{+}}\leqslant 1,K(l),y(l)\right)=\mathcal{O}(1)$, so that

$$\begin{eqnarray}&&G\left(r,{{K}_{0}},{{y}_{0}}\right)=\text{exp}\left(-\int_{0}^{l}\text{d}{{l}^{\prime}}\eta \left({{l}^{\prime}}\right)\right)\end{eqnarray} \tag{ 72 }$$

where $l=\text{ln}r$. The integral can be carried out by using $\eta (l)=1/4-x(l)/8$, as before, but with x(l) from (55).

$$\begin{eqnarray}&&\int_{0}^{l}\text{d}{{l}^{\prime}}\eta \left({{l}^{\prime}}\right)=\frac{l}{4}-\sqrt{|C|}\int_{0}^{l}\text{d}{{l}^{\prime}}\text{cot}\left(\sqrt{|C|}{{l}^{\prime}}+{{\psi}_{0}}\right), \\ &&\,\,\,{{\psi}_{0}}=\text{co}{{\text{t}}^{-1}}\frac{{{x}_{0}}}{\sqrt{|C|}}\approx \frac{\sqrt{|C|}}{{{x}_{0}}},\\ &&=\frac{l}{4}-\frac{1}{8}\text{ln}\left(\frac{\text{sin}\left(l\sqrt{|C|}+\sqrt{|C|}/{{x}_{0}}\right)}{\text{sin}\left(\sqrt{|C|}/{{x}_{0}}\right)}\right).\end{eqnarray} \tag{ 73 }$$

Thus, we obtain

$$\begin{eqnarray}&&G\left(r\ll {{\xi}_{+}}\right)=\text{exp}\left(-\int_{0}^{l}\text{d}{{l}^{\prime}}\eta \left({{l}^{\prime}}\right)\right)={{\text{e}}^{-l/4}}{{\left(\frac{\text{sin}\left(l\sqrt{|C|}+\sqrt{|C|}/{{x}_{0}}\right)}{\text{sin}\left(\sqrt{|C|}/{{x}_{0}}\right)}\right)}^{1/8}},\\ &&={{\text{e}}^{-l/4}}{{\left({{x}_{0}}l\right)}^{1/8}}={{r}^{-1/4}}{{\left(\text{ln}\,r\right)}^{1/8}},\,\,\,\,\,r={{\text{e}}^{l}}\ll {{\xi}_{-}}(T).\end{eqnarray} \tag{ 74 }$$

When ${{\xi}_{+}}>r>{{\xi}_{-}}$, $\pi >\sqrt{|C|}l>1/2$ where the behavior of the correlation function G(r) of (74) is not so clear. The factor in the brackets does not vary much for $\sqrt{|C|}l\sim \mathcal{O}(1)$ and the cross over to the behavior for $r\gg {{\xi}_{+}}$ is unclear.

The lowest order flow equations (37) yield a remarkable amount of information and appear to yield exact results for $T\leqslant {{T}_{\text{c}}}$ because y(l) is an irrelevant variable so that the system becomes Gaussian as $l\to \infty$. However, one sees that this marginally irrelevant variable is responsible for the ${{\left(\text{ln}r\right)}^{1/8}}$ factor in the two point correlation function G(r). Do higher order terms in the RG flow equations give a noticeable correction to the expressions for G(r)? This calculation has been performed [82, 83] who derived flow equations to next order in y(l) and found the correlation function at ${{T}_{\text{c}}}$ as

$$\begin{eqnarray}&&G(r)\sim {{r}^{-1/4}}{{\left(\text{ln}r\right)}^{1/8}}\left(1+\frac{1}{16}\frac{\text{ln}\,\text{ln}r}{\text{ln}r}\right).\end{eqnarray} \tag{ 75 }$$

Note that these results imply the finite size scaling form for the susceptibility at criticality for a $L\times L$ system

$$\begin{eqnarray}&&\chi (L)\sim {{L}^{7/4}}{{\left(\text{ln}L\right)}^{1/8}}.\end{eqnarray} \tag{ 76 }$$

This has been verified to great accuracy by numerical simulations on large systems up to L  =  2¹⁶ [84, 85].

I now turn to thermal properties, in particular the specific heat of the planar rotor model in 2D and of a thin film of ⁴He. In general, the specific heat is the most obvious quantity to measure since it usually displays a singularity, albeit a relatively weak one, at a phase transition. For example, the 2D Ising ferromagnet has a specific heat singularity of the form $c(T)\sim \text{ln}|t|$ while the susceptibility $\chi (T)\sim |t{{|}^{-1.75}}$. However, a measurement of the specific heat requires no special knowledge of the particular system such as the form of the order parameter or of the ordering field. For example, for an antiferromagnet the uniform susceptibility, or the response to a uniform magnetic field, has no dramatic divergence at ${{T}_{\text{c}}}$ but just a minor kink. The dramatic effect is in the staggered susceptibility which is difficult to measure as it is the response to a staggered magnetic field conjugate to the staggered magnetization. The specific heat of a ferromagnet and an antiferromagnet are the same and are measured by the same methods which require no knowledge of the symmetry of the system or the nature of the order parameter. In general, if a transition exists, a specific heat measurement will find it, except in cases such as the systems being considered here.

One needs the free energy $f\left({{K}_{0}},{{y}_{0}}\right)$ from which the specific heat is obtained by differentiation with respect to K₀ as in (38). Integrating the flow equation (37) for f(l) yields [86, 87]

$$\begin{eqnarray}&&f\left({{K}_{0}},{{y}_{0}}\right)={{\text{e}}^{-2l}}f\left(K(l),y(l)\right)-2\pi \int_{0}^{l}\text{d}{{l}^{\prime}}{{\text{e}}^{-2{{l}^{\prime}}}}{{y}^{2}}\left({{l}^{\prime}}\right).\end{eqnarray} \tag{ 77 }$$

When $T\leqslant {{T}_{\text{c}}}$, one can take the upper limit $l\to \infty$ because the RG equations (37) are valid for all $0\leqslant l\leqslant \infty$ and ${{\text{e}}^{-2l}}f\left(K(l),y(l)\right)\to 0$ so that

$$\begin{eqnarray}&&f\left({{K}_{0}},{{y}_{0}}\right)=2\pi \int_{0}^{\infty}\text{d}l\,{{\text{e}}^{-2l}}{{y}^{2}}(l).\end{eqnarray} \tag{ 78 }$$

When $T>{{T}_{\text{c}}}$, the fugacity y(l) is relevant and increases so that the RG equations become invalid for l  >  l^*, where l^* is chosen judiciously [86, 87] and

$$\begin{eqnarray}&&f\left({{K}_{0}},{{y}_{0}}\right)={{\text{e}}^{-2{{l}^{\ast}}}}f\left(K\left({{l}^{\ast}}\right),y\left({{l}^{\ast}}\right)\right)+2\pi \int_{0}^{{{l}^{\ast}}}\text{d}l{{\text{e}}^{-2l}}\,{{y}^{2}}(l),\end{eqnarray} \tag{ 79 }$$

The idea is to use the RG flow equations (37) for $l\leqslant {{l}^{\ast}}$ and estimate $f\left(K\left({{l}^{\ast}}\right),y\left({{l}^{\ast}}\right)\right)$ by an improved Debye–Hückel approximation. Detailed discussions of the matching point and calculations of $f\left(K\left({{l}^{\ast}}\right),y\left({{l}^{\ast}}\right)\right)$ are found in [86, 87]. The result is shown in figure 13. At the transition at $K_{\text{c}}^{-1}$, the specific heat has an unobservable essential singularity of the form $c(T)\sim \xi _{+}^{-2}(T)$ where ${{\xi}_{+}}(T)\sim {{\text{e}}^{b/\sqrt{|t|}}}$ at the critical point $K_{\text{c}}^{-1}$ with a non universal peak at much higher temperature. Since the specific heat is mainly due to the vortex degrees of freedom, one can understand the peak as being due to the entropy released by the unbinding of vortex pairs of decreasing size as the temperature increases. Since a very small density of very large pairs unbind at ${{T}_{\text{c}}}$, the specific heat is very small there. The density of vortex pairs increases as their size decreases and these unbind at higher T with a corresponding larger entropy release and larger specific heat. The maximum will be when the largest density of the smallest vortex pairs unbind, which is estimated to be $40\%$ above ${{T}_{\text{c}}}$ as shown in figure 13 [87].

Zoom In Zoom Out Reset image size

Rather surprisingly, as pointed out by Chui and Weeks [88, 89], the 2D planar rotor model is intimately related to models for the roughening of facets of crystals in thermodynamic equilibrium as temperature is increased. We represent the particle positions of an ideal flat facet as a 2D lattice of sites labelled by $\mathbf{r}$. Even if the equilibrium configuration of a particular crystal is a smooth planar facet which may appear in a crystal grown at a particular temperature, but it also may not because of competing out of equilibrium facets with faster growth rates. It is well known that a facet with screw dislocations terminating in it will grow very much faster than the ideal equilibrium facet. There is some hope that a quantum system such as a ⁴He crystal immersed in superfluid may equilibrate sufficiently fast to be seen.

The partition function of the 2D Coulomb gas representation of the planar rotor model of (31) is written as [56]

$$\begin{eqnarray}&&Z\left({{K}_{0}},{{y}_{0}}\right)=\underset{n\left(\mathbf{r}\right)}{\mathop \sum}\,\text{exp}\left(2{{\pi}^{2}}{{K}_{0}}\underset{\mathbf{r},{{\mathbf{r}}^{\prime}}}{\mathop \sum}\,G\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)n\left(\mathbf{r}\right)n\left({{\mathbf{r}}^{\prime}}\right)+\text{ln}{{y}_{0}}\underset{\mathbf{r}}{\mathop \sum}\,{{n}^{2}}\left(\mathbf{r}\right)\right) \\ &&=\underset{n\left(\mathbf{r}\right)}{\mathop \sum}\,\underset{\mathbf{r}}{\mathop \prod}\,\int_{-\infty}^{+\infty}\frac{\text{d}\phi \left(\mathbf{r}\right)}{2\pi}\text{exp}\left(-\frac{1}{2{{K}_{0}}}\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,{{\left(\phi \left(\mathbf{r}\right)-\phi \left({{\mathbf{r}}^{\prime}}\right)\right)}^{2}}\,\,\right. \\ &&\,\,\,\,\,\,\,\,\,\,\,+\underset{\mathbf{r}}{\mathop \sum}\,\left(2\pi \text{i}n\left(\mathbf{r}\right)\phi \left(\mathbf{r}\right)+\text{ln}\,{{y}_{0}}{{n}^{2}}\left(\mathbf{r}\right)\right),\\ &&G\left(\mathbf{r},{{\mathbf{r}}^{\prime}}\right)=\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{q}}{{{(2\pi )}^{2}}}\frac{{{\text{e}}^{\text{i}\mathbf{q}\centerdot \left(\mathbf{r}-{{\mathbf{r}}^{\prime}}\right)}}}{{{q}^{2}}}=\frac{1}{2\pi}\text{ln}\frac{|\mathbf{r}-{{\mathbf{r}}^{\prime}}|}{L},\end{eqnarray} \tag{ 80 }$$

with L the facet size. Consider two cases (i) when the vortex fugacity y₀  =  1, use the identity [90]

$$\begin{eqnarray}&&\underset{n=-\infty}{\overset{+\infty}{\mathop \sum}}\,{{\text{e}}^{2\pi \text{i}n\phi}}=\underset{h=-\infty}{\overset{+\infty}{\mathop \sum}}\,\delta (\phi -h)\end{eqnarray}$$

to obtain the discrete Gaussian model [56, 88, 89, 93]

$$\begin{eqnarray}&&Z\left({{K}_{0}}\right)=\underset{h\left(\mathbf{r}\right)=-\infty}{\overset{+\infty}{\mathop \sum}}\,\text{exp}\left(-\frac{1}{2{{K}_{0}}}\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,{{\left(h\left(\mathbf{r}\right)-h\left({{\mathbf{r}}^{\prime}}\right)\right)}^{2}}\right).\end{eqnarray} \tag{ 81 }$$

This is interpreted as a model of the facet in which the height of the facet at $\mathbf{r}$, relative to a reference plane, is $h\left(\mathbf{r}\right)$ and neighboring columns of atoms differ by an integer number of lattice spacings. The temperature of the roughening model is proportional to K₀ while the temperature of the dual planar rotor is $K_{0}^{-1}$. When K₀ is small, it is clear that all columns tend to have the same height $h\left(\mathbf{r}\right)$ and the facet is smooth and flat. On the other hand, at higher temperature with K₀ large, the column heights $h\left(\mathbf{r}\right)$ will vary with position leading to a rough facet with the correlation function $\langle {{\left(h\left(\mathbf{r}\right)-h(0)\right)}^{2}}\rangle \sim \text{ln}\,r$.

(ii) ${{y}_{0}}\ll 1$ when the partition function becomes

$$\begin{eqnarray}&&Z\left({{K}_{0}},{{y}_{0}}\right)=\mathop{\int}^{}\underset{\mathbf{r}}{\mathop \prod}\,\frac{\text{d}\phi \left(\mathbf{r}\right)}{2\pi}\left[1+2{{y}_{0}}\text{cos}\left(2\pi \phi \left(\mathbf{r}\right)\right)\right]\\ &&\,\,\,\,\,\,\,\times \text{exp}\left[-\frac{1}{2{{K}_{0}}}\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,\left(\phi \left(\mathbf{r}\right)-\phi {{\left({{\mathbf{r}}^{\prime}}\right)}^{2}}\right],\right. \\ &&=\mathop{\int}^{}\underset{\mathbf{r}}{\mathop \prod}\,\frac{\text{d}\phi \left(\mathbf{r}\right)}{2\pi}\text{exp}\left[-\frac{1}{2{{K}_{0}}}\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,{{\left(\phi \left(\mathbf{r}\right)-\phi \left({{\mathbf{r}}^{\prime}}\right)\right)}^{2}}\right. \\ &&\left. \,\,\,\,\,\,\,+2{{y}_{0}}\underset{\mathbf{r}}{\mathop \sum}\,\text{cos}\left(2\pi \phi \left(\mathbf{r}\right)\right)\right].\end{eqnarray} \tag{ 82 }$$

To convert this to the sine-Gordon roughening model [88, 89, 93], the height variable is defined by $h\left(\mathbf{r}\right)=b\phi \left(\mathbf{r}\right)$ where $h\left(\mathbf{r}\right)$ is the local height at $\mathbf{r}$ in steps of spacing b between lattice planes. The sine-Gordon representation of the roughening problem has the Hamiltonian

$$\begin{eqnarray}&&H(\gamma,u)={{H}_{0}}+{{H}_{u}}=\frac{\gamma}{2}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla h\left(\mathbf{r}\right)\right)}^{2}}-u\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{{{a}^{2}}}\text{cos}\left(\frac{2\pi h\left(\mathbf{r}\right)}{b}\right),\end{eqnarray} \tag{ 83 }$$

where u  =  2y₀ when ${{y}_{0}}\ll 1$ and a is the short distance cut off.

Flow equations for $\gamma (l)$ and u(l) [88, 93] are derived in a similar way to K(l) and y(l) in the vortex system by defining a renormalized stiffness ${{\gamma}_{\text{R}}}(T)$ as

$$\begin{eqnarray}&&F\left(\mathbf{v}\right)-F(0)=\frac{1}{2}\Omega {{\gamma}_{\text{R}}}(T){{\mathbf{v}}^{2}}\,\,\,\text{with}\,\,\,\,\mathbf{v}=\frac{1}{\Omega }\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\left(\nabla h\left(\mathbf{r}\right)\right),\end{eqnarray} \tag{ 84 }$$

where $F\left(\mathbf{v}\right)$ is the free energy of a surface of area $\Omega$ and $\mathbf{v}$ is the average slope of the surface. Define $h\left(\mathbf{r}\right)=\mathbf{v}\centerdot \mathbf{r}+{{h}^{\prime}}\left(\mathbf{r}\right)$ with $\mathbf{v}=\left({{v}_{1}},{{v}_{2}}\right)$ so that the sine-Gordon Hamiltonian is

$$\begin{eqnarray}&&H\left(\mathbf{v}\right)=\frac{1}{2}\Omega \gamma {{\mathbf{v}}^{2}}+{{H}_{0}}\left[{{h}^{\prime}}\left(\mathbf{r}\right)\right]+{{H}_{u}}\left[{{h}^{\prime}}\left(\mathbf{r}\right)+\mathbf{v}\centerdot \mathbf{r}\right].\end{eqnarray} \tag{ 85 }$$

Now expand $F\left(\mathbf{v}\right)=-{{k}_{\text{B}}}T\text{lnTr}\,\text{exp}\left[-H\left(\mathbf{v}\right)/{{k}_{\text{B}}}T\right]$ to $\mathcal{O}\left({{\mathbf{v}}^{2}}\right)$ to obtain

$$\begin{eqnarray}&&{{\gamma}_{\text{R}}}(T)=\gamma +\frac{1}{2\Omega }\underset{\alpha,\beta}{\mathop \sum}\,\frac{{{\partial}^{2}}}{\partial {{v}_{\alpha}}\partial {{v}_{\beta}}}\left({{\langle {{H}_{u}}\rangle}_{0}}-\frac{1}{2{{k}_{\text{B}}}T}\left[{{\langle H_{u}^{2}\rangle}_{0}}-\langle {{H}_{u}}\rangle _{0}^{2}\right]\right)+\cdots \\ &&=\gamma +\frac{{{k}_{\text{B}}}T}{8}{{\left(\frac{2\pi}{b}\right)}^{2}}{{u}^{2}}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{{{a}^{2}}}{{\left(\frac{r}{a}\right)}^{2}}{{\langle \text{cos}\left(\frac{2\pi}{b}\left[{{h}^{\prime}}\left(\mathbf{r}\right)-{{h}^{\prime}}(0)\right]\right)\rangle}_{0}},\\ &&K_{\text{R}}^{-1}={{K}^{-1}}+{{\pi}^{3}}{{u}^{2}}\int_{a}^{\infty}\frac{\text{d}r}{a}{{\left(\frac{r}{a}\right)}^{3-2\pi K}},\end{eqnarray} \tag{ 86 }$$

where $K\equiv {{k}_{\text{B}}}T/\left(\gamma {{b}^{2}}\right)$ and we have used $\langle \text{cos}\left[(2\pi /b)\left({{h}^{\prime}}\left(\mathbf{r}\right)-\right. \right.$$\left. \left. {{h}^{\prime}}(0)\right)\right]{{\rangle}_{0}}={{(r/a)}^{-2\pi K}}$.

The last line of (86) yields renormalization group flow equations exactly as for the planar rotor model by rescaling the cut off $a\to a\,{{\text{e}}^{l}}$,

$$\begin{eqnarray}&&\frac{\text{d}{{K}^{-1}}(l)}{\text{d}l}={{\pi}^{3}}{{u}^{2}}(l)+\mathcal{O}\left({{u}^{4}}\right)\\ &&\frac{\text{d}u(l)}{\text{d}l}=\left(2-\pi K(l)\right)u(l)+\mathcal{O}\left({{u}^{3}}\right).\end{eqnarray} \tag{ 87 }$$

This is of exactly the same form as the flow equations (37) except for the interpretation of the coupling constant K. For the planar rotor model, $K(T)={{(\hslash /m)}^{2}}\rho _{s}^{0}/\left({{k}_{\text{B}}}T\right)$ while, for the roughening model, $K(T)={{k}_{\text{B}}}T/\left(\gamma {{b}^{2}}\right)$. The T dependence is inverted, as expected, because the models are duals of each other [56] and the low temperature phase of one is the high temperature phase of the other. Thus, in the high T phase of the roughening model, the coefficient u(l) of $\text{cos}\left(2\pi h\left(\mathbf{r}\right)\right)/b$ is irrelevant and the roughening Hamiltonian becomes a simple Gaussian model at long length scales with

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{1}{2}{{K}_{\text{R}}}(T)\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla h\left(\mathbf{r}\right)\right)}^{2}}.\end{eqnarray} \tag{ 88 }$$

The integer spacing between the crystal planes becomes irrelevant and the interface becomes rough in the sense that the interface width $\langle {{h}^{2}}\left(\mathbf{r}\right)\rangle \sim \text{ln}L$ in a size L system.

One can immediately infer some exact results for the roughening model by exploiting the relationship with the planar rotor model. At the roughening transition, the universal relation of equation (46) becomes

$$\begin{eqnarray}&&{{K}_{\text{R}}}\left(T_{\text{R}}^{+}\right)=\frac{{{k}_{\text{B}}}{{T}_{\text{R}}}}{{{\gamma}_{\text{R}}}\left(T_{\text{R}}^{+}\right){{b}^{2}}}=\frac{2}{\pi}.\end{eqnarray} \tag{ 89 }$$

This follows directly from the planar rotor result and the duality relation and implies that facets with the largest lattice spacing b will have the highest roughening temperature ${{T}_{\text{R}}}$, assuming the surface stiffness γ does not depend much on b. However, to convert these theoretical predictions into experimentally measurable quantities, one has to remember that the facets are the faces of a 3D real crystal, and the equilibrium theory of a constant volume crystal is required.

The basic idea is that the crystal shape is determined by minimizing the total surface free energy at fixed volume [91]. A flat facet with $-{{x}_{0}}<x<+{{x}_{0}}$, in terms of reduced coordinates $\tilde{x}=x/L,\,\tilde{y}=y/L$ and $\tilde{h}=h/L$, obeys

$$\begin{eqnarray}\tilde{\mathop{h}}\,\left(\tilde{\mathop{x}}\,,0\right)=\left\{\begin{array}{c} {{\tilde{\mathop{h}}\,}_{0}}-A{{\left(|\tilde{\mathop{x}}\,|-{{\tilde{\mathop{x}}\,}_{0}}\right)}^{3/2}} & \,\,\,\,|\tilde{\mathop{x}}\,|>{{\tilde{\mathop{x}}\,}_{0}}, \\ {{\tilde{\mathop{h}}\,}_{0}} & \,\,\,\,|\tilde{\mathop{x}}\,|<{{\tilde{\mathop{x}}\,}_{0}}. \end{array}\right. \end{eqnarray} \tag{ 90 }$$

Here, the reduced facet size ${{\tilde{\mathop{x}}\,}_{0}}\sim b/{{\xi}_{+}}(T)\sim \text{exp}\left(-B/\sqrt{\left(T-{{T}_{\text{R}}}\right)/{{T}_{\text{R}}}}\right)$ and ${{\tilde{h}}_{0}}={{f}_{0}}/\gamma$, with f₀ the surface free energy density.

For $T>{{T}_{\text{R}}}$, the facet has vanished, ${{\tilde{x}}_{0}}=0$, and the reduced height $\tilde{h}\left(\tilde{x},0\right)$ is

$$\begin{eqnarray}&&\tilde{h}\left(\tilde{x},0\right)={{\tilde{h}}_{0}}-\frac{\gamma}{{{\gamma}_{\text{R}}}(T)}{{\tilde{x}}^{2}}.\end{eqnarray} \tag{ 91 }$$

The universal jump in the stiffness constant translates, for a crystal facet, to the reduced curvature of the rough surface ${{K}_{\text{R}}}(T)=(L/R){{k}_{\text{B}}}T/\left(\gamma {{b}^{2}}\right)=2/\pi$ at $T=T_{\text{R}}^{+}$ which jumps discontinuously to zero at $T_{\text{R}}^{-}$. Here, L is the facet size and R is the radius of curvature [91, 93]. Measuring the predicted jump in the curvature of a disappearing facet is not easy and, to the best of my knowledge, not all the universal predictions have been verified experimentally, particularly that the facet size vanishes as $\xi _{+}^{-1}(T)$ [92, 94]. Attempts to measure the reduced curvature ${{K}_{\text{R}}}(T)$ of (89) have been made on hcp crystals in contact with superfluid [95], but with mixed success. Despite being the system most likely to be in equilibrium, the measurements are consistent with theory but suffer seriously from lack of equilibration.

One motivation for developing a theory of transitions in 2D was the experimental results for very thin films of ⁴Helium adsorbed on a substrate [36, 37, 68]. This system is closely related to the 2D planar rotor model and the predictions from the theoretical model are directly applicable to the experiments provided one can argue that substrate effects are irrelevant in the renormalization group sense. Of course, the substrate is essential as a fluid film a few atomic lengths thick, a monolayer crystal or magnet cannot exist in the absence of a substrate. About the only systems which can exist without a substrate are suspended liquid crystal films [96–103] and a monolayer of graphene [104]. Therefore, it is essential to assess the effects of all possible perturbations due to a substrate on the ideal two dimensional system of section 2.2.

Before discussing some of the substrate effects on 2D planar rotor systems, we can see why superfluidity in thin ⁴Helium films can be understood quantitatively by the idealized 2D theory which ignores the effects of the far from perfect Mylar substrate. This can be modeled by a random potential coupled to the local mass density $|\psi \left(\mathbf{r}\right){{|}^{2}}$. The order parameter is the quantum mechanical condensate wave function $\psi \left(\mathbf{r}\right)=\,|\psi \left(\mathbf{r}\right)|{{\text{e}}^{\text{i}\theta \left(\mathbf{r}\right)}}$ and the random substrate potential cannot couple directly to the phase of the wave function but only to its amplitude and is irrelevant. This justifies using the ideal planar rotor model for a superfluid Helium film on a Mylar substrate.

Many other systems are influenced by the substrate particularly a planar rotor magnet in a symmetry breaking crystal field such as a planar ferromagnet on a periodic substrate described by a Hamiltonian

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{1}{2}K(T)\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla \theta \left(\mathbf{r}\right)\right)}^{2}}-{{h}_{p}}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\text{cos}\left(\,p\theta \left(\mathbf{r}\right)\right).\end{eqnarray} \tag{ 92 }$$

where h_p is the strength of the p-fold symmetry breaking field. When $|{{h}_{p}}|\gg 1$, the free energy is minimized when $\text{cos}\theta =+1\,\,\Rightarrow \,\,\theta \left(\mathbf{r}\right)=2\pi n/p$ for h_p  >  0, and, when h_p  <  0, $\theta \left(\mathbf{r}\right)=(2n+1)\pi /p$, with $n=0,1,\cdots,p-1$. These form the p-state clock models which correspond to the Ising ferromagnet for p  =  2 and to the 3-state Potts model for p  =  3.

For $p\geqslant 5$, there are many possibilities since higher harmonics of the site anisotropy are possible and also of the interaction terms. The most general Hamiltonian on a lattice with sites $\mathbf{r}$ with p-fold symmetry breaking is

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\underset{n=1}{\overset{[\,p/2]}{\mathop \sum}}\,\left[\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,{{K}_{n}}\left[1-\text{cos}\left[n\left(\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)\right)\right]\right]\right. \\ &&-\underset{\mathbf{r}}{\mathop \sum}\,{{h}_{np}}\text{cos}\left(np\theta \left(\mathbf{r}\right)\right].\end{eqnarray} \tag{ 93 }$$

This Hamiltonian represents many statistical mechanical models in two dimensions such as Potts models [105], the Ashkin–Teller model [106], the 8-vertex model [107], etc. These can all be written as generalized Coulomb gases [108–114]. For example, the 4-state Potts model is obtained when p  =  4, ${{h}_{4}}\to \infty$, h₈  =  0, so that $\theta \left(\mathbf{r}\right)=0,\pi /2,\pi,3\pi /2,$ and ${{K}_{1}}=2{{K}_{2}}$. Then the Hamiltonian becomes $H/\left({{k}_{\text{B}}}T\right)=-{{K}_{2}}\left(2{{\delta}_{s{{s}^{\prime}}}}+1\right)$ with $s=0,1\cdots,3$.

In the simplest case of n  =  1 only in (93), the symmetry breaking term in the partition function can be written as [55, 56]

$$\begin{eqnarray}&&{{\text{e}}^{{{h}_{p}}\text{cos}\,p\theta}}\to \underset{m=-\infty}{\overset{+\infty}{\mathop \sum}}\,{{\text{e}}^{\text{i}pm\theta \,}}y_{p}^{{{m}^{2}}},\end{eqnarray} \tag{ 94 }$$

with ${{y}_{p}}\approx {{h}_{p}}/2$. The Coulomb gas representation of the planar rotor in a p-fold symmetry breaking field becomes [56]

$$\begin{eqnarray}&&Z\propto {{Z}_{\text{c}}}\left(K,{{y}_{0}},{{y}_{p}}\right)=\underset{n\left(\mathbf{r}\right)}{\mathop \sum}\,\underset{m\left(\mathbf{R}\right)}{\mathop \sum}\,{{\text{e}}^{{{\mathcal{H}}_{\text{c}}}}},\\ &&{{\mathcal{H}}_{\text{c}}}=\pi {{K}_{0}}(T)\underset{\mathbf{R},{{\mathbf{R}}^{\prime}}}{\mathop \sum}\,n\left(\mathbf{R}\right)n\left({{\mathbf{R}}^{\prime}}\right)\text{ln}\left(\frac{|\mathbf{R}-{{\mathbf{R}}^{\prime}}|}{a}\right)+\text{ln}{{y}_{0}}\underset{\mathbf{R}}{\mathop \sum}\,{{n}^{2}}\left(\mathbf{R}\right)\\ &&+\frac{{{p}^{2}}}{4\pi {{K}_{0}}}\underset{\mathbf{r},{{\mathbf{r}}^{\prime}}}{\mathop \sum}\,m\left(\mathbf{r}\right)m\left({{\mathbf{r}}^{\prime}}\right)\text{ln}\left(\frac{|\mathbf{r}-{{\mathbf{r}}^{\prime}}|}{a}\right)+\text{ln}{{y}_{p}}\underset{\mathbf{r}}{\mathop \sum}\,{{m}^{2}}\left(\mathbf{r}\right)\\ &&+\text{i}p\underset{\mathbf{R},\mathbf{r}}{\mathop \sum}\,n\left(\mathbf{R}\right)m\left(\mathbf{r}\right)\Theta\left(\mathbf{R},\mathbf{r}\right),\\ &&\Theta\left(\mathbf{R},\mathbf{r}\right)=\text{ta}{{\text{n}}^{-1}}\left(\frac{Y-y}{X-x}\right).\end{eqnarray} \tag{ 95 }$$

The integer vortex charges $n\left(\mathbf{R}\right)$ lie on the sites $\mathbf{R}=(X,Y)$ of the dual lattice and the symmetry breaking charges $m\left(\mathbf{r}\right)$ lie on the original lattice with sites $\mathbf{r}=(x,y)$. Neutrality conditions $\sum_{\mathbf{R}}n\left(\mathbf{R}\right)=0=\sum_{\mathbf{r}}m\left(\mathbf{r}\right)$ are satisfied by both types of charges. The two Coulomb gases are coupled by the angular term $ip\sum_{\mathbf{R},\mathbf{r}}n\left(\mathbf{R}\right)m\left(\mathbf{r}\right)\Theta\left(\mathbf{R},\mathbf{r}\right)$. One can see immediately that equation (95) is self dual by relabeling $n\left(\mathbf{R}\right)\leftrightarrow m\left(\mathbf{r}\right)$, ${{y}_{0}}\leftrightarrow {{y}_{p}}$ and ${{K}_{0}}\leftrightarrow {{p}^{2}}/\left(4{{\pi}^{2}}{{K}_{0}}\right)$. The partition function is unchanged by these transformations [56], so that

$$\begin{eqnarray}&&Z\left(K,{{y}_{0}},{{y}_{p}}\right)=Z\left(\frac{{{p}^{2}}}{4{{\pi}^{2}}{{K}_{0}}},{{y}_{p}},{{y}_{0}}\right).\end{eqnarray} \tag{ 96 }$$

One can write down immediately the RG flow equations by rescaling the short distance cut off $a\to a{{\text{e}}^{l}}$ as expansions in the fugacities y₀ and y_p to lowest order [56]. These are consistent with the duality relation of equation (96)

$$\begin{eqnarray}&&\frac{\text{d}K_{0}^{-1}(l)}{\text{d}l}=4{{\pi}^{3}}y_{0}^{2}(l)-\pi {{p}^{2}}K_{0}^{-2}(l)y_{p}^{2}(l),\\ &&\frac{\text{d}{{y}_{0}}(l)}{\text{d}l}=\left(2-\pi {{K}_{0}}(l)\right){{y}_{0}}(l),\\ &&\frac{\text{d}{{y}_{p}}(l)}{\text{d}l}=\left(2-\frac{{{p}^{2}}}{4\pi {{K}_{0}}(l)}\right){{y}_{p}}(l).\end{eqnarray} \tag{ 97 }$$

From these equations, when p  >  4, there is a region of temperature when both fugacities y₀ and y_p are irrelevant given by

$$\begin{eqnarray}&&\frac{2}{\pi}\leqslant {{K}_{\text{R}}}(T)\leqslant \frac{{{p}^{2}}}{8\pi}.\end{eqnarray} \tag{ 98 }$$

In this range of T, the system is in a massless Gaussian phase. At low temperature $K_{0}^{-1}<8\pi /{{p}^{2}}$, y_p(l) is relevant and the system is in the low T phase where there is long range order with the system in one of the p ordered states. The phase diagram is shown in figure 14 left for p  =  6 and for p  =  4 in figure 14 right. When $p\leqslant 3$, there is no intermediate massless phase as one or both fugacities ${{y}_{0}},{{y}_{p}}$ are relevant for all K. This is interpreted as a direct transition from the large K ordered state to a low K disordered state. For p  =  2 and p  =  3, at ${{h}_{p}}=\infty$ it is easy to see that the Hamiltonian of (93) reduces to that of the Ising or the 3-state Potts model

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{K}{p-1}\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,\left(1-{{\delta}_{s\left(\mathbf{r}\right),s\left({{\mathbf{r}}^{\prime}}\right)}}\right),\,\,\,\,s\left(\mathbf{r}\right)=1,\cdots,p.\end{eqnarray} \tag{ 99 }$$

Zoom In Zoom Out Reset image size

I have shown that the planar rotor model, Z_p models and the discrete Gaussian model can be transformed into Coulomb gases which is useful for analyzing critical behavior. The transformation can be done exactly for Villain models [55, 108] because of the quadratic form of the Hamiltonian. However, many realistic models are not quadratic. For example, the solid on solid model in which the step energy is linear in $|h\left(\mathbf{r}\right)-h\left({{\mathbf{r}}^{\prime}}\right)|$ cannot be transformed simply into a Coulomb gas model. However, the Hamiltonian can be expanded in powers of $\nabla \theta$, the leading term being quadratic. One can then perform a more standard renormalization group analysis by momentum shell integration [75] and all non gaussian terms with more than two gradients are irrelevant. These RG arguments imply that SOS etc models with other than Villain interactions differ from a Coulomb gas by irrelevant operators only. Thus, provided one assumes that these models and Coulomb gases can be connected by RG flows with no fixed points intervening, they are asymptotically equivalent. Given this qualitative assumption, one can make quantitatively exact conclusions [83]. For the special case, p  =  4, the multicritical point, $\pi K=2$ and ${{y}_{0}}={{y}_{4}}=0$, is the meeting point of three critical lines of continuously varying exponents separating a high temperature disordered phase from two low T ordered phases with two different spin orientations with $\theta \left(\mathbf{r}\right)=n\pi /2$ or $\theta \left(\mathbf{r}\right)=(2n+1)\pi /4$, where n  =  0, 1, 2, 3. This point is identical to the critical point of the F model [115], and the three critical lines meeting there are in the universality class of the Baxter 8-vertex model [116].

Other interesting models of experimental relevance are obtained by including both the first and second harmonics of the p-fold anisotropy. The planar rotor model with both four and eight fold anisotropy

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=-K\underset{<\mathbf{r},{{\mathbf{r}}^{\prime}}>}{\mathop \sum}\,\text{cos}\left(\theta \left(\mathbf{r}\right)-\theta \left({{\mathbf{r}}^{\prime}}\right)\right)\\ &&-\underset{\mathbf{r}}{\mathop \sum}\,\left[{{h}_{4}}\text{cos}\left(4\theta \left(\mathbf{r}\right)\right)-{{h}_{8}}\text{cos}\left(8\theta \left(\mathbf{r}\right)\right)\right]\end{eqnarray} \tag{ 100 }$$

is an interesting theoretical model but can also describe the adsorption of hydrogen on the (1 0 0) surface of tungsten. When the hydrogen coverage is changed, h₈ does not change but h₄ can be driven through zero when a structural transition of the tungsten surface occurs [117, 118]. At low T, $\pi K\geqslant \pi /2$, equation (100) has two ordered phases with finite magnetization in one of the four directions $\theta \left(\mathbf{r}\right)=n\pi /2$ with n  =  0, 1, 2, 3 when h₄  >  0 and $\theta \left(\mathbf{r}\right)=(2n+1)\pi /4$ when h₄  <  0. The other anisotropy term ${{h}_{8}}\text{cos}8\theta \left(\mathbf{r}\right)$ is irrelevant for ${{K}^{-1}}>\pi /8$ and relevant for ${{K}^{-1}}\leqslant \pi /8$. In the temperature range $\pi /8\leqslant {{K}^{-1}}\leqslant \pi /2$ the ordered phases will be separated by a line of continuous transitions at h₄  =  0 At lower temperatures, ${{K}^{-1}}<\pi /8$, h₈ becomes relevant which has important consequences. There are two possibilities depending on the sign of h₈: (a) h₈  >  0 when the orderings favored by h₄ and by h₈ are compatible, independent of the sign of h₄. The transition between the two ordered phases at h₄  =  0 will be discontinuous first order [117, 118]. In case (b) h₈  <  0, the orderings favored by ${{h}_{4}},{{h}_{8}}$ are incompatible leading to an extra phase when $|{{h}_{4}}|$ is small with this phase separated from the other two ordered phases by continuous transitions in the 2D Ising universality class. Scenario (i) is compatible with experimental data [119, 120] which seems to favor a first order transition at low temperatures. Figure 15 illustrates schematic phase diagrams for h₈  >  0 (left) and h₈  <  0 (right).

Zoom In Zoom Out Reset image size

Other systems where the first two harmonics have important effects are tilted hexatic phases of liquid crystals. Experiments on the transitions between different tilted hexatic phases in thermotropic liquid crystals have observed a weakly first order transition from a hexatic-I to a hexatic-F phase [121]. In an analogous layered lyotropic liquid crystal the ${{L}_{\beta I}}$ and ${{L}_{\beta F}}$ phases, which correspond to the hexatic-I and hexatic-F phases respectively, were seen to be separated by a third ${{L}_{\beta L}}$ phase with continuous Ising-like transitions between the phases as in figure 15(b) [122]. These observations can be explained by modeling the system in terms of coupled planar rotor models for the bond angle, $\theta \left(\mathbf{r}\right)$, and the tilt azimuthal angle, $\phi \left(\mathbf{r}\right)$, of the director with Hamiltonian [123–125]

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\left(\frac{1}{2}{{K}_{6}}{{(\nabla \theta )}^{2}}+\frac{1}{2}{{K}_{1}}{{(\nabla \phi )}^{2}}+g\nabla \theta \centerdot \nabla \phi +V(\theta -\phi )\right),\\ &&V(\theta -\phi )=-{{h}_{6}}\text{cos}\left(6(\theta -\phi )\right)-{{h}_{12}}\text{cos}\left(12(\theta -\phi )\right)+\cdots,\end{eqnarray} \tag{ 101 }$$

which is very similar to (100). In the very low temperature regime of interest, but above the crystallization temperature ${{T}_{\text{m}}}$, one can write (101) as

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\left(\frac{1}{2}{{K}_{+}}{{\left(\nabla {{\theta}_{+}}\right)}^{2}}+\frac{1}{2}{{K}_{-}}{{\left(\nabla {{\theta}_{-}}\right)}^{2}}+V\left({{\theta}_{-}}\right)\right),\\ &&{{\theta}_{+}}=\alpha \theta +\beta \phi,\,\,\,\alpha =1-\beta =\frac{{{K}_{6}}+g}{{{K}_{6}}+{{K}_{1}}+2g},\\ &&{{K}_{+}}={{K}_{6}}+{{K}_{1}}+2g;\,\,\,{{K}_{-}}=\frac{{{K}_{1}}{{K}_{6}}-{{g}^{2}}}{{{K}_{+}}}.\end{eqnarray} \tag{ 102 }$$

A detailed analysis of this model system is carried out [125] with results agreeing with the experimental observations by adjusting the sign of h₁₂ in (101). The phase diagram for this model is essentially identical to figure 15 with ${{h}_{4}}\to {{h}_{6}}$ and ${{h}_{8}}\to {{h}_{12}}$ [125].

Dierker, Pindak and Meyer found a remarkable five-armed star defect sketched in figure 16 in thin tilted hexatic liquid crystal films at low temperatures [126]. This can be described by the Hamiltonian of (102) with h₁₂  >  0 so that there is a first order transition at h₆  =  0 between a hexatic I phase for h₆  >  0 and a hexatic F phase for h₆  <  0 [125]. It is easy to see that the Hamiltonian is minimized by ${{\theta}_{-}}\left(\mathbf{r}\right)=2\pi n/6$, $n=0,1,\cdots,5$, so that ${{\theta}_{-}}\left(\mathbf{r}\right)$ can have a discontinuity of $2\pi /6$. When a film is formed, it can have a vortex about which the tilt azimuthal angle $\phi \left(\mathbf{r}\right)$ goes through $2\pi$. When the system is cooled into the hexatic I phase, the bond orientational order increases and $\theta \left(\mathbf{r}\right)$ becomes locked to $\phi \left(\mathbf{r}\right)$ by the relevant h₆ term, so that the $2\pi$ vortex in $\phi \left(\mathbf{r}\right)$ becomes a $2\pi$ disclination in $\theta \left(\mathbf{r}\right)$. To reduce its energy, this defect breaks up into an N armed star structure with energy

$$\begin{eqnarray}&&\frac{\Delta E}{{{k}_{\text{B}}}T}=\frac{{{K}_{+}}N}{2}\int_{-\pi /N}^{\pi /N}\text{d} \Theta\int_{a}^{\infty}r\text{d}r\left(|\nabla \theta _{+}^{\text{star}}{{|}^{2}}-|\nabla \theta _{+}^{\text{no}\,\text{star}}{{|}^{2}}\right)+\epsilon RN,\end{eqnarray} \tag{ 103 }$$

where R is the length of an arm, ε is the energy per unit length in units of ${{k}_{\text{B}}}T$. As shown in [125, 126], $\Delta E/\left({{k}_{\text{B}}}T\right)$ is obtained by solving Laplace's equation ${{\nabla}^{2}}{{\theta}_{+}}(r,\Theta)=0$ in the wedge $-\pi /N<\Theta<+\pi /N$ subject to the boundary conditions

$$\begin{eqnarray}{{\theta}_{+}}(r,\pm \pi /N)=\left\{\begin{array}{c} \pm (\pi /N)(1-\alpha N/6) & \,\,\text{for}\,\,r<R \\ \pm \pi /N & \,\,\text{for}\,\,r>R \end{array}\right. \end{eqnarray} \tag{ 104 }$$

so that

$$\begin{eqnarray}&&\frac{\Delta E}{{{k}_{\text{B}}}T}=-N{{K}_{+}}g(N)\text{ln}\frac{R}{a}+N\epsilon \frac{R}{a}\,\,\text{where}\,\,g(N) \\ &&=\frac{\alpha \pi}{36}(12-\alpha -\alpha N).\end{eqnarray} \tag{ 105 }$$

Zoom In Zoom Out Reset image size

For the equilibrium arm length, minimize this with respect to R so that [125]

$$\begin{eqnarray}&&\frac{R}{a}=\frac{{{K}_{+}}g(N)}{\epsilon},\\ &&\frac{\Delta E}{{{k}_{\text{B}}}T}=-Ng(N){{K}_{+}}\left(\text{ln}\frac{{{K}_{+}}g(N)}{\epsilon}-1\right).\end{eqnarray} \tag{ 106 }$$

For reasonable physical values ${{K}_{6}}\gg {{K}_{1}},g$, so that the parameters of (101) become $\alpha +\beta =1$ with $\alpha \gg \beta$ and $\gamma \sim {{10}^{3}}$. For these parameters, the value N  =  5 minimizes the energy $\Delta E/{{k}_{\text{B}}}T$, in agreement with observation [126]. N  =  6 would relieve completely the bond orientational strain between the arms but the cost of the extra arm is prohibitive [126].

The expression for the arm length R can be used to test experimentally one prediction of the KTHNY theory of 2D melting [125, 126]. The equilibrium arm length of the arms of the star defect is given by (106) and is proportional to the temperature dependent coupling constant ${{K}_{+}}(T)={{K}_{6}}+{{K}_{1}}+2g$ where K₆ is the bond orientational stiffness called $K_{\text{A}}^{\text{R}}(T)$ in (153). In the liquid crystal system of [126], the coupling K₊ (T) is dominated by the bond orientational stiffness ${{K}_{\text{A}}}(T)$ whose temperature dependence is sketched in figure 19. Although the experiment of [126] is probably not in the asymptotic regime, the data for the arm length R(T) close to the crystal-hexatic transition was fitted to the form [126]

$$\begin{eqnarray}&&R(T)\sim \text{exp}\left(b{{\left(T-{{T}_{\text{m}}}\right)}^{-\tilde{\nu}}}\right).\end{eqnarray} \tag{ 107 }$$

This yields $\tilde{\nu}=0.38\pm 0.15$, which is consistent with the theoretical prediction $\tilde{\nu}=0.3696\cdots$ [127–130]. This is one of the earliest experiments to test the KTHNY theory of melting and the observations are completely consistent with all the theoretical predictions.

There have been a number of investigations of freely suspended thin liquid crystal films [100, 131–136] to look for the BKTHNY melting scenario and the elusive hexatic phase. This requires a material with the sequence of phases as T is lowered smectic-A—hexatic-B—crystal-B to avoid the extra complications due to tilted molecules. One of the aims was to investigate the predicted scaling relation between the harmonics of the bond orientational order parameter [137–139]

$$\begin{eqnarray}&&{{C}_{6n}}=\text{Re}\langle {{\text{e}}^{6\text{i}n\theta}}\rangle =C_{6}^{{{\sigma}_{n}}}\\ &&{{\sigma}_{n}}=n+\lambda n(n-1).\end{eqnarray} \tag{ 108 }$$

This has been verified in 3D [137] with $\lambda \approx 0.3$ and $\lambda =0.8$ for a thin smectic-C film [133, 134]. However, the theoretical value in 2D is $\lambda =1$ which follows from a low T Gaussian theory for the planar rotor model [138]. Finally, agreement was found between theory and experiment on a two layer film of smectic B 54COOBC which has bond but no herringbone order [139–141]. Interestingly, early measurements claimed disagreement with KT theory because of a specific heat anomaly at the Sm-A—Hex-B transition but later work [141] found a new Sm-A' phase with a different density and that the specific heat anomaly is at the Sm-A—Sm-A' transition while the Sm-A'—Hex-B—Cry-B phase transitions are consistent with the disclination/dislocation scenario.

Dislocations [148] are assumed to be responsible for the melting of a 2D crystal [18, 19] and, in analogy with vortices in the 2D planar rotor model, we want a correlation function which gives the renormalized elastic constants. Nelson and Halperin [128, 129] argued that inverse of the tensor of renormalized elastic constants $C_{ijkl}^{\text{R}}$ is

$$\begin{eqnarray}&&{{C}_{\text{R};ijkl}}={{\tilde{\mathop{\mu}}\,}_{\text{R}}}(T)\left({{\delta}_{ik}}{{\delta}_{jl}}+{{\delta}_{il}}{{\delta}_{jk}}\right)+{{\tilde{\mathop{\lambda}}\,}_{\text{R}}}(T){{\delta}_{ij}}{{\delta}_{kl}},\\ &&C_{\text{R};ijkl}^{-1}=\frac{\langle {{U}_{ij}}{{U}_{kl}}\rangle}{\Omega a_{0}^{2}},\\ &&{{U}_{ij}}=-\frac{1}{2}{{\mathop{\oint}^{}}_{\text{P}}}\text{d}l\left({{u}_{i}}{{\hat{n}}_{j}}+{{u}_{j}}{{\hat{n}}_{i}}\right),\end{eqnarray} \tag{ 124 }$$

where P is the perimeter of the crystal and $\mathbf{\hat{n}}$ is a unit vector normal to the perimeter pointing outwards from the crystal. Note that this definition of the renormalized compliance tensor or inverse elastic tensor, $C_{\text{R};ijkl}^{-1}$ in terms of fluctuations at the perimeter is exactly how the compliances are measured. Applying Green's theorem to the expression for U_ij in (124), one obtains [129]

$$\begin{eqnarray}&&{{U}_{ij}}=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{u}_{ij}}\left(\mathbf{r}\right)+\frac{1}{2}\underset{\mathbf{R}}{\mathop \sum}\,\left[{{b}_{i}}\left(\mathbf{R}\right){{\epsilon}_{jl}}{{R}_{l}}+{{b}_{j}}\left(\mathbf{R}\right){{\epsilon}_{il}}{{R}_{l}}\right],\end{eqnarray} \tag{ 125 }$$

where $\mathbf{r}$ denote the sites of the triangular crystal lattice and $\mathbf{R}$ the sites of the hexagonal dual lattice on which live the dislocations. It is not hard to see that the dislocation contribution to $\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{u}_{ij}}\left(\mathbf{r}\right)$ vanishes, so that the singular part of U_ij is

$$\begin{eqnarray}&&{{U}_{ij}}=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\phi}_{ij}}\left(\mathbf{r}\right)+U_{ij}^{s},\\ &&U_{ij}^{s}=\frac{1}{2}{{a}_{0}}\underset{\mathbf{R}}{\mathop \sum}\,\left[{{b}_{i}}\left(\mathbf{R}\right){{\epsilon}_{jl}}{{R}_{l}}+{{b}_{j}}\left(\mathbf{R}\right){{\epsilon}_{il}}{{R}_{l}}\right].\end{eqnarray} \tag{ 126 }$$

We have derived a suitable correlation function for the renormalized elastic constants, or response functions, from which we can deduce renormalization group flows, in exact analogy to the renormalized stiffness constant ${{K}_{\text{R}}}(T)={{\hslash}^{2}}\rho _{s}^{\text{R}}(T)/\left({{m}^{2}}{{k}_{\text{B}}}T\right)$ of equation (57). We have

$$\begin{eqnarray}&&C_{\text{R};ijkl}^{-1}=C_{ijkl}^{-1}+{{\Omega }^{-1}}\langle U_{ij}^{s}U_{kl}^{s}\rangle,\\ &&C_{ijkl}^{-1}=\frac{1}{\Omega a_{0}^{2}}\langle \mathop{\int}^{}{{\text{d}}^{2}}{{\mathbf{r}}_{1}}{{\phi}_{ij}}\left({{\mathbf{r}}_{1}}\right)\mathop{\int}^{}{{\text{d}}^{2}}{{\mathbf{r}}_{2}}{{\phi}_{kl}}\left({{\mathbf{r}}_{2}}\right)\rangle \\ &&=\frac{1}{4\tilde{\mathop{\mu}}\,}\left({{\delta}_{ik}}{{\delta}_{jl}}+{{\delta}_{il}}{{\delta}_{jk}}\right)-\frac{\tilde{\mathop{\lambda}}\,}{4\tilde{\mathop{\mu}}\,\left(\tilde{\mathop{\mu}}\,+\tilde{\mathop{\lambda}}\,\right)}{{\delta}_{ij}}{{\delta}_{kl}}.\end{eqnarray} \tag{ 127 }$$

The correlation function $\langle U_{ij}^{s}U_{kl}^{s}\rangle$ is evaluated by a perturbation expansion in the dislocation fugacity $y={{\text{e}}^{-{{E}_{\text{c}}}/{{k}_{\text{B}}}T}}$ to obtain the RG flow equations [128–130]

$$\begin{eqnarray}&&\frac{\text{d}{{\tilde{\mathop{\mu}}\,}^{-1}}(l)}{\text{d}l}=3\pi {{y}^{2}}(l){{\text{e}}^{\tilde{\mathop{K}}\,(l)/8\pi}}{{I}_{0}}\left(\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{d{{\left(\tilde{\mathop{\mu}}\,(l)+\tilde{\mathop{\lambda}}\,(l)\right)}^{-1}}}{\text{d}l}=3\pi {{y}^{2}}(l){{\text{e}}^{\tilde{\mathop{K}}\,(l)/8\pi}}\left[{{I}_{0}}\left(\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)-{{I}_{1}}\left(\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)\right]+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{\text{d}{{\tilde{\mathop{K}}\,}^{-1}}(l)}{\text{d}l}=\frac{3}{4}\pi {{y}^{2}}(l){{\text{e}}^{\tilde{\mathop{K}}\,(l)/8\pi}}\left(2{{I}_{0}}\left(\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)-{{I}_{1}}\left(\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)\right)+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{\text{d}y(l)}{\text{d}l}=\left(2-\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)y(l)+2\pi {{y}^{2}}(l){{\text{e}}^{\tilde{\mathop{K}}\,(l)/8\pi}}{{I}_{0}}\left(\frac{\tilde{\mathop{K}}\,(l)}{8\pi}\right)+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\tilde{\mathop{K}}\,(l)=\frac{4\tilde{\mathop{\mu}}\,(l)\left(\tilde{\mathop{\mu}}\,(l)+\tilde{\mathop{\lambda}}\,(l)\right)}{2\tilde{\mathop{\mu}}\,(l)+\tilde{\mathop{\lambda}}\,(l)}.\end{eqnarray} \tag{ 128 }$$

Here, $\tilde{K}(l)$ is the coupling constant in the dislocation Hamiltonian ${{H}_{\text{D}}}/\left({{k}_{\text{B}}}T\right)$ of (121). The derivation of the flow equations (128) is identical in spirit to the vortex flows of (37) but more complex because there are three different elementary dislocations in a triangular lattice and only one elementary vortex in the planar rotor model. The origin of the term $\mathcal{O}\left(\,{{y}^{2}}\right)$ in the flow equation for the dislocation fugacity y(l) is that two different elementary dislocations can combine to form a third elementary dislocation [128–130]. The RG flows for ${{\tilde{K}}^{-1}}(l)$ and y(l) are sketched in figure 18.

Zoom In Zoom Out Reset image size

The recursion relations are very similar to those for the planar rotor model. Above a temperature ${{T}_{\text{m}}}$, $\tilde{K}(l)<16\pi$ and the fugacity y(l) increases with l, implying that the dislocations unbind and become free to move under any small applied stress so that the elastic moduli vanish and the crystal has melted. The flow equations (128) up to a length scale ${{\xi}_{+}}(T)$ [129],

$$\begin{eqnarray}&&{{\xi}_{+}}(T)={{\text{e}}^{{{l}^{\ast}}}}=\text{exp}\left({{b}^{\prime}}{{t}^{-\nu}}\right)\,\,\,\,\,\text{with}\,\,\,\,\nu =0.3696\cdots.\end{eqnarray} \tag{ 129 }$$

with ${{b}^{\prime}}$ a non universal constant. The scale ${{\xi}_{+}}(T)$ can be interpreted as the scale below which dislocations are bound together in pairs or triplets of zero total Burger's vector. Equivalently, on length scales ${{\text{e}}^{{{l}^{\ast}}}}\leqslant {{\xi}_{+}}(T)$, the system responds elastically, while on length scales ${{\text{e}}^{{{l}^{\ast}}}}>{{\xi}_{+}}(T)$, the system responds like a fluid.

When $T\leqslant {{T}_{\text{m}}}$, $\tilde{K}(l)\geqslant 16\pi$ and the fugacity $y(l)\to 0$ as $l\to \infty$ so that, at all length scales, the system responds elastically. The renormalization group transformation is derived from the perturbation expansion for the renormalized elastic constants ${{\tilde{\mu}}_{\text{R}}}(T)$, ${{\tilde{\lambda}}_{\text{R}}}(T)$ and ${{\tilde{K}}_{\text{R}}}(T)$ so that [129]

$$\begin{eqnarray}&&{{\tilde{\mathop{\mu}}\,}_{\text{R}}}\left(\tilde{\mathop{\mu}}\,,\tilde{\mathop{\lambda}}\,,y\right)={{\tilde{\mathop{\mu}}\,}_{\text{R}}}\left(\tilde{\mathop{\mu}}\,(l),\tilde{\mathop{\lambda}}\,(l),y(l)\right),\\ &&{{\tilde{\mathop{\lambda}}\,}_{\text{R}}}\left(\tilde{\mathop{\mu}}\,,\tilde{\mathop{\lambda}}\,,y\right)={{\tilde{\mathop{\lambda}}\,}_{\text{R}}}\left(\tilde{\mathop{\mu}}\,(l),\tilde{\mathop{\lambda}}\,(l),y(l)\right),\\ &&{{\tilde{\mathop{K}}\,}_{\text{R}}}\left(\tilde{\mathop{K}}\,,y\right)={{\tilde{\mathop{K}}\,}_{\text{R}}}\left(\tilde{\mathop{K}}\,(l),y(l)\right).\end{eqnarray} \tag{ 130 }$$

The Young's modulus, ${{\tilde{K}}_{\text{R}}}(T)$ is analogous to the renormalized stiffness ${{K}_{\text{R}}}(T,y)$ for a ⁴He film. At low T, corresponding to $\tilde{K}(l)\geqslant 16\pi$, from (128), it is seen that the fugacity y(l) is irrelevant so that the system has no free dislocations at large length scales and the stiffness

$$\begin{eqnarray}&&{{\tilde{\mathop{\mu}}\,}_{\text{R}}}\left(\tilde{\mathop{\mu}}\,,\tilde{\mathop{\lambda}}\,,y\right)=\underset{l\to \infty}{\mathop{\lim}}\,\tilde{\mathop{\mu}}\,(l),\\ &&{{\tilde{\mathop{\lambda}}\,}_{\text{R}}}\left(\tilde{\mathop{\mu}}\,,\tilde{\mathop{\lambda}}\,,y\right)=\underset{l\to \infty}{\mathop{\lim}}\,\tilde{\mathop{\lambda}}\,(l),\\ &&{{\tilde{\mathop{K}}\,}_{\text{R}}}(K,y)=\underset{l\to \infty}{\mathop{\lim}}\,\tilde{\mathop{K}}\,(l).\end{eqnarray} \tag{ 131 }$$

For $T\leqslant {{T}_{\text{m}}}$, the values of ${{\tilde{\mu}}_{\text{R}}}(T)$ and ${{\tilde{\lambda}}_{\text{R}}}(T)$ can, in principle, be obtained by integrating (128) and they depend on the initial values ${{\mu}_{0}},{{\lambda}_{0}}$ and are non universal. However, the combination ${{\tilde{K}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)=16\pi$ does have a universal value from (128), just like the universal stiffness of the 2D planar rotor model [129],

$$\begin{eqnarray}&&{{\tilde{K}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)=\underset{T\to T_{\text{m}}^{-}}{\mathop{\lim}}\,\frac{4{{{\tilde{\mu}}}_{\text{R}}}(T)\left({{{\tilde{\mu}}}_{\text{R}}}(T)+{{{\tilde{\lambda}}}_{\text{R}}}(T)\right)}{2{{{\tilde{\mu}}}_{\text{R}}}(T)+{{{\tilde{\lambda}}}_{\text{R}}}(T)}=16\pi.\end{eqnarray} \tag{ 132 }$$

To obtain more detailed behavior near the transition at $16\pi {{K}^{-1}}=1$, assume the dislocation fugacity $y\ll 1$, write $16\pi {{K}^{-1}}(l)=1+x(l)$ and expand the flow equations (128) to lowest order in x(l) to obtain [128–130]

$$\begin{eqnarray}&&\frac{\text{d}x(l)}{\text{d}l}=12{{\pi}^{2}}A{{y}^{2}}(l)+\mathcal{O}\left(\,{{y}^{3}}\right)={{Y}^{2}}(l)+\mathcal{O}\left({{Y}^{3}}\right),\,\,\,\,\\ &&Y(l)=\pi \sqrt{12A}y(l)\\ &&\frac{\text{d}Y(l)}{\text{d}l}=2x(l)Y(l)+2\alpha {{Y}^{2}}(l)+\mathcal{O}\left({{Y}^{3}}\right),\\ &&A={{e}^{2}}\left(2{{I}_{0}}(2)-{{I}_{1}}(2)\right)\approx 21.94,\\ &&B={{e}^{1}}{{I}_{0}}(2)\approx 6.20,\\ &&\alpha =\frac{B}{\sqrt{12A}}.\end{eqnarray} \tag{ 133 }$$

The flows in the x, Y plane are obtained from

$$\begin{eqnarray}&&\frac{\text{d}Y}{\text{d}x}=2\frac{x}{Y}+2\alpha,\end{eqnarray} \tag{ 134 }$$

which yields the renormalization group invariant

$$\begin{eqnarray}&&|Y(l)-\sqrt{{{\alpha}^{2}}+2}(1+b)x(l){{|}^{(1+b)/2}}|Y(l)+\sqrt{{{\alpha}^{2}}+2}(1-b)x(l){{|}^{(1-b)/2}}=C(t),\\ &&b=\frac{\alpha}{\sqrt{{{\alpha}^{2}}+2}},\\ &&C(t)=|{{Y}_{0}}-\sqrt{{{\alpha}^{2}}+2}(1+b){{X}_{0}}{{|}^{(1+b)/2}}|{{Y}_{0}}+\sqrt{{{\alpha}^{2}}+2}(1-b){{X}_{0}}{{|}^{(1-b)/2}}.\end{eqnarray} \tag{ 135 }$$

Note that, when the parameter B  =  0, this reduces to the planar rotor problem since then (135) is just the square root of (49). The RG flows from (135) are shown in figure 18.

One can extract some more detailed behavior of the system in the vicinity of the melting transition at $T={{T}_{\text{m}}}$ by solving the recursion relations of equation (133) in detail [127, 129, 130]. To simplify the algebra, it is convenient to redefine the fugacity as $Y=\pi \sqrt{12A}y$ so that the recursion relations for x and Y become [130]

$$\begin{eqnarray}&&\frac{\text{d}x}{\text{d}l}={{Y}^{2}},\\ &&\frac{\text{d}Y}{\text{d}l}=2xY+2\alpha {{Y}^{2}},\,\,\,\,\,\alpha =\frac{B}{\sqrt{12A}}.\end{eqnarray} \tag{ 136 }$$

On the critical line, Y(l)  =  −m₀x(l) with ${{m}_{0}}=\sqrt{{{\alpha}^{2}}+2}-\alpha$ and

$$\begin{eqnarray}&&\frac{\text{d}|x(l)|}{\text{d}l}=-m_{0}^{2}|x(l){{|}^{2}},\,\,\,\,\Rightarrow \,\,\,\,x(l)=-\frac{|{{x}_{0}}|}{1+|{{x}_{0}}|m_{0}^{2}l}.\end{eqnarray} \tag{ 137 }$$

For x₀ slightly away from the critical separatrix, write [129, 130] Y(l)  =  −m₀x(l)  +  D(l) where

$$\begin{eqnarray}&&\frac{\text{d}D(l)}{\text{d}l}=-2x(l)D(l)+\mathcal{O}\left({{D}^{2}}\right)\,\,\,\Rightarrow \,\,\,D(l)={{D}_{0}}{{\left(1\,+\,|{{x}_{0}}|m_{0}^{2}l\right)}^{2/m_{0}^{2}}}.\end{eqnarray} \tag{ 138 }$$

The initial value ${{D}_{0}}\propto t=\left(T-{{T}_{\text{m}}}\right)/{{T}_{\text{m}}}\ll 1$. When t  <  0, the flow follows the critical separatrix for some distance and then Y(l) breaks away and plunges rapidly to for l  >  l^*. One can calculate l^* from the condition $D\left({{l}^{\ast}}\right)\approx Y\left({{l}^{\ast}}\right)$ which gives

$$\begin{eqnarray}&&{{D}_{0}}{{\left(1\,+\,|{{x}_{0}}|m_{0}^{2}{{l}^{\ast}}\right)}^{2/m_{0}^{2}}}\approx \frac{{{m}_{0}}|{{x}_{0}}|}{1+|{{x}_{0}}|m_{0}^{2}{{l}^{\ast}}},\\ &&{{l}^{\ast}}\sim D_{0}^{-m_{0}^{2}/\left(m_{0}^{2}+2\right)}\sim |t{{|}^{-\tilde{\mathop{\nu}}\,}}.\end{eqnarray} \tag{ 139 }$$

From this we can define the length scale ${{\xi}_{-}}(T)$ as the scale at which the deviation from the critical separatrix becomes significant [128–130]

$$\begin{eqnarray}&&{{\xi}_{-}}(T)={{\text{e}}^{{{l}^{\ast}}}}={{\text{e}}^{b|t{{|}^{-\tilde{\nu}}}}}.\end{eqnarray} \tag{ 140 }$$

This length scale exists both above and below ${{T}_{\text{m}}}$ and has the same interpretation as for the planar rotor model. When $T\leqslant {{T}_{\text{m}}}$, ${{\xi}_{-}}(T)$ is the separation of the largest bound pair of primitive dislocations and there are no larger pairs. For $T>{{T}_{\text{m}}}$, up to this scale, the dislocations can be regarded as bound in neutral pairs and triplets while at larger length scales the dislocations are unbound and freely mobile.

Since the dislocation fugacity y(l) scales to zero for $T\leqslant {{T}_{\text{m}}}$, the system is elastic and the renormalized elastic constants ${{\tilde{\mu}}_{\text{R}}}(T)$ and ${{\tilde{\lambda}}_{\text{R}}}(T)$ are finite but non universal and are given by equation (131). The Young's modulus ${{\tilde{K}}_{\text{R}}}(T)$ has a universal value at $T_{m}^{-}$ of ${{\tilde{K}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)=16\pi$ with [129]

$$\begin{eqnarray}&&{{\tilde{K}}_{\text{R}}}(T)=16\pi \left(1+c|t{{|}^{{\tilde{\nu}}}}\right),\\ &&{{\tilde{\mu}}_{\text{R}}}(T)={{\tilde{\mu}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)\left(1+c|t{{|}^{{\tilde{\nu}}}}\right),\\ &&{{\tilde{\lambda}}_{\text{R}}}(T)={{\tilde{\lambda}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)\left(1+c|t{{|}^{{\tilde{\nu}}}}\right)\end{eqnarray} \tag{ 141 }$$

where c is a constant $\mathcal{O}(1)$. This behavior has been verified in detail by experiment and simulation [144]. The x-ray structure function

$$\begin{eqnarray}&&S\left(\mathbf{q}\right)=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\text{e}}^{\text{i}\mathbf{q}\centerdot \mathbf{r}}}\langle {{\text{e}}^{\text{i}\mathbf{q}\centerdot \left[\mathbf{u}\left(\mathbf{r}\right)-\mathbf{u}(0)\right]}}\rangle \end{eqnarray} \tag{ 142 }$$

where the Debye Waller factor for $\mathbf{q}=\mathbf{G}$ is

$$\begin{eqnarray}&&{{C}_{\mathbf{G}}}\left(\mathbf{r}\right)=\langle {{\rho}_{\mathbf{G}}}\left(\mathbf{r}\right){{\rho}_{-\mathbf{G}}}(0)\rangle =\langle {{\text{e}}^{\text{i}\mathbf{G}\centerdot \left(\mathbf{u}\left(\mathbf{r}\right)-\mathbf{u}(0)\right)}}\rangle \sim {{r}^{-{{\eta}_{\mathbf{G}}}(T)}}\\ &&{{\eta}_{\mathbf{G}}}(T)=\frac{|\mathbf{G}{{|}^{2}}a_{0}^{2}\left(3{{{\tilde{\mu}}}_{\text{R}}}(T)+{{{\tilde{\lambda}}}_{\text{R}}}(T)\right)}{4\pi {{{\tilde{\mu}}}_{\text{R}}}(T)\left({{{\tilde{\mu}}}_{\text{R}}}(T)+{{{\tilde{\lambda}}}_{\text{R}}}(T)\right)}.\end{eqnarray} \tag{ 143 }$$

The structure function for $|\mathbf{q}-\mathbf{G}|\ll 1/a$ is the Fourier transform of the Debye–Waller factor ${{C}_{\mathbf{G}}}\left(\mathbf{r}\right)$ so that $S\left(\mathbf{q}\right)\sim |\mathbf{q}-\mathbf{G}{{|}^{{{\eta}_{\text{G}}}(T)-2}}$. This diverges for the smaller values of $\mathbf{G}$ when ${{\eta}_{\text{G}}}<2$ as sketched in figure 17. Note that for a hexagonal lattice with lattice spacing a₀, the smallest reciprocal lattice vector $|{{\mathbf{G}}_{0}}|\,=4\pi /\left({{a}_{0}}\sqrt{3}\right)$ so that

$$\begin{eqnarray}&&{{\eta}_{{{\mathbf{G}}_{0}}}}(T)=\frac{4\pi}{3}\frac{3{{{\tilde{\mu}}}_{\text{R}}}(T)+{{{\tilde{\lambda}}}_{\text{R}}}(T)}{{{{\tilde{\mu}}}_{\text{R}}}(T)\left(2{{{\tilde{\mu}}}_{\text{R}}}(T)+{{{\tilde{\lambda}}}_{\text{R}}}(T)\right)}.\end{eqnarray} \tag{ 144 }$$

For systems with long range repulsive interactions, such as electrons trapped above the surface of superfluid ⁴He [150] and colloidal particles at a water/air interface [144], the ratio ${{\tilde{\lambda}}_{\text{R}}}(T)/{{\tilde{\mu}}_{\text{R}}}(T)\gg 1$ so that

$$\begin{eqnarray}&&{{\tilde{\mu}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)=4\pi \,\,\,\,\Rightarrow \,\,\,\,{{\eta}_{{{\mathbf{G}}_{0}}}}\left(T_{\text{m}}^{-}\right)=\frac{4\pi}{3{{{\tilde{\mu}}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)}=\frac{1}{3}.\end{eqnarray} \tag{ 145 }$$

This is not a very good estimate for a detailed comparison with theory because, the initial value of ${{\lambda}_{0}}=\infty$ is renormalized by dislocations to a finite value ${{\tilde{\lambda}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)/{{\tilde{\mu}}_{\text{R}}}\left(T_{\text{m}}^{-}\right)\approx 10$ for colloidal particles interacting by a 1/r³ repulsive interaction [143] and is $\approx 23$ for electrons on ⁴He [151]. The quantity measured by experiment on the electron system is $\Gamma(T)=\sqrt{\pi {{n}_{s}}}{{e}^{2}}/{{k}_{\text{B}}}T=137\pm 15$ at $T={{T}_{\text{m}}}$ [150] and a remarkably good fit has been obtained by combining a molecular dynamics simulation with the renormalization group [151–154].

Above ${{T}_{\text{m}}}$, the dislocation fugacity is relevant and, just as in the planar rotor model, there is a correlation length ${{\xi}_{+}}(T)$ which is the length below which dislocations can be considered as bound and above which they are free [129, 130] in analogy with vortices in the planar rotor model [19, 20]. The positional correlation function will decay exponentially for $r\gg {{\xi}_{+}}$

$$\begin{eqnarray}&&{{C}_{\mathbf{G}}}\left(\mathbf{r}\right)\,\,\sim \,\,{{\text{e}}^{-r/{{\xi}_{+}}(T)}},\\ &&{{\xi}_{+}}(T)\sim {{\text{e}}^{b{{t}^{-\tilde{\nu}}}}},\,\,\,\text{where}\,\,b=\text{constant}>0.\end{eqnarray} \tag{ 146 }$$

The scale dependent elastic constants ${{\tilde{\mu}}_{\text{R}}}(l),\,{{\tilde{\lambda}}_{\text{R}}}(l)=0$ for ${{\text{e}}^{l}}>{{\xi}_{+}}(T)$ but are finite for ${{\text{e}}^{l}}<{{\xi}_{+}}(T)$. This is an example of a system which behaves as an elastic solid at short length scales ${{\text{e}}^{l}}<{{\xi}_{+}}$ and as a fluid for longer scales ${{\text{e}}^{l}}>{{\xi}_{+}}$. The x-ray structure function

$$\begin{eqnarray}{{S}_{\text{G}}}\left(\mathbf{q}\right)\equiv \mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\text{e}}^{-\text{i}\left(\mathbf{q}-\mathbf{G}\right)\centerdot \mathbf{r}}}{{C}_{\text{G}}}\left(\mathbf{r}\right)\sim \left\{\begin{array}{c} |\mathbf{q}-\mathbf{G}{{|}^{-2+{{\eta}_{\text{G}}}(T)}} & \left(T<{{T}_{\text{c}}}\right), \\ {{\left(|\mathbf{q}-\mathbf{G}{{|}^{2}}+\xi _{+}^{-2}\right)}^{-1}} & \left(T>{{T}_{\text{c}}}\right), \end{array}\right. \end{eqnarray} \tag{ 147 }$$

as sketched in figure 17.

At first sight, this phase where dislocations are all unbound and free can be interpreted as a conventional isotropic fluid [18, 19, 127] as the translational order decays exponentially with distance. However, in their seminal work, Halperin and Nelson [128, 129] pointed out that this fluid phase is not the expected isotropic fluid but has remnants of the six-fold orientational order of the underlying hexagonal lattice of the low temperature solid phase. These possible orientational correlations may be described by a finite renormalized Frank constant ${{K}_{\text{A}}}(T)$ of the anisotropic fluid

$$\begin{eqnarray}&&K_{\text{A}}^{-1}(T)=\underset{q\to 0}{\mathop{\lim}}\,\frac{{{q}^{2}}}{\Omega }\langle \theta \left(\mathbf{q}\right)\theta \left(-\mathbf{q}\right)\rangle,\end{eqnarray} \tag{ 148 }$$

where θ is the angle of a bond between nearest neighbor particles relative to an arbitrary reference direction.

One expects that the distribution of bond angles $\theta \left(\mathbf{r}\right)$ are controlled by a free energy [128, 129],

$$\begin{eqnarray}&&\frac{{{F}_{\text{A}}}}{{{k}_{\text{B}}}T}=\frac{1}{2}{{K}_{\text{A}}}(T)\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla \theta \left(\mathbf{r}\right)\right)}^{2}}\end{eqnarray} \tag{ 149 }$$

where ${{K}_{\text{A}}}(T)$ must be finite for the existence of orientational correlations. In the underlying solid phase, the bond angle ${{\theta}_{s}}\left(\mathbf{r}\right)$ due to a set of dislocations is

$$\begin{eqnarray}&&{{\theta}_{s}}\left(\mathbf{r}\right)=-\frac{{{a}_{0}}}{2\pi}\underset{{{\mathbf{r}}^{\prime}}}{\mathop \sum}\,\frac{\mathbf{b}\left({{\mathbf{r}}^{\prime}}\right)\centerdot \left(\mathbf{r}-{{\mathbf{r}}^{\prime}}\right)}{|\mathbf{r}-{{\mathbf{r}}^{\prime}}{{|}^{2}}},\end{eqnarray} \tag{ 150 }$$

from which one can write the stiffness constant as a correlation function,

$$\begin{eqnarray}&&K_{\text{A}}^{-1}(T)=\underset{q\to 0}{\mathop{\lim}}\,\frac{a_{0}^{2}}{\Omega }\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{2}}}\langle {{b}_{i}}\left(\mathbf{q}\right){{b}_{j}}\left(-\mathbf{q}\right)\rangle.\end{eqnarray} \tag{ 151 }$$

The dislocation correlation function is to be calculated from the dislocation free energy of equation (121) written in Fourier space as [128, 129]

$$\begin{eqnarray}&&\frac{{{H}_{\text{D}}}}{{{k}_{\text{B}}}T}=\frac{1}{2\Omega }\underset{q}{\mathop \sum}\,\left[\frac{{\tilde{K}}}{{{q}^{2}}}\left({{\delta}_{ij}}-\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{2}}}\right)+\frac{2{{E}_{\text{c}}}a_{0}^{2}}{{{k}_{\text{B}}}T}{{\delta}_{ij}}\right]{{b}_{i}}\left(\mathbf{q}\right){{b}_{j}}\left(-\mathbf{q}\right).\end{eqnarray} \tag{ 152 }$$

At this point, one sees why the sign of the angular term in the dislocation interaction in (121) is so important and why ignoring the apparently shorter ranged angular term leads to incorrect results [19, 127]. The existence of a finite Franck constant, ${{K}_{\text{A}}}(T)$, depends on the transverse projection operator, ${{T}_{ij}}(q)={{\delta}_{ij}}-{{q}_{i}}{{q}_{j}}/{{q}^{2}}$, in (152). One can understand the behavior of the orientational stiffness constant ${{K}_{\text{A}}}(T)$ by using the recursion relations (128) up to the limit of their validity, ${{l}^{\ast}}=\text{ln}{{\xi}_{+}}(T)$, to obtain

$$\begin{eqnarray}&&{{K}_{\text{A}}}(T)={{K}_{\text{A}}}\left(\tilde{\mu}(0),\tilde{\lambda}(0),{{y}_{0}}\right)={{\text{e}}^{2l}}{{K}_{\text{A}}}\left(\tilde{\mu}(l),\tilde{\lambda}(l),y(l)\right).\end{eqnarray} \tag{ 153 }$$

The scaling factor ${{\text{e}}^{2l}}$ arises from the two extra factors of q in the definition of ${{K}_{\text{A}}}$ in equation (151) compared to $\tilde{\mu}$ and $\tilde{\lambda}$. At length scale ${{\text{e}}^{{{l}^{\ast}}}}={{\xi}_{+}}(T)$, it is reasonable to assume that there is a large density of free dislocations in the system and a sensible approximation is to regard the dislocation density $\mathbf{b}\left(\mathbf{r}\right)$ as a continuous variable as all values of $\mathbf{b}\left(\mathbf{r}\right)$ are present at scale l^*. This Debye–Hückel approximation allows the correlation function $\langle {{b}_{i}}\left(\mathbf{q}\right){{b}_{j}}\left(-\mathbf{q}\right)\rangle$ in (151) to be evaluated because the transverse part of the dislocation interaction in equation (152) does not contribute. One obtains

$$\begin{eqnarray}&&K_{\text{A}}^{-1}(T)\approx \frac{{{k}_{\text{B}}}T}{2{{E}_{\text{c}}}a_{0}^{2}}>0\end{eqnarray} \tag{ 154 }$$

which establishes that the fluid caused by free dislocations above ${{T}_{\text{m}}}$ does have orientational stiffness [128, 129] and is $not$ isotropic as assumed by Kosterlitz and Thouless [19].

Now one is left with a problem of a system with short range exponentially decaying translational order and algebraic orientational order described by phenomenological free energy of equation (149) with a finite bare orientational stiffness constant ${{K}_{\text{A}}}(T)$ of (154). This is the planar rotor problem discussed earlier in detail, the only difference that the topological defects are $\pi /3$ disclinations which interact logarithmically like vortices. The bond angle field $\theta \left(\mathbf{r}\right)$ is decomposed into a smooth part $\phi \left(\mathbf{r}\right)$ and a singular part ${{\theta}_{s}}\left(\mathbf{r}\right)$ which obeys

$$\begin{eqnarray}&&{{\mathop{\oint}^{}}_{\text{C}}}\text{d}l{{\theta}_{s}}\left(\mathbf{r}\right)=\frac{m\pi}{3}\,\,\,m=0\pm 1,\pm 2,\cdots,\end{eqnarray} \tag{ 155 }$$

where m is the total disinclination strength enclosed by the contour C. Thus one obtains the disclination free energy [128, 129] which is identical to the planar rotor problem of equation (31)

$$\begin{eqnarray}&&\left. \frac{{{H}_{\text{D}}}}{{{k}_{\text{B}}}T}=\frac{1}{2}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\right){{\left(\nabla \phi \left(\mathbf{r}\right)\right)}^{2}}-\frac{\pi {{K}_{\text{A}}}}{36}\underset{\mathbf{r},{{\mathbf{r}}^{\prime}}}{\mathop \sum}\,m\left(\mathbf{r}\right)m\left({{\mathbf{r}}^{\prime}}\right)\text{ln}\left(\frac{|\mathbf{r}-{{\mathbf{r}}^{\prime}}|}{a}\right)+\frac{{{E}_{\text{c}}}}{{{k}_{\text{B}}}T}\underset{\mathbf{r}}{\mathop \sum}\,{{m}^{2}}\left(\mathbf{r}\right).\end{eqnarray} \tag{ 156 }$$

Here $m\left(\mathbf{r}\right)=0,\pm 1,\cdots$ is a measure of the strength of the disinclination at $\mathbf{r}$, a is the core size of a disinclination, ${{E}_{\text{c}}}$ is the disclination core energy and the disclinations are subject to the neutrality condition, $\sum_{\mathbf{r}}m\left(\mathbf{r}\right)=0$.

Halperin and Nelson [128, 129] pointed out that one can immediately take over all the results for the superfluid film and showed there is a disinclination unbinding transition at T  =  T_I to the high temperature isotropic fluid. The renormalized orientational stiffness constant just below T_I is given by

$$\begin{eqnarray}&&\frac{36}{2\pi K_{\text{A}}^{\text{R}}\left(T_{I}^{-}\right)}={{\eta}_{6}}\left(T_{I}^{-}\right)=\frac{1}{4}.\end{eqnarray} \tag{ 157 }$$

The renormalized Franck constant, $K_{\text{A}}^{\text{R}}(T)$, behaves precisely like the renormalized stiffness constant of a 2D superfluid film except that, slightly above the melting temperature,

$$\scriptsize{\begin{eqnarray}&&K_{\text{A}}^{\text{R}}(T)\sim \xi _{+}^{2}(T)\sim {{\text{e}}^{c{{t}^{-\tilde{\mathop{\nu}}\,}}}},\,\,\,t=\frac{T-{{T}_{\text{m}}}}{{{T}_{\text{m}}}},c=\text{nonuniversal}\,\,\,\text{constant},\text{and}\,\,\,\text{jumps}\,\,\,\text{discontinuously}\,\,\,\text{to}\,\,\,\text{zero}\,\,\,\text{at}\,\,\,T_{I}^{+},\\ &&K_{\text{A}}^{\text{R}}(T)=\frac{72}{\pi}\left(1+b\sqrt{|t|}\right),t=\frac{T-{{T}_{I}}}{{{T}_{I}}}\leqslant 0,\\ &&K_{\text{A}}^{\text{R}}(T)=0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,t>0.\end{eqnarray}} \tag{ 158 }$$

When T  >  T_I, the orientational order is short ranged,

$$\begin{eqnarray}&&\langle {{\psi}_{6}}\left(\mathbf{r}\right)\psi _{6}^{\ast}\left(\mathbf{0}\right)\rangle \sim {{\text{e}}^{-r/{{\xi}_{6}}(T)}},\end{eqnarray} \tag{ 159 }$$

where the orientational correlation length ${{\xi}_{6}}(T)$ diverges exponentially as $T\to T_{I}^{+}$,

$$\begin{eqnarray}&&{{\xi}_{6}}(T)=\text{exp}\left(\frac{{{b}_{+}}}{{{\left(T/{{T}_{I}}-1\right)}^{1/2}}}\right).\end{eqnarray} \tag{ 160 }$$

Assuming that the initial value of the bare stiffness ${{K}_{0A}}\left({{l}^{\ast}}\right)$ is larger than the critical value $K_{0A}^{c}$, the renormalized value of the orientational stiffness will flow to a value $K_{\text{A}}^{\text{R}}(T)\geqslant 72/\pi$. For T slightly above the melting temperature ${{T}_{\text{m}}}$, the theory predicts that $K_{\text{A}}^{\text{R}}(T)\sim {{\text{e}}^{2{{l}^{\ast}}}}=\xi _{+}^{2}(T)$ and, at $T_{m}^{-}$, $K_{\text{A}}^{\text{R}}\left(T_{I}^{-}\right)=72/\pi$ [129]. Experimental measurements of $K_{\text{A}}^{\text{R}}(T)$ [158] agree well with these theoretical predictions. If, for some reason, the dislocation core energy ${{E}_{\text{c}}}$ of equation (154) is very small, it is possible that $K_{\text{A}}^{0}\left({{l}^{\ast}}\right)<{{K}_{0c}}$ and will flow immediately to zero. A hexatic phase will not exist, but this seems unlikely. It is also possible that the hexatic phase exists only over a very small range of ${{T}_{\text{m}}}<T\leqslant {{T}_{I}}$. Below ${{T}_{\text{m}}}$, the translational correlation length ${{\xi}_{+}}(T)=\infty$, so that $K_{\text{A}}^{\text{R}}(T)=\infty$, and there is true long range orientational order, as was first pointed out by Mermin [24]. The behavior of the renormalised orientational stiffness ${{K}_{\text{A}}}(T)$ is shown in figure 19.

Zoom In Zoom Out Reset image size

The most important, unexpected, and controversial prediction of the defect theory of melting is that, in 2D, the melting process is not a first order transition from a periodic solid to an isotropic fluid but proceeds by two successive continuous transitions with an orientationally ordered hexatic fluid between the low temperature periodic crystal and the high temperature isotropic fluid. In fact, this scenario has been verified quantitatively by careful experiments on a monolayer of paramagnetic colloidal particles trapped at a water/air interface [155–158]. A recent very large scale simulation of hard disks [159–161] verified the two step melting scenario with an intermediate hexatic phase with a continuous melting transition but a first order hexatic to isotropic fluid transition. However, in a very recent study [162] in the presence of pinning disorder, no signs of any first order character at either transition was observed. Since the original theoretical prediction of continuous two stage melting [129], there have been a number of experimental and numerical studies performed to test the heretical prediction of continuous melting by two continuous transitions with a hexatic fluid intervening between the solid and the isotropic fluid. At the time, conventional wisdom was that melting is a first order transition just as in 3D. The theoretical behavior of the orientational stiffness ${{K}_{\text{A}}}(T)$ is sketched in figure 19 which agrees fairly well with experiment.

Most experiments on 2D melting are performed on monolayer of adsorbed molecules on a supporting substrate such as graphite which has relatively large well oriented domains of binding sites arranged in a periodic array. One must therefore investigate the effect of such a periodic substrate on the possible states of the over layer and on the transitions from the most ordered to the most disordered states of the combined system. Halperin and Nelson investigated this in their seminal papers [128, 129], and the results are summarized here. Both the graphite substrate and the adsorbate overlayer have the same six-fold orientational symmetry but a generic substrate has a different periodicity to the natural periodicity of the adsorbed overlayer. This leads to the elastic free energy [128, 129]

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{1}{2{{k}_{\text{B}}}T}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{q}}{{{(2\pi )}^{2}}}{{u}_{i}}\left(\mathbf{q}\right){{D}_{ij}}\left(\mathbf{q}\right){{u}_{j}}\left(-\mathbf{q}\right)\\ &&+\underset{\mathbf{r},\mathbf{K}}{\mathop \sum}\,{{h}_{\mathbf{K}}}{{\text{e}}^{\text{i}\mathbf{K}\centerdot \mathbf{r}}}\left(1+\text{i}\mathbf{K}\centerdot \mathbf{u}\left(\mathbf{r}\right)+\cdots \right),\\ &&{{D}_{ij}}\left(\mathbf{q}\right)={{\mu}_{\text{R}}}{{q}^{2}}{{\delta}_{ij}}+\left({{\mu}_{\text{R}}}+{{\lambda}_{\text{R}}}\right){{q}_{i}}{{q}_{j}}.\end{eqnarray} \tag{ 161 }$$

One has assumed that the substrate potential is weak and can be written as $V\left(\mathbf{r}\right)=\sum_{\mathbf{K}}{{h}_{\mathbf{K}}}{{\text{e}}^{\text{i}\mathbf{K}\centerdot \left(\mathbf{r}+\mathbf{u}\left(\mathbf{r}\right)\right)}}$ where $\left\{\mathbf{K}\right\}$ are the reciprocal lattice vectors of the periodic substrate potential. The $\left\{\mathbf{K}\right\}$ are assumed incommensurate with the reciprocal lattice vectors $\left\{\mathbf{G}\right\}$ of the adsorbate, assumed to be in a floating solid phase. This can be described by a harmonic free energy with renormalized elastic constants ${{\mu}_{\text{R}}}$ and ${{\lambda}_{\text{R}}}$. Redefine the displacement field [129]

$$\begin{eqnarray}{{u}_{i}}\left(\mathbf{q}\right)\to u_{i}^{\prime}\left(\mathbf{q}\right)&=&{{u}_{i}}\left(\mathbf{q}\right)+\text{i}D_{ij}^{-1}\left(\mathbf{q}\right)\underset{\mathbf{K}}{\mathop \sum}\,{{h}_{\mathbf{K}}}{{K}_{j}}{{\Delta }_{\mathbf{K},\mathbf{q}}},\\ {{\Delta }_{\mathbf{K},\mathbf{q}}}&=&\left\{\begin{array}{c} 1 & \,\,\,\text{if}\,\,\,\mathbf{K}=\mathbf{G}+\mathbf{q} \\ 0 & \,\,\,\text{otherwise}, \end{array}\right. \\ \frac{H}{{{k}_{\text{B}}}T}&=&\frac{1}{2{{k}_{\text{B}}}T}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{q}}{{{(2\pi )}^{2}}}{{u}^{\prime}}\left(\mathbf{q}\right){{D}_{ij}}\left(\mathbf{q}\right)u_{j}^{\prime}\left(-\mathbf{q}\right)-\frac{\Omega }{2}\underset{\mathbf{K}}{\mathop \sum}\,h_{\mathbf{K}}^{2}{{K}_{i}}D_{ij}^{-1}\left(\mathbf{K}\right){{K}_{j}},\end{eqnarray} \tag{ 162 }$$

where $\mathbf{G}$ is a reciprocal lattice vector of the adsorbed layer. The second term of equation (162) can be rewritten as

$$\begin{eqnarray}&&C(\theta )=-\frac{\Omega }{2}\underset{\mathbf{K}}{\mathop \sum}\,\underset{s=1}{\overset{2}{\mathop \sum}}\,h_{\mathbf{K}}^{2}{{\left(\frac{\mathbf{K}\centerdot {{\widehat{\boldsymbol{\epsilon }}}_{s}}\left(\mathbf{K}\right)}{{{\omega}_{s}}\left(\mathbf{K}\right)}\right)}^{2}},\end{eqnarray} \tag{ 163 }$$

where ${{\widehat{\boldsymbol{\epsilon }}}_{s}}\left(\mathbf{K}\right)$ is the ${{s}^{\text{th}}}$ polarization vector and ${{\omega}_{s}}\left(\mathbf{K}\right)$ is the eigenfrequency of the matrix ${{D}_{ij}}\left(\mathbf{K}\right)$.

For a hexagonal over layer with lattice vectors $\left\{\mathbf{G}\right\}$ on a hexagonal substrate with incommensurate lattice vectors $\left\{\mathbf{K}\right\}$, perfect alignment of the crystal axes with the substrate corresponds to $\theta =0$, and one can write the Fourier expansion of $C(\theta )$ as

$$\begin{eqnarray}&&C(\theta )=\Omega \underset{k=0}{\overset{\infty}{\mathop \sum}}\,{{c}_{k}}\text{cos}(6k\theta ),\end{eqnarray} \tag{ 164 }$$

and, for a hexagonal adsorbate on a square substrate

$$\begin{eqnarray}&&C(\theta )=\Omega \underset{k=0}{\overset{\infty}{\mathop \sum}}\,{{c}_{k}}\text{cos}(12k\theta ).\end{eqnarray} \tag{ 165 }$$

Assuming that the relative orientations of the substrate and adsorbed layer vary slowly in space, θ is the bond angle $\theta \left(\mathbf{r}\right)=\frac{1}{2}\left(\partial u_{x}^{\prime}/\partial y-\partial u_{y}^{\prime}/\partial x\right)$ so that, for small deviations θ from perfect alignment, one obtains the effective elastic Hamiltonian [129, 130]

$$\begin{eqnarray}&&\frac{{{H}_{\text{E}}}}{{{k}_{\text{B}}}T}=\frac{1}{2}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{{{a}_{0}}}\left[2\tilde{\mathop{\mu}}\,u_{ij}^{\prime 2}+\tilde{\mathop{\lambda}}\,u_{kk}^{\prime 2}+\tilde{\mathop{\gamma}}\,{{\left({{\partial}_{y}}u_{x}^{\prime}-{{\partial}_{x}}u_{y}^{\prime}\right)}^{2}}\right], \\ &&\,\,\,\,\,\,\,\,\,\,\,\,=\frac{1}{2}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{a_{0}^{2}}\left[2\tilde{\mathop{\mu}}\,\phi _{ij}^{2}+\tilde{\mathop{\lambda}}\,\phi _{kk}^{2}+\tilde{\mathop{\gamma}}\,{{\left({{\partial}_{y}}{{\phi}_{x}}-{{\partial}_{x}}{{\phi}_{y}}\right)}^{2}}\right]+\frac{{{H}_{\text{D}}}}{{{k}_{\text{B}}}T},\\ &&\frac{{{H}_{\text{D}}}}{{{k}_{\text{B}}}T}=-\frac{1}{8\pi}\underset{\mathbf{r}\ne {{\mathbf{r}}^{\prime}}}{\mathop \sum}\,{{b}_{i}}\left(\mathbf{r}\right){{b}_{j}}\left({{\mathbf{r}}^{\prime}}\right)\left({{K}_{1}}\text{ln}\frac{|\mathbf{r}-{{\mathbf{r}}^{\prime}}|}{a}{{\delta}_{ij}}\\ &&-{{K}_{2}}\frac{{{\left(\mathbf{r}-{{\mathbf{r}}^{\prime}}\right)}_{i}}{{\left(\mathbf{r}-{{\mathbf{r}}^{\prime}}\right)}_{j}}}{|\mathbf{r}-{{\mathbf{r}}^{\prime}}{{|}^{2}}}\right)+\frac{{{E}_{\text{c}}}}{{{k}_{\text{B}}}T}\underset{\mathbf{r}}{\mathop \sum}\,{{\mathbf{b}}^{2}}\left(\mathbf{r}\right)\\ &&=\frac{1}{2}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{q}}{{{(2\pi )}^{2}}}\left[\frac{{{K}_{1}}+{{K}_{2}}}{2}\frac{1}{{{q}^{2}}}\left({{\delta}_{ij}}-\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{2}}}\right)\\ &&+\frac{{{K}_{1}}-{{K}_{2}}}{2}\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{4}}}+\frac{2{{E}_{\text{c}}}a_{0}^{2}}{{{k}_{\text{B}}}T}{{\delta}_{ij}}\right]{{b}_{i}}\left(\mathbf{q}\right){{b}_{j}}\left(-\mathbf{q}\right),\\ &&{{K}_{1}}=\frac{4\tilde{\mathop{\mu}}\,\left(\tilde{\mathop{\mu}}\,+\tilde{\mathop{\lambda}}\,\right)}{2\tilde{\mathop{\mu}}\,+\tilde{\mathop{\lambda}}\,}+\frac{4\tilde{\mathop{\mu}}\,\tilde{\mathop{\gamma}}\,}{\tilde{\mathop{\mu}}\,+\tilde{\mathop{\gamma}}\,},\\ &&{{K}_{2}}=\frac{4\tilde{\mathop{\mu}}\,\left(\tilde{\mathop{\mu}}\,+\tilde{\mathop{\lambda}}\,\right)}{2\tilde{\mathop{\mu}}\,+\tilde{\mathop{\lambda}}\,}-\frac{4\tilde{\mathop{\mu}}\,\tilde{\mathop{\gamma}}\,}{\tilde{\mathop{\mu}}\,+\tilde{\mathop{\gamma}}\,}.\end{eqnarray} \tag{ 166 }$$

To translate the coupling constants K₁ and K₂ into the language of [130] where the recursion relations were first worked out for this general dislocation Hamiltonian, one has

$$\begin{eqnarray}&&{{K}_{1}}=8{{\pi}^{2}}{{K}^{r}},\\ &&{{K}_{2}}=8{{\pi}^{2}}{{K}^{\theta}}.\end{eqnarray} \tag{ 167 }$$

From the recursion relations of [130], it follows that

$$\begin{eqnarray}&&\frac{\text{d}{{K}_{1}}(l)}{\text{d}l}=-\frac{3}{4}\pi {{y}^{2}}{{\text{e}}^{{{K}_{2}}/8\pi}}\left[\left(K_{1}^{2}+K_{2}^{2}\right){{I}_{0}}\left(\frac{{{K}_{2}}}{8\pi}\right)-{{K}_{1}}{{K}_{2}}{{I}_{1}}\left(\frac{{{K}_{2}}}{8\pi}\right)\right]+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{\text{d}{{K}_{2}}(l)}{\text{d}l}=-\frac{3}{4}\pi {{y}^{2}}\left[2{{K}_{1}}{{K}_{2}}{{I}_{0}}\left(\frac{{{K}_{2}}}{8\pi}\right)-\frac{1}{2}\left(K_{1}^{2}+K_{2}^{2}\right){{I}_{1}}\left(\frac{{{K}_{2}}}{8\pi}\right)\right]+\mathcal{O}\left(\,{{y}^{3}}\right),\\ &&\frac{\text{d}y(l)}{\text{d}l}=\left(2-\frac{{{K}_{1}}}{8\pi}\right)y+2\pi {{y}^{2}}{{\text{e}}^{{{K}_{2}}/16\pi}}{{I}_{0}}\left(\frac{{{K}_{2}}}{8\pi}\right)+\mathcal{O}\left(\,{{y}^{3}}\right).\end{eqnarray} \tag{ 168 }$$

From this, it follows that dislocation unbinding melting takes place when the fugacity ${{y}_{0}}={{\text{e}}^{-{{E}_{\text{c}}}/{{k}_{\text{B}}}T}}$ becomes relevant, which is controlled by K₁. The renormalized stiffness coefficient at $T_{\text{m}}^{-}$ is $K_{1}^{\text{R}}\left(T_{\text{m}}^{-}\right)=16\pi$ and the value of the exponent $\tilde{\nu}$ depends on $K_{2}^{\text{R}}\left(T_{M}^{-}\right)/K_{1}^{\text{R}}\left(T_{M}^{-}\right)$. From (168) it is clear that K₂(l) flows to some non-universal finite value at $l=\infty$ and, when K₂  =  0, $\tilde{\nu}=2/5$ as found by Nelson [127] and melting on an incommensurate substrate is very similar to a smooth substrate. As pointed out by Nelson [77], the situation above ${{T}_{\text{m}}}$ is rather different. The orientational bias can be modeled by adding a term ${{h}_{p}}\text{cos}p\theta$ to the hexatic free energy

$$\begin{eqnarray}&&\frac{{{F}_{\text{A}}}}{{{k}_{\text{B}}}T}=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\left(\frac{{{K}_{\text{A}}}(T)}{2}{{(\nabla \theta )}^{2}}-{{h}_{p}}\text{cos}\,p\theta \right).\end{eqnarray} \tag{ 169 }$$

Clearly, for a hexagonal overlayer with orientational order parameter ${{\psi}_{6}}={{\text{e}}^{6\text{i}\theta}}$ on a hexagonal substrate with p  =  6, h₆ acts like a uniform magnetic field inducing long range orientational order for all values of $T>{{T}_{\text{m}}}$, wiping out the hexatic to isotropic fluid transition on a smooth substrate at T_I.

When the substrate and adsorbed layer reciprocal lattice vectors $\left\{\mathbf{G}\right\}$ and $\left\{\mathbf{K}\right\}$ have a set $\left\{\mathbf{M}\right\}$ in common, the over layer may become locked to the substrate as a commensurate phase. This situation can be studied by considering the smallest common reciprocal lattice vector of length ${{M}_{0}}=|{{\mathbf{M}}_{0}}|$ and writing the free energy as [129]

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\frac{{{H}_{\text{E}}}}{{{k}_{\text{B}}}T}+\underset{|\mathbf{M}|={{M}_{0}}}{\mathop \sum}\,{{h}_{\mathbf{M}}}\underset{\mathbf{r}}{\mathop \sum}\,{{\text{e}}^{\text{i}\mathbf{M}\centerdot \mathbf{u}\left(\mathbf{r}\right)}},\end{eqnarray} \tag{ 170 }$$

where it is assumed that the effects of all incommensurate lattice vectors have been incorporated into the elastic part ${{H}_{\text{E}}}$. This can be done because these and terms in the substrate potential

$$\begin{eqnarray}&&\frac{{{H}_{S}}}{{{k}_{\text{B}}}T}=\underset{\mathbf{K}}{\mathop \sum}\,{{h}_{\mathbf{K}}}\underset{\mathbf{r}}{\mathop \sum}\,{{\text{e}}^{\text{i}\mathbf{K}\centerdot \mathbf{u}\left(\mathbf{r}\right)}}\end{eqnarray} \tag{ 171 }$$

with $|\mathbf{K}|>{{M}_{0}}$ are less relevant in the renormalization group sense than those with $|\mathbf{K}|={{M}_{0}}$.

The relevance of the ${{h}_{\mathbf{M}}}$ is obtained from the correlation function computed from the renormalized elastic free energy ${{H}_{\text{E}}}/\left({{k}_{\text{B}}}T\right)$ with ${{h}_{\mathbf{M}}}=0$,

$$\begin{eqnarray}&&{{\langle {{\text{e}}^{\text{i}\mathbf{M}\centerdot \left[\mathbf{u}\left(\mathbf{r}\right)-\mathbf{u}(0)\right]}}\rangle}_{{{H}_{\text{E}}}}}\sim {{r}^{-{{\eta}_{\mathbf{M}}}(T)}},\\ &&{{\eta}_{\mathbf{M}}}(T)=\frac{|\mathbf{M}{{|}^{2}}{{k}_{\text{B}}}T}{4\pi}\frac{3{{\mu}_{\text{R}}}(T)+{{\lambda}_{\text{R}}}(T)+{{\gamma}_{\text{R}}}(T)}{\left({{\mu}_{\text{R}}}(T)+{{\gamma}_{\text{R}}}(T)\right)\left(2{{\mu}_{\text{R}}}(T)+{{\lambda}_{\text{R}}}(T)\right)}.\end{eqnarray} \tag{ 172 }$$

From this, the renormalization group eigenvalue of ${{h}_{{{M}_{0}}}}(l)={{h}_{{{M}_{0}}}}(0){{\text{e}}^{{{\lambda}_{{{M}_{0}}}}l}}$ is

$$\begin{eqnarray}&&{{\lambda}_{{{M}_{0}}}}(T)=2-\frac{1}{2}{{\eta}_{{{M}_{0}}}}(T)\\ &&=2-\frac{M_{0}^{2}{{k}_{\text{B}}}T}{8\pi}\frac{3{{\mu}_{\text{R}}}(T)+{{\lambda}_{\text{R}}}(T)+{{\gamma}_{\text{R}}}(T)}{\left({{\mu}_{\text{R}}}(T)+{{\gamma}_{\text{R}}}(T)\right)\left(2{{\mu}_{\text{R}}}(T)+{{\lambda}_{\text{R}}}(T)\right)}.\end{eqnarray} \tag{ 173 }$$

From (173), at sufficiently low T, all the ${{h}_{\mathbf{M}}}$ are relevant so that the over layer is locked to the substrate, implying that a lattice gas description is more appropriate [163]. If the substrate is sufficiently coarse, ${{\lambda}_{{{M}_{0}}}}\left({{T}_{\text{m}}}\right)>0$, so that the adsorbate is locked to the substrate up to the melting temperature ${{T}_{\text{m}}}$, for ${{h}_{{{\mathbf{M}}_{0}}}}\to 0$ so that dislocation melting is inappropriate at any T. On the other hand, for a fine substrate mesh, M₀ is large and ${{\lambda}_{{{M}_{0}}}}\left({{T}_{\text{m}}}\right)<0$, so that the overlayer is unlocked from the substrate in a temperature range ${{T}_{l}}<T<{{T}_{\text{m}}}$, leading to a floating incommensurate solid which may be accidentally commensurate with the substrate. Halperin and Nelson [129] have derived a criterion for the unlocking temperature ${{T}_{\text{L}}}<{{T}_{\text{m}}}$. The consequence of this is that ideal 2D melting is still possible on a periodic crystal substrate provided the adsorbate is a floating incommensurate solid.

The most dramatic effect of the substrate potential is on the hexatic fluid phase $T>{{T}_{\text{m}}}$. One performs an analysis for the stiffness constant ${{K}_{\text{A}}}(T)$, defined in equation (151) in exactly the same way, using a Debye–Hückel [78] approximation and the dislocation Hamiltonian of equation (166) to obtain

$$\begin{eqnarray}K_{\text{A}}^{-1}(T)&=&\underset{\mathbf{q}\to 0}{\mathop{\lim}}\,\frac{a_{0}^{2}}{\Omega }\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{2}}}\langle {{b}_{i}}\left(\mathbf{q}\right){{b}_{j}}\left(-\mathbf{q}\right)\rangle \\ \langle {{b}_{i}}\left(\mathbf{q}\right){{b}_{j}}\left(-\mathbf{q}\right)\rangle &=&\frac{2{{q}^{2}}}{{{K}_{1}}+{{K}_{2}}+\frac{4{{E}_{\text{c}}}{{q}^{2}}a_{0}^{2}}{{{k}_{\text{B}}}T}}\left({{\delta}_{ij}}-\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{2}}}\right)+\frac{2{{q}^{2}}}{{{K}_{1}}-{{K}_{2}}+\frac{4{{E}_{\text{c}}}{{q}^{2}}a_{0}^{2}}{{{k}_{\text{B}}}T}}\frac{{{q}_{i}}{{q}_{j}}}{{{q}^{2}}},\\ K_{\text{A}}^{-1}(T)&=&\underset{\mathbf{q}\to 0}{\mathop{\lim}}\,\frac{2{{q}^{2}}}{{{K}_{1}}-{{K}_{2}}+\frac{4{{E}_{\text{c}}}{{q}^{2}}a_{0}^{2}}{{{k}_{\text{B}}}T}}=\left\{\begin{array}{c} 0, & \,\,\,\left({{K}_{1}}\ne {{K}_{2}}\right) \\ \frac{{{k}_{\text{B}}}T}{2{{E}_{\text{c}}}a_{0}^{2}}, & \,\,\,\left({{K}_{1}}={{K}_{2}}\right) \end{array}\right. \end{eqnarray} \tag{ 174 }$$

This means that the orientational stiffness constant ${{K}_{\text{A}}}(T)=\infty$ for an adsorbate on an incommensurate substrate, which, in turn, means that there is true long range orientational order for all $T>{{T}_{\text{m}}}$ and that there is no high temperature transition to the expected isotropic fluid. Thus, the experimental investigation seems to be a very difficult proposition and there is an excellent discussion in the review article by Strandburg [142] with extensive references to many aspects of the melting problem.

This difficulty was finally overcome in 1999 [157] by trapping a monolayer of colloidal particles at an air/water interface, doped with paramagnetic ions in a magnetic field B normal to the interface. The interaction between the particles is characterized by $\Gamma=\left({{\mu}_{0}}/4\pi \right)\left(\chi {{B}^{2}}\right){{(\pi n)}^{3/2}}/\left({{k}_{\text{B}}}T\right)$ where n is the areal number density and χ the magnetic susceptibility. This system is as close to the conditions of the theory as it is possible to get as the substrate is completely isotropic and free of any destructive perturbations. The particle positions are recorded by video camera so that all correlation functions can be readily calculated and compared with theory. The outcome is that, in this system, measurements agree quantitatively with the KTHNY theory in all respects with no adjustable parameters [155, 156].

One of the major experimental difficulties is equilibrating the system, despite it being only about 10⁵ particles. In the high field crystalline phase, the equilibration time is several days [157] even with the system being jiggled by applying a small ac ${{B}_{||}}$ in the plane of the system to anneal out the large out of equilibrium concentration of defects such as grain boundaries, interstitials and vacancies. Eventually, the system consists of a single domain with less than 0.1% of dislocations. To melt this 2D crystal, the temperature ${{\Gamma}^{-1}}$ is increased by decreasing B in small steps equilibrating for about one hour after each step [157]. This procedure results in ${{C}_{\text{G}}}\left(\mathbf{r}\right)$ decaying as a power law for $T\leqslant {{T}_{\text{m}}}$ and exponentially for $T>{{T}_{\text{m}}}$, and similarly for ${{C}_{6}}\left(\mathbf{r}\right)$, as shown in figure 20

$$\begin{eqnarray}{{C}_{\text{G}}}\left(\mathbf{r}\right)\sim \left\{\begin{array}{c} {{r}^{-{{\eta}_{\text{G}}}(T)}}, & \,\,\,T\leqslant {{T}_{\text{m}}}, \\ {{\text{e}}^{-r/{{\xi}_{+}}(T)}}, & \,\,T>{{T}_{\text{m}}} \end{array}\right. \\ {{C}_{6}}\left(\mathbf{r}\right)\sim \left\{\begin{array}{c} {{r}^{-{{\eta}_{6}}(T)}}, & \,\,\,{{T}_{\text{m}}}<T\leqslant {{T}_{I}}, \\ {{\text{e}}^{-r/{{\xi}_{6}}(T)}}, & \,\,T>{{T}_{I}} \end{array},\right. \end{eqnarray} \tag{ 175 }$$

where ${{\eta}_{\mathbf{G}}}(T)$ is given in (143), ${{\xi}_{+}}(T)$ in (146), ${{\eta}_{6}}(T)$ in (157) and ${{\xi}_{6}}(T)$ in (160). The experimental results for ${{C}_{\mathbf{G}}}\left(\mathbf{r}\right)$ and ${{C}_{6}}\left(\mathbf{r}\right)$ are shown in figure 20 and are consistent with the theoretical predictions. Note that ${{C}_{\mathbf{G}}}\left(\mathbf{r}\right)$ in the solid phase decays as a power law ${{r}^{-{{\eta}_{\mathbf{G}}}(T)}}$ with ${{\eta}_{\mathbf{G}}}\left({{T}_{\text{m}}}\right)=1/3$. This agrees with theory because the measured power law is actually ${{r}^{-{{\eta}_{{{\text{G}}_{0}}}}}}$ where G₀ is the magnitude of the smallest reciprocal lattice vector for a hexagonal lattice ${{G}_{0}}=4\pi /\left({{a}_{0}}\sqrt{3}\right)$.

Zoom In Zoom Out Reset image size

Measured values of the Young's modulus ${{\tilde{K}}_{\text{R}}}(T)$ compared with theory are shown in the upper panel of figure 21 and for the elastic constants ${{\tilde{\mu}}_{\text{R}}}(T)$ and ${{\tilde{\lambda}}_{\text{R}}}(T)+2{{\tilde{\mu}}_{\text{R}}}(T)$ in the lower panel which also agree well with theory. The experimental systems are electrons on ⁴He [150] and paramagnetic colloidal particles at an air/water interface [144] as these are the closest to the conditions of theory and are free from unwanted substrate effects. As shown in (145), for an incompressible lattice, ${{\eta}_{{{\text{G}}_{0}}}}=1/3$ which is close to the theoretical value for particles interacting by a 1/r³ repulsive potential. As shown in [144], the exponent

$$\begin{eqnarray}&&{{\eta}_{{{\text{G}}_{0}}}}(T)=\frac{{{{\tilde{\lambda}}}_{\text{R}}}(T)+3{{{\tilde{\mu}}}_{\text{R}}}(T)}{\pi \Gamma{{{\tilde{\mu}}}_{\text{R}}}(T)\left({{{\tilde{\lambda}}}_{\text{R}}}(T)+2{{{\tilde{\mu}}}_{\text{R}}}(T)\right)}{{\to}_{T=T_{m}^{-}}}\frac{0.985}{4{{\pi}^{2}}},\end{eqnarray} \tag{ 176 }$$

so that, at $T_{\text{m}}^{-}$, ${{\Gamma}_{m}}{{\eta}_{{{\text{G}}_{0}}}}\left(T_{\text{m}}^{-}\right)=1.496$, as shown in figure 22.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the main interests in melting in 2D is the order of the transition. It is generally agreed that the solid/liquid transition in 3D is first order characterized by a latent heat. In fact, the melting transition in 3D is considered as the paradigm of a discontinuous first order transition. The KTHNY theory of melting in 2D predicting successive continuous transitions from the low T solid to the high T isotropic liquid came as a heretical claim which needed to be demonstrated by experiment to be false. However, it turned out to be very difficult to either verify or falsify the KTHNY theory because there were no experimental systems available which obeyed the very stringent conditions of the theory which requires the adsorbate layer to be constrained by a featureless substrate potential of the form $V(x,y,z)=-{{V}_{0}}\delta \left(z-{{z}_{0}}\right)$ with ${{V}_{0}}\gg {{k}_{\text{B}}}T$ which restricts the adsorbate to a flat layer at z  =  z₀ where the substrate surface is at z  =  0. The theory requires that the potential V(x, y, z) is independent of the substrate coordinates (x, y) and this can be realized only in two situations: (i) electrons trapped in a layer above the surface of ⁴He [150, 153, 154] and (ii) paramagnetic colloidal particles trapped at a water/air interface [157].

Experiments on orientational order in the hexatic phase have also been performed on a layer of paramagnetic colloidal particles interacting by a 1/r³ repulsive energy [155, 157, 158] with very good agreement with the theoretical predictions [129]. Excellent reviews summarizing the present status of the transitions between the three phases of particles in two dimensions with 1/r and 1/r³ repulsive interactions are in [144, 156]. Unfortunately, orientational order and the hexatic phase cannot be studied in the electrons on helium system because measurements can be made only on the frequency spectrum of resonances of the coupled electron—ripplon modes [95, 150, 151]. The comparison between theory and experiment for the orientational stiffness ${{\tilde{K}}_{\text{A}}}(T)$ is shown in figure 24. The theoretical renormalization group predictions of the KTHNY theory are the solid line and the experimental measurements are the points. To the author's eye, the agreement is very good, considering that there are no adjustable parameters in the fits.

The orientational correlation function ${{G}_{6}}\left(\mathbf{r}\right)$ is investigated and the results are summarized in figure 23 where it is seen that ${{G}_{6}}(r)\sim {{r}^{-{{\eta}_{6}}(T)}}$ for $r/a\leqslant 30$, which is limited by the sample size. ${{G}_{6}}\left(\mathbf{r}\right)$ is consistent with the theoretical predictions with a power law decay with ${{\eta}_{6}}(T)\leqslant 1/4$ for ${{T}_{\text{m}}}<T\leqslant {{T}_{I}}$ and exponential decay for T  >  T_I. The correlation length ${{\xi}_{6}}(T)$ is shown in figure 24(a) and is consistent with the theoretical prediction ${{\xi}_{6}}(T)\sim {{\text{e}}^{b{{t}^{-1/2}}}}$ although the experimental points are scattered about the theoretical solid line. The orientational stiffness constant, ${{\tilde{K}}_{\text{A}}}(T)$, is shown in figure 24(b) where agreement with theory, $\tilde{K}_{\text{A}}^{\text{R}}\left(T_{I}^{+}\right)=72/\pi$ and $\tilde{K}_{\text{A}}^{\text{R}}\left(T_{m}^{+}\right)\sim \text{exp}\left(b{{t}^{-1/2}}\right)$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On balance, the agreement with theory verifies the KTHNY theory of melting at least for this particular system as there are no contradictions. However, it must be mentioned that the actual critical region in which the exponential divergence of the correlation length ${{\xi}_{+}}(t)\sim {{\text{e}}^{b{{t}^{-\tilde{\nu}}}}}$ sets in is estimated to be $\mathcal{O}\left({{10}^{8}}\right)$ lattice spacings [164]. This is much larger than accessible experimental sizes which casts grave doubt about the significance of the fitting of data to the theoretical form of ${{\xi}_{+}}$! In particular, because of the limited sample sizes, these experiments are unable to exclude weak first order transitions seen in large scale simulations of hard disks [159–161].

The transition the 2D planar rotor and melting models are argued to be driven by vortex and dislocation unbinding [19] which lead to continuous transitions. However, symmetry arguments do not exclude the possibility that, in some systems, these defect driven transitions are first order. The earliest study of melting in 2D which gave convincing evidence of a transition was carried out by Alder and Wainwright [31] and there have been many others since then.

However, a recent very large scale simulation of a hard disk system with N  =  1024² disks [160, 161] obtains results which agree with the theoretical predictions for the solid—hexatic transition but show that the hexatic—isotropic fluid transition is weakly first order with a discontinuity in the density of about 0.015. Although this scenario sounds unlikely, it is entirely possible. It is known that the KT transition can be made first order by reducing the vortex core energy ${{E}_{\text{c}}}$ to below some critical value [165–167]. In [165], a modified planar rotor Hamiltonian on a square lattice was proposed

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}={{K}_{0}}(T)\underset{<ij>}{\mathop \sum}\,\left(1-{{\left[\text{co}{{\text{s}}^{2}}\left(\left({{\theta}_{i}}-{{\theta}_{j}}\right)/2\right)\right]}^{{{p}^{2}}}}\right).\end{eqnarray} \tag{ 177 }$$

This is precisely the planar rotor model when p  =  1 which has a continuous transition driven by vortex unbinding while, for large p, $V(\theta )\approx 2J$ for $\theta \geqslant \pi /p$ and $V(\theta )\approx \frac{1}{2}J{{p}^{2}}{{\theta}^{2}}$ for $\theta \ll \pi /p$ [166]. The physical argument for a first order transition for $p\gg 1$ is that the Hamiltonian of (177) resembles that of a Potts model [105] with a large number of states which is known to undergo a first order transition [107]. Vortex unbinding is at ${{k}_{\text{B}}}{{T}_{\text{KT}}}\approx J{{p}^{2}}$ and vacancy condensation at ${{k}_{\text{B}}}{{T}_{\text{D}}}\approx V(\pi )=2J$. For p sufficiently large, one expects that ${{T}_{\text{D}}}<{{T}_{\text{KT}}}$ so that the continuous KT transition is preempted by a first order vacancy condensation [166]. These physical arguments are verified by simulations [165–167] although the system sizes are somewhat limited and proved rigorously in [168]. Also, the vortex core energy ${{E}_{\text{c}}}$ is reduced as p increases which indicates that a defect description with an appropriate defect fugacity y₀ is sufficient to account for both types of transition. In particular, this can explain the first order hexatic-isotropic fluid transition observed in simulations of the hard disk system [159–161]. It is worth noting that this simulation of up to N  =  2²⁰ hard discs is, in fact, larger than any 2D experimental system. Of course, it differs somewhat from a real particle system since these have interactions beyond the hard core repulsion.

In some of the earliest Monte Carlo (MC) simulations of hard spheres and discs interacting by a Lennard-Jones potential $V(r)=4\epsilon \left[{{(\sigma /r)}^{12}}-{{(\sigma /r)}^{6}}\right]$ [26], it was found that a system of N  =  56 particles in 2D and up to 256 in 3D with periodic boundary conditions showed indications of a gas/liquid transition in 3D but not in 2D. Studies of the 3D hard sphere system by MC [27] and by molecular dynamics (MD) [28–30] agreed on the equation of state and on a first order transition in 3D with small finite size effects. Later, a MD simulation of N  =  870 hard disks in 2D was done [31] where a first order melting/freezing transition was observed. Since these early studies, there has been over fifty years of numerical study on the apparently simple hard disk system culminating in the massive simulations by Krauth and coworkers [159–161]. It was discovered that 2D melting is neither first order nor KTHNY but has the low T crystal melts to a hexatic fluid by a continuous KTHNY transition followed by a very weak first order transition between the hexatic and isotropic fluids.

Since Alder and Wainwright's simulation of the hard disk system [31], many simulations of 2D hard disk and Lennard-Jones systems have been performed, and excellent reviews of the situation before 1985 are by [142, 169]. The great majority of simulations concluded that melting in 2D is by a single first order transition, with few exceptions e.g. [170] who pointed out that his simulations were unable to distinguish between a single weak first order and the KTHNY scenario. It is of interest to note that the existence of the hexatic fluid phase has been confirmed by simulation only very recently [159, 161]. Even relatively large simulations have not been able to yield an unambiguous confirmation of the existence of a hexatic phase [171–178] and conclude that 2D melting is a single first order transition with the exception of [176] where the conclusion was that 2D melting takes place by a single continuous transition. This is mainly due to size limitations, the largest system simulated was N  =  2¹⁴ disks, but also the underlying expectation that 2D melting should be similar to 3D first order melting. Thus, usually only orientational order was followed and translational order assumed to be slaved to this but, even when the two orders were considered, the hexatic phase was not seen unambiguously. In a simulation of 2²⁰ hard disks, Jaster [179] obtained evidence of the hexatic phase by finite size scaling but the nature of the hexatic—isotropic fluid transition could not be determined [180]. This was finally settled very recently by Krauth and coworkers [159, 161], who found a very narrow hexatic phase reached from the solid by KTHNY melting followed by a very weak first order transition to the isotropic fluid. The first order nature is verified by the finite size scaling of the interface free energy $\Delta f(N)\sim {{N}^{-1/2}}$ for particle number N  =  2ⁿ, n  =  14, 16, 18, 20, provided the correlation length $\xi \ll \sqrt{N}$ [181]. The hard disk simulations [160, 161] obey this criterion thus establishing the first order nature of the hexatic-isotropic fluid transition. The crystal-hexatic transition was also checked to be KTHNY by the behavior of the density–density correlation function, $\langle {{\rho}_{{{\text{G}}_{0}}}}\mathbf{r}{{\rho}_{{{\text{G}}_{0}}}}(0)\rangle \sim {{r}^{-\eta}}$ in the solid phase and $\sim {{\text{e}}^{-r/\xi}}$ in the hexatic phase with $\eta =1/3$ at the hexatic-crystal transition [160, 161].

An elastic solid can be described by interacting dislocations and disclinations, provided the concentration of point defects such as vacancies and interstitials is not too large [19, 128, 129]. These point defects are irrelevant in the renormalization group sense and can be absorbed into the elastic constants $\mu,\lambda$. Thus, a natural way to study melting numerically is to simulate the dislocation free energy of (121) with parameters K₀ and dislocation core energy ${{E}_{\text{c}}}$. When ${{E}_{\text{c}}}\to 0$, dislocations are expected to form grain boundaries between crystallites of different orientation and Chui [182, 183] constructed a theory of grain boundary melting which gives a first order transition where both translational and orientational order vanish simultaneously. Saito [184, 185] simulated a system of dislocations of (121) for different values of the core energy ${{E}_{\text{c}}}$. It is found that, for large ${{E}_{\text{c}}}>E_{\text{c}}^{\ast}\approx 1.22{{K}_{0}}$ a continuous KYHNY melting transition occurs at ${{K}_{\text{R}}}=16\pi$ [129] but when ${{E}_{\text{c}}}<E_{\text{c}}^{\ast}$, there is first order melting directly to an isotropic fluid [184, 185]. When ${{E}_{\text{c}}}<E_{\text{c}}^{\ast}$, the dislocations organize themselves into grain boundary loops and first order melting results [182, 183]. Of course, these point defects are important for the dynamics of the collective defects like dislocations. Dislocation glide conserves particle number but climb or motion parallel to the Burgers vector requires diffusion of vacancies and interstitials and is therefore very slow [186]. Our purely elastic description does not take such effects into consideration and to do this would need more of a microscopic description.

In the simpler system of superfluid ⁴He with some non-superfluid defects, ³He atoms, there is a continuous transition for low impurity concentration which becomes first order for larger concentrations [86]. In this work, the impurities are incorporated into the planar rotor model on a lattice by introducing a t_i  =  0, 1 to represent a ³He atom by t_i  =  0 and a ⁴He atom by t_i  =  1. The resulting Hamiltonian is [86]

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}={{K}_{1}}\underset{<ij>}{\mathop \sum}\,{{t}_{i}}{{t}_{j}}\left[\left(1-\text{cos}\left({{\theta}_{i}}-{{\theta}_{j}}\right)\right)\right]-{{K}_{2}}\underset{<ij>}{\mathop \sum}\,{{t}_{i}}{{t}_{j}}+\Delta \underset{i}{\mathop \sum}\,{{t}_{i}}.\end{eqnarray}$$

This was studied on a triangular lattice by an approximate Migdal [187] type bond moving renormalization group procedure due to Kadanoff [188]. As expected, both continuous and discontinuous transitions appeared. Although never carried out with impurities, a calculation along these lines has been performed by Nelson [127]. It seems that this would be a relatively simple way to study the effects of point defects (vacancies and interstitials) on the continuous dislocation and disclination unbinding transitions and to explore possible melting scenarios. This goes beyond the purely elastic theories of 2D melting which require the introduction of disclinations to induce a first order transition. One outstanding problem with this approach is that the first order transition induced by disclinations [189, 190] requires going beyond linear elasticity and causes the simultaneous destruction of both translational and orientational order. This seems rather restrictive and certainly does not always happen in real systems. A typical phase diagram from purely elastic theory is sketched in figure 25.

Zoom In Zoom Out Reset image size

In a series of papers [189–196], Kleinert and coworkers developed and simulated a model for 2D melting including an important additional parameter which allows for first order melting within a theory based on elasticity. By allowing for higher order elastic effects [195] one can incorporate local rotations

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\mathop{\int}^{}{{\text{d}}^{2}}r\left[2\tilde{\mathop{\mu}}\,u_{ij}^{2}+\tilde{\mathop{\lambda}}\,u_{kk}^{2}+\left(2\tilde{\mathop{\mu}}\,+\tilde{\mathop{\lambda}}\,\right){{l}^{\prime 2}}{{\left({{\partial}_{i}}{{u}_{kk}}\right)}^{2}}+4\tilde{\mathop{\mu}}\,{{l}^{2}}{{\left({{\partial}_{i}}\omega \right)}^{2}}+\cdots \right],\end{eqnarray} \tag{ 178 }$$

where u_ij is the strain tensor and $\omega =\frac{1}{2}\left({{\partial}_{1}}{{u}_{2}}-{{\partial}_{2}}{{u}_{1}}\right)$ is the local rotation field. The parameters l² and ${{l}^{\prime 2}}$ appear in $\omega \left(\mathbf{q}\right)$ relations as [195]

$$\begin{eqnarray}&&\omega _{\text{T}}^{2}=\mu \left({{q}^{2}}+{{l}^{2}}{{q}^{4}}+\cdots \right)\\ &&\omega _{\text{L}}^{2}=(2\mu +\lambda )\left({{q}^{2}}+{{l}^{\prime 2}}{{q}^{4}}+\cdots \right)\end{eqnarray} \tag{ 179 }$$

However, Kleinert [195, 196] has shown that the parameter ${{l}^{\prime 2}}$ is irrelevant for melting and considers only l², which is a measure of rotational stiffness. A model combining both mechanisms leading to first order melting [192] and [86] might be more realistic than existing models.

Most simulations of melting of an elastic solid have been performed in the dual roughening representation, which is an exact transformation of the interacting dislocations and disclinations with long range interactions into a roughening model on the dual lattice with short range interactions [142, 193, 197, 198]. These models have a phase diagram of figure 25 with both continuous melting or a first order transition directly to the isotropic fluid and have a definite computational advantage of nearest neighbor interactions but incorporating local defects is more awkward.

The KTHNY melting theory inspired a large experimental effort to verify or falsify the two continuous transition scenario, summarized by Strandburg in a comprehensive review [142]. There have been several studies of liquid crystal films [96–103, 121, 199–201]. In some of these, melting seems to proceed according to the KTHNY scenario [121] but these freely suspended systems are much richer than theory suggests with many different possible phases depending on temperature and film thickness [201]. One interesting observation is that, as the film reaches its minimum thickness of two layers, a hexatic phase is observed which seems to behave as predicted by theory [123, 202], although the evidence for a separate hexatic phase is not conclusive.

Adsorbed gases on substrates provide a huge variety of 2D systems but the interplay of the periodic elastic adsorbed layer with reciprocal lattice vectors $\mathbf{G}$ and the periodic substrate with reciprocal lattice vectors $\mathbf{M}$ complicates things dramatically. The adsorbed over layer can be commensurate or incommensurate with the substrate, depending on the relationship of $\mathbf{G}$ and $\mathbf{M}$ [128, 129]. An early but comprehensive theoretical review on commensurate/incommensurate transitions is [163] and one on experiments on rare gases on graphite [203]. Graphite is one of the most widely used substrates because it can provide fairly large atomically flat surfaces of linear size  ∼1600 Å. Favorite adsorbates are rare gases such as Ar, Kr, Ne and Xe which are believed to interact by a simple Lennard-Jones interaction and differ only in size and monolayers of these have been extensively studied by diffraction of synchrotron radiation [204–206]. Xenon adsorbed on graphite has been extensively studied [207] where both first order melting of the incommensurate solid at lower coverage and continuous melting at higher coverage was seen. Krypton on graphite has been studied extensively [206, 208–214] and a tongue of fluid found separating the commensurate and incommensurate phases as shown in figure 26. High resolution scattering experiments are consistent with a smooth continuous F-IC transition [208] but the resolution is insufficient to make detailed fits to theory [215]. However, the scattering data seems to be inconsistent with KTHNY melting [208].

Zoom In Zoom Out Reset image size

Argon is very similar to Krypton except it has a slightly smaller atomic size and always forms an incommensurate solid overlayer on graphite with a phase diagram of figure 27 left. Diffraction studies of the melting transition around monolayer coverage are consistent with KTHNY melting but the correlation length is equally well fitted by a conventional power law $\xi (t)={{\xi}_{0}}{{t}^{-\nu}}$ [216]. It must be remembered that length scales are limited to $\xi =\mathcal{O}\left({{10}^{3}}\right)$ Å. These results agree with an early specific heat measurement [217] but disagree with a later specific heat measurement [218] and earlier diffraction experiments [204]. A later investigation [219] is also consistent with KTHNY melting and also with the prediction that the solid over layer is rotated from the substrate axes [220, 221] but again detailed comparison with theory is difficult because of insufficient experimental resolution. A study using high resolution vapor pressure isotherms [222] yields evidence in favor of continuous melting by a two step process which is constant with KTHNY melting. Although not conclusive, this study does resolve some discrepancies between the older measurements, which makes the KTHNY scenario more likely for monolayer Argon on graphite.

Zoom In Zoom Out Reset image size

Xenon adsorbed on graphite has a phase diagram of figure 27 right and has also been studied by diffraction [207, 219, 223–225] and by thermodynamic measurements [226–228]. The synchrotron diffraction experiments [207, 223, 224] are consistent with first order melting for low coverages which crosses over to KTHNY continuous melting at higher coverages. Thermodynamic measurements of the compressibility of the overlayer [227] are consistent with this. However, a high precision heat capacity and compressibility measurements [228] are not consistent with a crossover to continuous KTHNY melting and claim that melting is always first order. The scattering experiment of Rosenbaum et al [225] showed that melting at high coverages is well described by KTHNY melting in the presence of a weak ${{h}_{6}}\text{cos}6\theta$ substrate potential. A study of melting of Xenon adsorbed on Ag(1 1 1) [229] agrees with the KTHNY scenario. In fact, the 6-fold substrate potential of the silver substrate is extremely small [229] and this model fits the scattering data very well. The observed hexatic order of the fluid seems entirely due to the Xenon overlayer. The only effect of the substrate is to align the 6-fold hexatic axes by steps in the substrate [229] and the hexatic orientational order is intrinsic to the Xenon over layer. The remaining rare gas, Neon, on graphite has also been studied [230–232] but in much less detail than the others. In view of the recent simulations of the hard disk system [159] it seems that the conclusions from experiments must be viewed with some caution because of the very limited system sizes and the difficulty of distinguishing continuous from weak first order transitions.

Many other adsorbates on graphite have been studied by thermodynamic measurements with many complicated phase diagrams proposed. A neutron scattering study of melting of ethylene adsorbed on graphite has yielded evidence of a continuous transition [233] and compared to methane on the same substrate. Both films melt from incommensurate solids which rules out substrate heterogeneity being responsible for the apparent continuity in ethylene. Methane films show no signs of this broadening. Thus, the broadening of the transition is taken to be intrinsic to the overlayer melting by successive continuous transitions but the resolution does not permit verification that this is KTHNY melting [233].

Another relevant system in 2D is a freely suspended film of liquid crystal [96, 100, 132] on which measurements are possible. There are a number of theoretical papers generalizing earlier work on melting of arrays of point particles to 2D arrays of anisotropic particles [123, 202, 234]. Ostlund and Halperin [202] discuss dislocation mediated melting of 2D arrays of rods tilted from the normal to the layer. Note that, if the tilt angle is zero, the system can be described by an array of point particles so the theory is the same as the previous section. However, for tilted rods as in figure 28(a), the projection of the rods into the plane is a set of directed arrows arranged in a distorted triangular lattice. Two obvious physical cases considered are (i) the projection of the molecular axes tend to point along one of the six nearest neighbor bond directions as sketched in figure 28(b) and (ii) tend to point intermediate between two bonds as in figure 28(c). A very important consequence of the molecular tilt is that, even if in the absence of tilt, the molecules would form a regular hexagonal lattice, the molecules will form a uniaxial hexagonal lattice [123].

Zoom In Zoom Out Reset image size

To study the dislocation mediated melting of a 2D solid composed of rodlike molecules tilted at angle θ from the plane and the subsequent transitions which is relevant for freely suspended films of smectic liquid crystals [96, 100, 121, 131, 132], we must go through the same steps as for the melting of a hexagonal lattice. Unfortunately, even the first step is rather complicated because the energy of a set of dislocations in a uniaxial hexagonal lattice was not known in the late 1970's. The most explicit calculation is in [235] but, as pointed out by Ostlund and Halperin [202] the result depended on the choice of a cutoff and must be incorrect. In appendix A of [202] the dislocation energy is calculated without a cutoff and is more likely to be correct. In this review, I will try to focus on the salient points and leave all the gory technical details to the references [202].

The first complication is to compute the energy of a set of dislocations embedded in the elastic medium with the uniaxial symmetry of the underlying distorted hexagonal lattice shown in figure 29. In such a system, there are four different nonzero elastic constants or, equivalently, four distinct finite compliances. The compliance tensor S_ijkl is the inverse of the elastic tensor C_ijkl or ${{C}_{ijkl}}{{S}_{klmn}}=\frac{1}{2}\left({{\delta}_{im}}{{\delta}_{jn}}+{{\delta}_{in}}{{\delta}_{jm}}\right)$. The finite compliances are ${{S}_{1111}},{{S}_{2222}},{{S}_{1122}},{{S}_{1212}}$ where the axes are shown in figure 29.

$$\begin{eqnarray}&&H\left({{\mathbf{r}}_{1}}\cdots {{\mathbf{r}}_{n}}\right)={{H}_{\text{E}}}+{{H}_{\text{D}}},\\ &&\frac{{{H}_{\text{E}}}}{{{k}_{\text{B}}}T}=\frac{1}{2}\mathop{\int}^{}\frac{{{\text{d}}^{2}}\mathbf{r}}{a_{0}^{2}}{{\phi}_{ij}}{{C}_{ijkl}}{{\phi}_{kl}},\\ &&\frac{{{H}_{\text{D}}}}{{{k}_{\text{B}}}T}=\frac{1}{2}\underset{\alpha \ne \beta =1}{\overset{6}{\mathop \sum}}\,b_{i}^{\alpha}b_{j}^{\beta}{{E}_{ij}}\left({{\mathbf{r}}^{\alpha}}-{{\mathbf{r}}^{\beta}}\right)+\underset{\alpha}{\mathop \sum}\,E_{ij}^{c}b_{i}^{\alpha}b_{j}^{\alpha},\end{eqnarray} \tag{ 180 }$$

where $\alpha,\beta$ denote one of the six dislocations of figure 29 in an anisotropic hexagonal lattice, i, j  =  1, 2 are the two independent directions in 2D and ${{\mathbf{r}}^{\alpha}}$ is the position of the dislocation of type α. There are two independent Burger's vectors ${{\mathbf{b}}_{\text{I}}}$ and ${{\mathbf{b}}_{\text{II}}}$. The interaction energy of a pair of  ±  type α dislocations separated by $\mathbf{r}$ is ${{K}_{\alpha}}=b_{i}^{\alpha}b_{j}^{\alpha}{{E}_{ij}}\left(\mathbf{r}\right)$, where ${{E}_{ij}}\left(\mathbf{r}\right)={{K}_{ij}}\text{ln}(r/a)+{{V}_{ij}}(\theta )$ is the interaction energy which is a very, complicated expression given in appendix A of [202].

Zoom In Zoom Out Reset image size

The solid phase has dislocations which are tightly bound in pairs and when the temperature is raised, the pairs unbind and the crystalline order is destroyed. This melting occurs when $\Delta F=\Delta E-T\Delta S$ where $\Delta S={{k}_{\text{B}}}\text{ln}\left({{L}^{2}}/{{a}^{2}}\right)$ and $\Delta E=\pi K\text{ln}(L/a)$ where K is determined by the lattice spacing and elastic constants, L is the size of the system and a the lattice constant [19]. The melting temperature is given by $\Delta F=0$ or $\pi K=2{{k}_{\text{B}}}T$. There are two types of dislocation and two values ${{K}_{\text{I}}}$ and ${{K}_{\text{II}}}$ in an anisotropic solid, as shown in figure 29. Thus, melting of the solid will be controlled by the smaller of ${{K}_{\text{I}}},\,{{K}_{\text{II}}}$ which is, in turn, controlled by the shorter of the elementary Burger's vectors ${{\mathbf{b}}_{\text{I}}}$ and ${{\mathbf{b}}_{\text{II}}}$. Following the discussion of [123, 202, 234], one can study the melting of the anisotropic solid by the renormalization group procedure of section 3.1. Following [202], the compliance tensor S_ijkl is the inverse of the elastic tensor C_ijkl so that ${{C}_{ijkl}}{{S}_{klmn}}=\frac{1}{2}\left({{\delta}_{im}}{{\delta}_{jn}}+{{\delta}_{in}}{{\delta}_{jm}}\right)$ and, to $\mathcal{O}\left(\,{{y}^{2}}\right)$, and S_ijkl is given by

$$\begin{eqnarray}&&{{S}_{ijkl}}=S_{ijkl}^{0}+\frac{1}{4}\underset{\alpha =1}{\overset{3}{\mathop \sum}}\,y_{\alpha}^{2}\int_{0}^{2\pi}\text{d}\theta {{\text{e}}^{-{{V}_{\alpha}}(\theta )}}\left(b_{i}^{\alpha}{{\epsilon}_{jm}}{{\hat{r}}_{m}}+b_{j}^{\alpha}{{\epsilon}_{im}}{{\hat{r}}_{m}}\right) \\ &&\times \left(b_{k}^{\alpha}{{\epsilon}_{ln}}{{\hat{r}}_{n}}+b_{l}^{\alpha}{{\epsilon}_{kn}}{{\hat{r}}_{n}}\right)\int_{a}^{\infty}\frac{\text{d}r}{a}{{\left(\frac{r}{a}\right)}^{3-{{K}_{\alpha}}}}+\mathcal{O}\left(\,{{y}^{3}}\right),\end{eqnarray}$$

where the interaction energy of a  ±  pair of type α dislocations separated by $\mathbf{r}$ is

$$\begin{eqnarray}&&{{E}_{\alpha}}\left(\mathbf{r}\right)={{K}_{\alpha}}\text{ln}(r/a)+{{V}_{\alpha}}\left(\mathbf{\hat{r}}\right).\end{eqnarray}$$

The RG equations are derived from this by the standard method of rescaling $a\to a{{\text{e}}^{\delta l}}$ to obtain [202]

$$\begin{eqnarray}&&\frac{\text{d}{{S}_{1111}}}{\text{d}l}=y_{\text{I}}^{2}F_{1}^{22}+2y_{\text{II}}^{2}\text{co}{{\text{s}}^{2}}{{\phi}_{0}}F_{2}^{22},\\ &&\frac{\text{d}{{S}_{2222}}}{\text{d}l}=2y_{\text{II}}^{2}\text{si}{{\text{n}}^{2}}{{\phi}_{0}}F_{2}^{11},\\ &&\frac{\text{d}{{S}_{1122}}}{\text{d}l}=2y_{\text{II}}^{2}\text{sin}{{\phi}_{0}}\text{cos}{{\phi}_{0}}F_{2}^{12},\\ &&\frac{\text{d}{{S}_{1212}}}{\text{d}l}=\frac{1}{4}y_{\text{I}}^{2}F_{1}^{11}+\frac{1}{2}y_{\text{II}}^{2}\left(\text{co}{{\text{s}}^{2}}{{\phi}_{0}}+\text{si}{{\text{n}}^{2}}{{\phi}_{0}}F_{2}^{22}-2\text{sin}{{\phi}_{0}}\text{cos}{{\phi}_{0}}F_{2}^{12}\right),\\ &&\frac{\text{d}{{y}_{\text{I}}}}{\text{d}l}=\left(2-\frac{{{K}_{\text{I}}}}{2}\right){{y}_{\text{I}}}+2{{A}_{23}}y_{\text{II}}^{2},\\ &&\frac{\text{d}{{y}_{\text{II}}}}{\text{d}l}=\left(2-\frac{{{K}_{\text{II}}}}{2}\right){{y}_{\text{II}}}+2{{A}_{12}}{{y}_{\text{I}}}{{y}_{\text{II}}},\\ &&\text{where}\,\,F_{\alpha}^{ij}={{\left({{\mathbf{b}}^{\alpha}}\right)}^{2}}\int_{0}^{2\pi}\text{d}\theta {{\hat{r}}_{i}}{{\hat{r}}_{j}}{{\text{e}}^{-{{V}_{\alpha}}(\theta )}},\,\,\,{{A}_{\alpha \beta}}=\int_{0}^{2\pi}\text{d}\theta {{\text{e}}^{b_{i}^{\alpha}b_{j}^{\beta}{{V}_{ij}}(\theta )}}.\end{eqnarray} \tag{ 181 }$$

An expression for ${{V}_{ij}}(\theta )$ is given in appendix A of [202]. The quantities $F_{\alpha}^{ij}$ and ${{A}_{\alpha \beta}}$ depend on the compliances which have finite fixed point values and differ from these by $\mathcal{O}\left(\,{{y}^{2}}\right)$. Thus, one can consider these as l-independent constants [202].

For type I melting, ${{K}_{\text{II}}}>4$ and ${{K}_{\text{I}}}\approx 4$ so that the fugacity ${{y}_{\text{II}}}\to 0$ because $\left(\frac{{{K}_{\text{II}}}}{2}-2\right)>0$ and the flow equations involving ${{y}_{\text{I}}}$ become

$$\begin{eqnarray}&&\frac{\text{d}{{S}_{1111}}}{\text{d}l}=y_{\text{I}}^{2}F_{1}^{22},\\ &&\frac{\text{d}{{S}_{1212}}}{\text{d}l}=\frac{1}{4}y_{\text{I}}^{2}F_{1}^{11},\\ &&\frac{\text{d}{{y}_{\text{I}}}}{\text{d}l}=\frac{1}{2}\left(4-{{K}_{\text{I}}}\right){{y}_{\text{I}}}.\end{eqnarray} \tag{ 182 }$$

Near the melting temperature, the compliances S₁₁₂₂ and S₂₂₂₂ are dominated by analytic terms since their RG flows are not controlled by the fugacity ${{y}_{\text{I}}}$. The compliances ${{S}_{1111}},{{S}_{1212}}$ are expanded as

$$\begin{eqnarray}&&{{S}_{1111}}(l)=S_{1111}^{\ast}+m\left(\,{{y}_{\text{I}}}(l)+D(l)\right).\end{eqnarray} \tag{ 183 }$$

Here, $S_{1111}^{\ast}$ is the fixed point value of S₁₁₁₁ at melting and is equal to its renormalized value at ${{T}_{\text{m}}}$. D(l) is the deviation from the incoming separatrix with $D(0)=\,|t|$ and ${{y}_{\text{I}}}$ is the value of the fugacity on the incoming separatrix.

Below melting, $T\leqslant {{T}_{\text{m}}}$, after some further algebra [202], one can derive the flow equations for ${{y}_{\text{I}}}(l)$ and D(l)

$$\begin{eqnarray}&&\frac{\text{d}D}{\text{d}l}=-\frac{F_{1}^{22}}{{{m}_{-}}}{{y}_{\text{I}}}D,\\ &&\frac{\text{d}{{y}_{\text{I}}}}{\text{d}l}=\frac{{{F}_{1}}}{{{m}_{-}}}y_{\text{I}}^{2}\end{eqnarray} \tag{ 184 }$$

where the slope m₋  <  0 [202]. For $|t|\ll 1$, ${{y}_{\text{I}}}(l)\sim {{l}^{-1}}$ for $l\gg 0$. Also $D(l){{y}_{\text{I}}}(l)=D(0){{y}_{\text{I}}}(0)$ which allows one to define a value l^* such that $D\left({{l}^{\ast}}\right)={{y}_{\text{I}}}\left({{l}^{\ast}}\right)$ such that for l  >  l^*, the fugacity ${{y}_{\text{I}}}(l)\to 0$ so that one can write $|t|\propto D(0)\propto {{\left({{l}^{\ast}}\right)}^{2}}$. In turn, this observation, allows one to define a length scale ${{\xi}_{-}}(t)$ and the singularities in the compliances

$$\begin{eqnarray}&&{{\xi}_{-}}(t)=\text{exp}\left(|t{{|}^{-1/2}}\right),\\ &&S_{ijkl}^{\text{R}}(t)=S_{ijkl}^{\ast}\left(1+b\sqrt{|t|}\right),\,\,\,\,\,\,\,\,(ijkl)=(1111),(1212),\\ &&S_{ijkl}^{\text{R}}(t)=S_{ijkl}^{\ast}\left(1+{{b}^{\prime}}\xi _{-}^{-2p}\right)\,\,\,2p={{K}_{\text{II}}}-4,\,\,\,(ijkl)=(1122),(2222).\end{eqnarray} \tag{ 185 }$$

where b  >  0 and ${{b}^{\prime}}>0$ are nonuniversal constants $\mathcal{O}(1)$. Thus for $T\leqslant {{T}_{\text{m}}}$, the compliances S₁₁₁₁ and S₁₂₁₂ have $\sqrt{|t|}$ singularities while S₁₁₂₂ and S₂₂₂₂ have undetectable essential singularities $\xi _{S}^{-2p}$ at ${{T}_{\text{m}}}$, since ${{K}_{\text{II}}}>4$ [202].

Above melting, $T>{{T}_{\text{m}}}$, the RG trajectories for ${{y}_{\text{I}}}(l)$ follow the incoming separatrix up to $l=l_{1}^{\ast}$, break away and join the outgoing separatrix ${{y}_{\text{I}}}(l)={{m}_{+}}\left({{S}_{1111}}(l)-S_{1111}^{\ast}\right)$ for another $l_{2}^{\ast}$ RG iterations. ${{y}_{\text{I}}}\left(l_{2}^{\ast}\right)=\mathcal{O}(1)$ but is still small. Thus, we see that the length scale ${{\xi}_{S}}\sim \text{exp}\left(l_{1}^{\ast}+l_{2}^{\ast}\right)$ is the typical separation of free dislocations of type I. The fugacity ${{y}_{\text{II}}}(l)$ is decreasing exponentially for $l\leqslant l_{1}^{\ast}+l_{2}^{\ast}$. Thus, we can integrate out the type I dislocations by a Debye–Hückel approximation to obtain an effective Hamiltonian for the type II dislocations at scale $l_{S}^{\ast}=l_{1}^{\ast}+l_{2}^{\ast}$, which is the Hamiltonian for a 2D smectic liquid crystal [202, 234]

$$\begin{eqnarray}&&\frac{H_{\text{D}}^{S}}{{{k}_{\text{B}}}T}=\frac{B}{2}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{k}\left(\frac{{{\left(\lambda {{k}_{x}}\right)}^{2}}}{k_{y}^{2}+{{\lambda}^{2}}k_{x}^{4}}+\frac{2{{E}_{s}}\xi _{S}^{2}}{B}\right)|{{b}_{y}}\left(\mathbf{k}\right){{|}^{2}},\\ &&B={{\left(S_{2222}^{\ast}\right)}^{-1}},\,\,\,\,{{\lambda}^{2}}=2\xi _{S}^{2}E_{11}^{c}S_{2222}^{\ast}=2\xi _{S}^{2}S_{2222}^{\ast},\\ &&{{E}_{s}}=E_{22}^{c}\left(l_{S}^{\ast}\right)=-\text{ln}{{y}_{\text{II}}}\left(l_{S}^{\ast}\right)=p\left(l_{1}^{\ast}+l_{2}^{\ast}\right).\end{eqnarray} \tag{ 186 }$$

At first sight, it seems that at this length scale one can treat both type I and type II dislocations by a Debye–Hückel approximation but ${{y}_{\text{II}}}\left(l_{S}^{\ast}\right)\ll 1$ so that the discrete nature of ${{\mathbf{b}}_{\text{II}}}$ must still be accounted for [234]. Thus, a further $l_{3}^{\ast}$ iterations are performed so that

$$\begin{eqnarray}&&{{y}_{\text{II}}}\left(l_{1}^{\ast}+l_{2}^{\ast}+l_{3}^{\ast}\right)={{y}_{\text{II}}}(0)\,\,\Rightarrow \,\,{{\text{e}}^{-p\left(l_{1}^{\ast}+l_{2}^{\ast}\right)+l_{3}^{\ast}}}=1.\end{eqnarray} \tag{ 187 }$$

At this length scale, the system becomes a gas of free dislocations of all types and one can use a Debye–Hückel approximation. At this scale, the Frank constants

$$\begin{eqnarray}&&{{K}_{x}}\sim \xi _{S}^{2}\sim \xi _{S}^{2}\sim \text{exp}\left(2{{t}^{-1/2}}\right),\\ &&{{K}_{y}}\sim \xi _{N}^{2}\sim \xi _{S}^{2(\,p+1)}.\end{eqnarray} \tag{ 188 }$$

At length scale L with $L<{{\xi}_{S}}(t)$, all dislocations are bound and the system responds like a solid while, if ${{\xi}_{S}}<L<{{\xi}_{N}}$, it has a smectic like response. At still larger scales, ${{\xi}_{N}}<L<{{\xi}_{I}}$, it behaves like a nematic with two Frank constants ${{K}_{y}}/{{K}_{x}}\propto \xi _{S}^{2p}$. For still larger scales, higher order terms like ${{\left(\nabla \centerdot \mathbf{n}\right)}^{4}}$ in the Hamiltonian will renormalize ${{K}_{x}},{{K}_{y}}$ to equality for $L>{{\xi}_{I}}\sim \text{exp}\left(c\xi _{S}^{2}\right)$ [236]. This melting scenario is called type I melting [202] with a phase diagram sketched in figure 27, left. When ${{K}_{\text{II}}}>{{K}_{\text{I}}}$, the smectic-like phase is absent and the phase diagram is sketched in figure 27, right.

Another system which, at first sight, seems that it should be a good example of experimentally accessible Kosterlitz–Thouless physics is a thin film of superconductor. In this system, it is clear that the physics will be controlled by vortices in the phase of the superconducting order parameter $\Psi\left(\mathbf{r}\right)=\,|\Psi\left(\mathbf{r}\right)|{{\text{e}}^{\text{i}\phi \left(\mathbf{r}\right)}}$ where ϕ is the phase of the condensate wavefunction. This is exactly analogous to the superfluid ⁴He situation. The main difference is that the superconducting condensate is charged which implies that the all important vortex–vortex interaction is screened from the $\text{ln}r$ form of the neutral superfluid to 1/r due to screening currents flowing round the vortices [237, 238].

This is best discussed from the Ginzburg–Landau form of the free energy for a charged condensate of Cooper pairs [239]

$$\begin{eqnarray}&&\mathcal{F}\,[\Psi]=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\text{d}z\left\{\frac{1}{2{{m}^{\ast}}}|\left(-\text{i}\hslash \nabla -\frac{{{e}^{\ast}}}{c}\mathbf{A}\right)\Psi{{|}^{2}}\right. \\ &&\left. +\frac{1}{2}r(T)|\Psi{{|}^{2}}+\frac{1}{4}u|\Psi{{|}^{4}}+\frac{{{B}^{2}}}{8\pi}-\frac{\mathbf{H}\centerdot \mathbf{B}}{4\pi}\right\}.\end{eqnarray} \tag{ 189 }$$

$\Psi\left(\mathbf{r}\right)$ is the coarse grained Cooper pair wave function which is non zero only in the plane of the film z  =  0, $\mathbf{A}\left(\mathbf{r},z\right)$ is the vector potential which is finite outside the plane of the film and $\mathbf{B}\left(\mathbf{r},z\right)=\nabla \times \mathbf{A}$. The parameter r(T)  <  0 for T  <  T₀ the BCS temperature for the onset of the formation of Cooper pairs and u is positive and T independent. Here $\mathbf{r}=(x,y)$ is a point in the film in the z  =  0 plane, the Cooper pair mass ${{m}^{\ast}}=2{{m}_{e}}$ and e^*  =  2e. The essential difference between the neutral ⁴He superfluid and the charged superfluid of a superconductor is that a vortex in a superconductor has a circulating electric current producing an associated magnetic field $\mathbf{B}\left(\mathbf{r},z\right)$ which extends outside the film at $\left(\mathbf{r},0\right)$.

The behavior of a superconducting film was first discussed by Pearl [237] who found that the circulating current of a vortex at $\mathbf{r}=0$ falls off as 1/r² instead of exponentially as in a bulk superconductor. Also, the screening length in a film is ${{\lambda}_{\text{eff}}}=\lambda _{0}^{2}/d$, where d is the film thickness and ${{\lambda}_{0}}=\sqrt{{{m}^{\ast}}{{c}^{2}}/\left(4\pi n_{s}^{\ast}{{\left({{e}^{\ast}}\right)}^{2}}\right)}=\mathcal{O}\left({{10}^{4}}~\overset{\circ}{\mathop{\text{A}}}\,\right)$ is the London penetration depth. The Ginzburg–Landau coherence length is

$$\begin{eqnarray}&&{{\xi}_{0}}=\frac{\hslash}{{{m}^{\ast}}}\frac{1}{\sqrt{2|r(T)|}}\ll d\end{eqnarray} \tag{ 190 }$$

which is the length scale over which $|\Psi\left(\mathbf{r}\right)|$ varies so that $\Psi\left(\mathbf{r},z\right)$ is constant over the film thickness. Since one is interested in the temperature range over which phase fluctuations are important, it is clear that the parameter $r(T)\ll 0$ so that $|\Psi\left(\mathbf{r}\right)|\,=\mathcal{O}(1)$ since the phase is not defined otherwise. Thus, one can assume that the temperature $T\ll {{T}_{0}}$, the Ginzburg–Landau critical temperature. The system is effectively two dimensional when the film thickness $d\ll {{\xi}_{0}}$ since we can integrate (189) over z to obtain an effective 2D theory with $|{{\Psi}_{0}}{{|}^{2}}=n_{s}^{0}$, the Cooper pair areal number density, which is taken to be constant except at vortex cores. The vortex core has radius $\mathcal{O}\left({{\xi}_{0}}\right)$ which is the smallest length scale of the system.

When $\mathbf{H}=0$, these considerations lead to the free energy $\mathcal{F}\left(\theta,\mathbf{A}\right)$

$$\begin{eqnarray}&&\mathcal{F}\left(\theta,\mathbf{A}\right)=\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\text{d}z\left\{\frac{n_{s}^{0}}{2{{m}^{\ast}}}{{\left(\hslash \nabla \theta -\frac{{{e}^{\ast}}}{c}\mathbf{A}\right)}^{2}}+\frac{1}{8\pi}{{\left(\nabla \times \mathbf{A}\right)}^{2}}\right\}\end{eqnarray} \tag{ 191 }$$

From the extremal equation $\delta \mathcal{F}/\delta \mathbf{A}=0$, we obtain

$$\begin{eqnarray}&&\nabla \times \left(\nabla \times \mathbf{A}\right)=\nabla \times \mathbf{B}=\frac{4\pi \mathbf{J}}{c}=4\pi |{{\Psi}_{0}}{{|}^{2}}\frac{{{e}^{\ast}}\hslash}{{{m}^{\ast}}c}\left(\nabla \theta -\frac{{{e}^{\ast}}}{\hslash c}\mathbf{A}\right),\\ &&\Rightarrow \mathbf{J}\left(\mathbf{r},z\right)={{e}^{\ast}}|{{\Psi}_{0}}{{|}^{2}}\frac{\hslash}{{{m}^{\ast}}}\left(\nabla \theta -\frac{{{e}^{\ast}}}{\hslash c}\mathbf{A}\right),\end{eqnarray} \tag{ 192 }$$

where e^*  =  2e and ${{m}^{\ast}}=2{{m}_{e}}$ are the charge and mass of a Cooper pair. To see the effect of the field $\mathbf{B}\left(\mathbf{r},z\right)$ due to a set of vortices in the superconducting plane at z  =  0, one follows the analysis of Pearl [237] and Nelson [238] by assuming the superconducting current is only in the z  =  0 plane so that

$$\begin{eqnarray}&&\mathbf{J}\left(\mathbf{r},z\right)=\frac{c{{\phi}_{0}}}{8{{\pi}^{2}}{{\lambda}_{\bot}}}\left(\nabla \theta \left(\mathbf{r}\right)-\frac{2\pi}{{{\phi}_{0}}}{{\mathbf{A}}_{2d}}\left(\mathbf{r}\right)\right)\delta (z),\end{eqnarray} \tag{ 193 }$$

where the flux quantum ${{\phi}_{0}}=\pi \hslash c/e$, ${{\mathbf{A}}_{2d}}\left(\mathbf{r}\right)=\mathbf{A}\left(\mathbf{r},0\right)$, and the transverse penetration depth ${{\lambda}_{\bot}}={{m}_{e}}{{c}^{2}}/\left(8\pi {{n}_{s}}{{e}^{2}}\right)$. Now one can use the standard equations of magnetostatics to solve for the current due to a vortex of unit strength at $\mathbf{r}=0$ so that

$$\begin{eqnarray}&&-{{\nabla}^{2}}\mathbf{A}\left(\mathbf{r},z\right)+\frac{1}{{{\lambda}_{\bot}}}{{\mathbf{A}}_{2d}}\left(\mathbf{r}\right)=\frac{{{\phi}_{0}}}{2\pi {{\lambda}_{\bot}}}\frac{\mathbf{\hat{z}}\times \mathbf{r}}{{{r}^{2}}}\delta (z),\end{eqnarray} \tag{ 194 }$$

since $\nabla \theta \left(\mathbf{r}\right)=2\pi \left(\mathbf{\hat{z}}\times \mathbf{r}\right)/{{r}^{2}}$ due to the vortex at $\mathbf{r}=0$ and ${{\mathbf{A}}_{2d}}\left(\mathbf{r}\right)\equiv \mathbf{A}\left(\mathbf{r},z=0\right)$. This can be solved by taking Fourier transforms, with the result [237, 238]

$$\begin{eqnarray}&&{{\mathbf{J}}_{2\text{d}}}(r,\theta )=\frac{c{{\phi}_{0}}}{8{{\pi}^{2}}\text{d}{{\lambda}_{\bot}}r}\widehat{\boldsymbol{\theta }},\,\,\,\,r\ll {{\lambda}_{\bot}},\\ &&=\frac{c{{\phi}_{0}}}{4{{\pi}^{2}}\text{d}{{r}^{2}}}\widehat{\boldsymbol{\theta }},\,\,\,\,\,\,\,\,\,\,r\gg {{\lambda}_{\bot}},\end{eqnarray} \tag{ 195 }$$

which translates into the vortex–vortex interaction behaving as $\text{ln}r$ for $r<{{\lambda}_{\bot}}$ and as 1/r² for $r\gg {{\lambda}_{\bot}}$.

It is very important to note that the onset of superconductivity in 2D is due to vortex type phase fluctuations which requires a pre-existing condensate or pairing of electrons. This phenomenon is described by a Ginzburg–Landau free energy functional of (189) in which the condensate forms at the mean field temperature T₀ when r(T₀)  =  0. Even at T  <  T₀, although there is a finite condensate or BCS pairing, $|{{\Psi}_{0}}{{|}^{2}}>0$, there is not necessarily any superconductivity. This requires phase coherence or a finite renormalized superfluid density, just as in the neutral superfluid ⁴He which is controlled by vortices. The onset of phase coherence or true superconductivity in the thermodynamic limit of infinite system size $L\to \infty$ occurs when the experimentally measurable kinetic inductance ${{L}_{K\Box }}(T)={{m}_{e}}/\left(2{{n}_{s}}{{e}^{2}}\right)$ which is related to the transverse penetration depth

$$\begin{eqnarray}&&{{\lambda}_{\bot}}(T)=\frac{{{m}_{e}}{{c}^{2}}}{8\pi {{e}^{2}}{{n}_{s}}}=\frac{{{L}_{K\Box }}(T){{c}^{2}}}{4\pi}.\end{eqnarray} \tag{ 196 }$$

In the absence of an external magnetic field $\mathbf{H}=0$, the condition for the onset of superconductivity is given by the analogue of (46), [240, 241]

$$\begin{eqnarray}&&{{k}_{\text{B}}}{{T}_{\text{c}}}{{\lambda}_{\bot}}\left({{T}_{\text{c}}}\right)=\frac{\phi _{0}^{2}}{32{{\pi}^{2}}},\end{eqnarray} \tag{ 197 }$$

where ${{T}_{\text{c}}}$ is the temperature at which phase coherence first appears

From this discussion, we see that there can be no true phase transition in a superconducting film because of the finite penetration depth ${{\lambda}_{\bot}}(T)$ for $T\leqslant {{T}_{\text{c}}}$ [19]. However, since ${{\lambda}_{\bot}}(T)$ can be $\mathcal{O}\left({{10}^{-2}}\right)$m which is larger than most experimental systems so that we can regard most superconducting films at finite temperature when quantum effects are unimportant as described by the class of Hamiltonians

$$\begin{eqnarray}&&\frac{\mathcal{F}\left(\theta,\mathbf{A}\right)}{{{k}_{\text{B}}}T}={{K}_{0}}(T)\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}\left[1-\text{cos}\left(\nabla \theta \left(\mathbf{r}\right)-\frac{2e}{\hslash c}\mathbf{A}\left(\mathbf{r}\right)\right)\right]\\ &&=\frac{{{K}_{0}}(T)}{2}\mathop{\int}^{}{{\text{d}}^{2}}\mathbf{r}{{\left(\nabla \theta \left(\mathbf{r}\right)-\frac{2e}{\hslash c}\mathbf{A}\left(\mathbf{r}\right)\right)}^{2}},\\ &&{{K}_{0}}(T)=\frac{{{\hslash}^{2}}n_{s}^{0}(T)}{2m{{k}_{\text{B}}}T},\end{eqnarray} \tag{ 198 }$$

for a uniform 2D superconducting film, where $n_{s}^{0}(T)$ is the bare superfluid number density and 2m the Cooper pair mass. In a clean uniform film, the mean field Cooper pair formation temperature T₀ and the phase coherence onset temperature ${{T}_{\text{c}}}<{{T}_{0}}$ are rather close and effects due to these tend to become mixed together. It is known that Cooper pairing alone does not imply superconductivity which requires global phase coherence. Thus, to study superconductivity in a thin film, it is important that T₀ and ${{T}_{\text{c}}}$ are well separated. This can be accomplished by a granular film with a large resistance R_N in the normal state [242–247] or by coupling the grains by weak links [248] as in arrays of Josephson junctions [249].

When the external magnetic field H_z  =  0, the superconducting film should be like a finite size 2D planar rotor model with $L={{\lambda}_{\bot}}$. It was pointed out [240, 241] that the vortex unbinding mechanism [19] could be very important in thin superconducting films provided the transverse penetration depth ${{\lambda}_{\bot}}$ was very large which could be realized in very thin films with a large normal sheet resistance ${{\lambda}_{\bot}}\propto {{R}_{\Box }}$. When ${{R}_{\Box }}>{{10}^{3}}\Omega$, there is a reasonable separation between the pair formation and KT temperatures ${{T}_{0}}>{{T}_{\text{KT}}}$ which is essential for the experimental observation of vortex unbinding. They proposed [240] that ${{T}_{\text{KT}}}\approx 2.18{{T}_{0}}{{R}_{\text{c}}}/{{R}_{\Box }}$ but, on a more fundamental level [249, 250],

$$\begin{eqnarray}&&{{k}_{\text{B}}}{{T}_{\text{KT}}}\approx \frac{\pi {{\hslash}^{2}}}{8{{e}^{2}}{{L}_{\Box }}}=\frac{{{i}_{\text{c}}}{{\phi}_{0}}}{4}=\frac{\pi {{E}_{\text{J}}}}{2},\end{eqnarray} \tag{ 199 }$$

where the vortex fugacity is ignored, which fits the data much better.

In a uniform superconducting film in a finite magnetic field B normal to the plane of the film there will be a hexagonal lattice of vortices each containing one flux quantum ${{\phi}_{0}}=\hslash c/2e$ which form a stable triangular lattice of lattice spacing a₀ where [251, 252]

$$\begin{eqnarray}&&\frac{\sqrt{3}a_{0}^{2}}{2}=\frac{1}{n}=\frac{{{\phi}_{0}}}{B},\end{eqnarray} \tag{ 200 }$$

where n is the areal vortex density due to the field. In the presence of arbitrarily weak pinning of vortices, the flux lattice will be pinned for $T<{{T}_{\text{m}}}$, the lattice melting temperature, and vortices will be free to move when $T>{{T}_{\text{m}}}$ so that the onset of superconductivity can be identified with the vortex lattice melting temperature ${{T}_{\text{m}}}$ [252].

$$\begin{eqnarray}&&\frac{a_{0}^{2}\mu \left({{T}_{\text{m}}}\right)\left[\mu \left({{T}_{\text{m}}}\right)+\lambda \left({{T}_{\text{m}}}\right)\right]}{{{k}_{\text{B}}}{{T}_{\text{m}}}\left[2\mu \left({{T}_{\text{m}}}\right)+\lambda \left({{T}_{\text{m}}}\right)\right]}=4\pi,\end{eqnarray} \tag{ 201 }$$

and $\mu (T)=0$ when $T>{{T}_{\text{m}}}$. In (201), the elastic constants should be interpreted as the renormalized ${{\mu}^{\text{R}}}(T)$ and ${{\lambda}^{\text{R}}}(T)$. The phase diagram of a 2D superconductor in the (H, T) plane is sketched in figure 31 assuming that the transition is caused by vortex lattice melting [252].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the most important field regime, $\xi _{\text{GL}}^{-2}\gg B/{{\phi}_{0}}\gg \lambda _{\bot}^{-2}$, the bare ${{\lambda}_{0}}=\infty$ because of the long range vortex–vortex interactions and the bare shear modulus [253]

$$\begin{eqnarray}&&{{\mu}_{0}}(T)=\frac{{{\phi}_{0}}B}{32{{\pi}^{2}}{{\lambda}_{\bot}}(T)},\\ &&{{k}_{\text{B}}}{{T}_{\text{m}}}=C\frac{\phi _{0}^{2}}{64\sqrt{3}{{\pi}^{2}}{{\lambda}_{\bot}}\left({{T}_{\text{m}}}\right)},\end{eqnarray} \tag{ 202 }$$

where $0.4\leqslant C\leqslant 0.75$ [252]. When $T>{{T}_{\text{m}}}$, theory says that there should be a hexatic fluid up to T_i but this is not relevant as there is no known probe coupling to a vortex hexatic. Of course, in any real system, there will be pinning of vortices by some random pinning sites [252] which will not be discussed in this review. However, it should be noted that the assumed KTHNY melting of the vortex lattice is controversial. Even the existence of a 2D flux lattice has been questioned [254]. Later work [255–261] showed that, although numerical estimates of ${{\mu}_{\text{R}}}(T)$ are consistent with the KTHNY scenario, vortex lattice melting proceeds by a weakly first order transition so that ${{\mu}_{\text{R}}}(T)$ is very close to the KTHNY value.

There has been much interest in trying to apply the theory of 2D vortex lattice melting to thin films of YBCCO in a transverse magnetic field. This has met with only partial success [262–264] as vortex lattice melting does not appear to be the mechanism for the destruction of superconductivity [265, 266]. Vortex pinning by random impurities dominates [267] so that the system is better described as a moving glass phase of a current driven moving elastic lattice [268–271]. However, the basic ideas of dislocations and disclinations in a 2D vortex lattice have been verified by direct observation [272].

A regular array of superconducting grains coupled by almost identical SNS or SIS junctions can be made with up to about 10⁶ grains with modern lithographic deposition techniques [273] and almost any desired geometry and coupling constant arrangement can be constructed. Also, when an external magnetic field is applied, a frustrated XY model can be obtained, described by Coulomb gas Hamiltonian [274]

$$\begin{eqnarray}&&H=\frac{1}{2}\underset{ij}{\mathop \sum}\,\left({{q}_{i}}-{{Q}_{i}}\right)C_{ij}^{-1}\left({{q}_{j}}-{{Q}_{j}}\right)+{{E}_{\text{J}}}\left(1-\text{cos}\left({{\theta}_{1}}-{{\theta}_{j}}-{{A}_{ij}}\right)\right),\end{eqnarray} \tag{ 203 }$$

where the first term is the charging energy with C_ij the capacitance matrix, ${{q}_{i}}=2e{{n}_{i}}$ with n_i the number of Cooper pairs on grain i, and ${{Q}_{i}}=2e{{q}_{\text{ex}}}=\sum_{j}{{C}_{ij}}V_{j}^{\text{ex}}$ is the charge on grain i induced by the external gate voltage. The capacitance matrix C_ij has diagonal elements ${{C}_{ii}}={{C}_{0}}+zC$ where z is the number of nearest neighbors and C_ij  =  −C for i, j nearest neighbors as shown in figure 32. ${{E}_{\text{J}}}$ is the Josephson coupling energy between adjacent grains and ${{A}_{ij}}=(2e/\hslash c)\int_{i}^{j}\mathbf{A}\centerdot \text{d}\mathbf{l}$. The magnetic frustration is the flux f penetrating a plaquette given by $f=(1/2\pi )\sum_{\text{P}}{{A}_{ij}}$ where $\sum_{\text{P}}$ means the sum over the bonds round the plaquette in a counter clockwise direction. Quantum mechanical effects result from the non commuting operators $\left[{{\theta}_{i}},{{n}_{j}}\right]=\text{i}{{\delta}_{ij}}$ which become important at $T=\mathcal{O}{{\left({{10}^{-3}}\right)}^{{}^\circ}}K$ when the capacitances C₀, C are small so the grain charging energy ${{E}_{\text{c}}}$ is comparable to the Josephson energy ${{E}_{\text{J}}}$.

Zoom In Zoom Out Reset image size

Many experiments are done on arrays of large capacity superconducting grains on a lattice coupled either by the proximity effect [248] or arrays of Josephson junctions with ${{E}_{\text{C}}}\ll {{E}_{\text{J}}}$ [249]. These are best described by a lattice model of the form

$$\begin{eqnarray}&&\frac{\mathcal{F}(\theta,A)}{{{k}_{\text{B}}}T}=\frac{{{K}_{0}}(T)}{2}\underset{\langle ij\rangle}{\mathop \sum}\,\left[1-\text{cos}\left({{\theta}_{i}}-{{\theta}_{j}}-\frac{2\pi}{{{\phi}_{0}}}{{A}_{ij}}\right)\right],\end{eqnarray} \tag{ 204 }$$

where ${{A}_{ij}}=\int_{{{\mathbf{r}}_{i}}}^{{{\mathbf{r}}_{j}}}\text{d}\mathbf{r}\centerdot \mathbf{A}\left(\mathbf{r}\right)$ and charging effects can be ignored.

A junction array in a uniform magnetic field is a realization of an interesting class of models which, in the Coulomb gas representation is

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=\pi {{K}_{0}}(T)\underset{i,j}{\mathop \sum}\,\left({{n}_{i}}+{{f}_{i}}\right){{G}_{ij}}\left({{n}_{j}}+{{f}_{j}}\right)\end{eqnarray} \tag{ 205 }$$

where i labels the center of the ith plaquette on the dual lattice, ${{n}_{i}}=0,\pm 1,\cdots$ and $0\leqslant {{f}_{i}}\leqslant 1/2$ the magnetic flux ${{f}_{i}}=2\pi /{{\phi}_{0}}\mathop{\oint}^{}\text{d}\mathbf{l}\centerdot \mathbf{A}$ penetrating the ith plaquette. The interaction G_ij is the screened Coulomb interaction, ${{G}_{ij}}\sim \text{ln}\left(|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}|/a\right)$ for $|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}|<{{\lambda}_{\bot}}$ and ${{G}_{ij}}\sim \left({{a}^{2}}/|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}{{|}^{2}}\right)$ when $|{{\mathbf{R}}_{i}}-{{\mathbf{R}}_{j}}|>{{\lambda}_{\bot}}$. Even in this classical limit, these systems provide realizations of a great variety of models in 2D and have been the subject of intense activity in the 1980's. The models become even more interesting, and difficult, at very low T when quantum effects play a vital role [274–282].

The systems with a Hamiltonian of (205) represent a large class of statistical mechanical models. When the frustration f_i  =  0 in zero applied magnetic field, the system on a square lattice becomes a 2D planar rotor model when the screening length ${{\lambda}_{\bot}}\to \infty$. In a junction array one can make this larger than the array size, so the junction array is a realization of the planar rotor model. In the phase representation the Hamiltonian of an array is [283],

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}={{K}_{0}}\underset{<ij>}{\mathop \sum}\,\left[1-\text{cos}\left({{\theta}_{i}}-{{\theta}_{j}}-{{A}_{ij}}\right)\right],\\ &&{{A}_{ij}}+{{A}_{jk}}+{{A}_{kl}}+{{A}_{li}}=2\pi f\left(\mathbf{R}\right).\end{eqnarray} \tag{ 206 }$$

Because of the arbitrariness in the vector potential, for the fully frustrated case of f  =  1/2, one can choose A_ij  =  0 on horizontal bonds and ${{A}_{ij}}=\pi$ on alternating vertical bonds. The system is mapped into one with ferromagnetic bonds of strength J on horizontal bonds and alternating vertical strength J ferromagnetic and strength $\eta J$ antiferromagnetic bonds [284]. The phase diagram in the $\eta,T$ plane from MC simulations shows that, away from the isotropic point $\eta =1$, there is an Ising transition followed by an XY transition at higher T. The situation near the isotropic point $\eta =1$ is unclear but they seemed to occur simultaneously [284].

About the same time as this investigation was being done, there was much interest in the isotropic fully frustrated array with f  =  1/2 [283–303]. It was well understood that the isotropic fully frustrated system has both chiral order and phase coherence at sufficiently low temperature and that phase coherence required chiral, Ising, long range order, but not vice-versa. Thus the main question was whether the XY transition happened at a lower temperature than the Ising transition or whether the two transitions happened simultaneously and, if so, what is the nature of this transition? The problem is analytically intractable and there have been many simulations over the years, mostly inconclusive with a few exceptions.

The problem was discussed in [283] where it was explained that the renormalized helicity modulus in the XY ordered phase must obey $\Upsilon (T)/{{k}_{\text{B}}}T\geqslant 2/\pi$ and will drop discontinuously to zero at $T={{T}_{\text{KT}}}$ with a jump to zero $\Delta \Upsilon /{{k}_{\text{B}}}{{T}_{\text{KT}}}\geqslant 2/\pi$. There are two possibilities: (i) $\Delta \Upsilon /{{k}_{\text{B}}}{{T}_{\text{KT}}}=2/\pi$ at some ${{T}_{\text{KT}}}<{{T}_{I}}$ or (ii) as $T\to T_{I}^{-}$, there is a non-universal jump $\Delta \Upsilon /{{k}_{\text{B}}}T>2/\pi$. At first sight, these alternatives could be readily distinguished by simulation, but this hope was rapidly dashed because either the system sizes were too small to resolve the two transitions or they could be resolved with ${{T}_{\text{KT}}}<{{T}_{I}}$ but the jump in $\Upsilon (T)$ at ${{T}_{\text{KT}}}$ we larger than the universal value of $2/\pi$. This situation led to the suggestion that the XY-Ising model [287] would be easier to simulate where the Hamiltonian is

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=-\underset{<ij>}{\mathop \sum}\,\left[\left(A+B{{s}_{i}}{{s}_{j}}\right)\text{cos}\left({{\theta}_{i}}-{{\theta}_{j}}\right)+C{{s}_{i}}{{s}_{j}}\right],\end{eqnarray} \tag{ 207 }$$

where ${{s}_{i}}=\pm 1$ represents the chirality of a plaquette. However, in the end it turns out that despite the short range interactions of the XY-Ising model, in the interesting region of $C\leqslant 0$ and A  =  B, simulations of the model of (207) turn out to be just as difficult as the original representation.

The derivation of (207) for frustrated systems was done in a series of papers [287, 288, 304–310] by a set of approximate steps based on symmetry and renormalization group arguments and missed some of the short distance physics. For example, for an isotropic fully frustrated junction array [309], A  =  B so that there is no phase stiffness across an Ising domain wall with ${{s}_{i}}{{s}_{j}}=-1$ when i, j represent sites on opposite sides of domains of opposite chirality. In the original system, described by (206) with f_i  =  f  =  1/2, there is a finite phase stiffness across domain walls. However, in an important paper, Korshunov [311] showed that domain wall kinks unbind at ${{T}_{\text{K}}}<{{T}_{\text{V}}}$, the vortex unbinding or phase coherence temperature. It was also shown that ${{T}_{\text{V}}}<{{T}_{\text{DW}}}$, the Ising critical temperature. These results were verified numerically [312]. This establishes, in an isotropic fully frustrated system, (i) phase coherence is destroyed at ${{T}_{\text{V}}}<{{T}_{I}}$ and (ii) the domain walls roughen at ${{T}_{\text{K}}}<{{T}_{\text{V}}}$ so that the XY-Ising model of (207) describes the long distance behavior near continuous transitions of (206). It is a reasonable model of a fully frustrated isotropic junction array on a square or triangular lattice. This is contrary to the findings in [289] where a difference was found between a square and triangular lattice when f  =  1/2.

There is yet another representation obtained by taking the dual of the XY part which yields the RSOS-Ising model on a square lattice with

$$\begin{eqnarray}&&\frac{H[h,s]}{{{k}_{\text{B}}}T}=\tilde{A}\underset{<\text{I},\text{J}>}{\mathop \sum}\,|{{h}_{I}}-{{h}_{\text{J}}}|+\tilde{C}\underset{<i,j>}{\mathop \sum}\,{{s}_{i}}{{s}_{j}},\end{eqnarray} \tag{ 208 }$$

where the nearest neighbor sites  <i, j >  are on the original square lattice and $<\text{I},\text{J}>$ are the corresponding sites on the dual lattice. ${{h}_{I}}=0,\pm 1$, ${{s}_{i}}=\pm 1$ and $\tilde{A}$ and $\tilde{C}$ are given in terms of A and C of (207), [313]. There is a constraint on the height variables h_I such that a step in h is forbidden to cross a domain wall in s [313]. Another duality transformation on the s_i gives exactly the coarse grained model for CsCl(001) facets [314]. The XY-Ising model of (207) describes the critical behavior of a surprising number of models, including the 4-state Potts model and also contains the much discussed points where Ising and XY degrees of freedom are simultaneously critical at $A-B=\pm {{x}_{\text{c}}}$ where the Ising and XY transition lines cross, as shown in figure 32, center. There has been some discussion that new critical behavior might be found in isotropic systems at a transition in a new universality class between a fully ordered to a disordered phase [296, 309, 315–317]. This was later shown not to happen by Olsson by large scale simulations which were sufficiently accurate to distinguish the two transitions [318, 319]. The phase diagrams of figure 32 summarise my understanding of these fully frustrated planbar rotor models.

Finally, it seems that some issues still remain in the isotropic fully frustrated system despite all the effort expended. In a relatively recent study of the FFXY model and the coupled Ising XY model (207) on the A  =  B line by large scale simulations with up to 10⁶ sites on a square lattice [320], there is agreement with [313] for a wide range of the parameter A and C but disagreement when C  <  −5 when the Ising and XY transition lines seem to merge into a single line at a multicritical point at a critical value C^* with  −5  >  C^*  >  −7 [320]. In the RSOS Ising representation, the two critical lines do not merge into a single line, but just become very close together [313]. This disagreement is now purely philosophical and might be resolved by studying the full XYI model of (207) in the space of all parameters A, B, C. It is unlikely, but not impossible, that a single multicritical point will emerge for $C\ll 0$ with $A-B\to 0$ but such a scenario will require massive computations. The issue of systems with nearby Ising and XY transitions also arises in the competition between surface reconstruction and roughening in (1 1 0) facets of fcc crystals which are described by a two component body centered solid on solid model or a staggered 6-vertex model [321–324]. Unfortunately, this is also insoluble but is amenable to accurate transfer matrix studies.

There are few experiments on these fully frustrated systems because the XY and Ising transitions are too close to resolve in an isotropic junction array. However, they become quite well separated in an anisotropic array [284, 305] but there are some technical issues in fabricating an array with the required anisotropy in the coupling strengths. The helicity modulus in an isotropic array with Hamiltonian

$$\begin{eqnarray}&&\frac{H}{{{k}_{\text{B}}}T}=-\underset{\mathbf{r},\boldsymbol{\mu }}{\mathop \sum}\,\left({{A}_{\mu}}+{{B}_{\mu}}\right)s\left(\mathbf{r}\right)s\left(\mathbf{r}+\boldsymbol{\mu }\right)\text{cos}\left(\theta \left(\mathbf{r}\right)-\theta \left(\mathbf{r}+\boldsymbol{\mu }\right)\right)\\ &&-{{C}_{\mu}}s\left(\mathbf{r}\right)s\left(\mathbf{r}+\boldsymbol{\mu }\right),\end{eqnarray}$$

where $\mu =x,y$ for a square lattice. The measurable helicity moduli are [305] near the low T Ising transition reflects the Ising singularity [305]

$$\begin{eqnarray}&&\frac{{{\gamma}_{\mu}}}{{{k}_{\text{B}}}T}\sim \left(T-{{T}_{I}}\right)\text{ln}|T-{{T}_{I}}|.\end{eqnarray} \tag{ 209 }$$

This very weak singularity was observed in a custom built junction array [325] by measuring the complex sheet impedance $Z(\omega T)=R(\omega,T)+\text{i}\omega L(\omega,T)$ and observing a small peak in the resistivity at ${{T}_{I}}<{{T}_{\text{KT}}}$ well below the BKT temperature obtained from the universal relation $L_{\Box }^{-1}\left({{T}_{\text{KT}}}\right)=\left(8\pi /{{\phi}_{0}}\right){{k}_{\text{B}}}{{T}_{\text{KT}}}$. The custom built junction array and the peak in the resistivity are shown in figure 33.

Zoom In Zoom Out Reset image size

Quantum effects become important at very low temperatures and when the charging and Josephson energies are comparable ${{E}_{\text{c}}}/{{E}_{\text{J}}}=\mathcal{O}(1)$. This requires very small junctions as discussed in a comprehensive recent review [274]. The partition function corresponding to the Hamiltonian of (203) is [275]

$$\begin{eqnarray}&&Z=\mathop{\int}^{}{{\Pi}_{i}}\mathcal{D}{{\phi}_{i}}(\tau )\text{exp}(-S[\phi ]),\\ &&S[\phi ]=\int_{0}^{\beta}\text{d}\tau \left[\frac{{{C}_{0}}}{8{{e}^{2}}}\underset{i}{\mathop \sum}\,\overset{\centerdot}{\mathop{\phi}}\,_{i}^{2}+\underset{<ij>}{\mathop \sum}\,\left(\frac{C}{8{{e}^{2}}}{{\left({{\overset{\centerdot}{\mathop{\phi}}\,}_{i}}-{{\overset{\centerdot}{\mathop{\phi}}\,}_{j}}\right)}^{2}}-{{E}_{\text{J}}}\text{cos}\left({{\phi}_{i}}-{{\phi}_{j}}\right)\right)\right],\end{eqnarray} \tag{ 210 }$$

where $\beta =1/{{k}_{\text{B}}}T$ and $\overset{\centerdot}{\mathop{\phi}}\,=\text{d}\phi /\text{d}\tau$. The first and second terms are charging energies in terms of the voltages at the lattice sites. To include dissipative tunneling, the action includes a term

$$\begin{eqnarray}&&{{S}_{\text{D}}}[\phi ]=\frac{1}{2}\int_{0}^{\beta}\text{d}\tau \text{d}{{\tau}^{\prime}}\underset{<ij>}{\mathop \sum}\,\alpha \left(\tau -{{\tau}^{\prime}}\right)F\left({{\phi}_{ij}}(\tau )-{{\phi}_{ij}}\left({{\tau}^{\prime}}\right)\right),\end{eqnarray} \tag{ 211 }$$

where, for Ohmic dissipation

$$\begin{eqnarray}&&\alpha (\tau )=\frac{\pi}{2{{e}^{2}}{{\beta}^{2}}{{R}_{N}}}\frac{1}{\text{si}{{\text{n}}^{2}}(\pi \tau /\beta )}.\end{eqnarray} \tag{ 212 }$$

The dissipation mechanism might be normal electron tunneling by discrete charge transfer and

$$\begin{eqnarray}&&{{F}_{\text{QP}}}\left[{{\phi}_{ij}}\right]=1-\text{cos}\left(\frac{{{\phi}_{ij}}(\tau )-{{\phi}_{ij}}\left({{\tau}^{\prime}}\right)}{2}\right)\end{eqnarray} \tag{ 213 }$$

which is $4\pi$ periodic. Other mechanisms which yield quadratic forms of $F[\phi ]$ [274] and may force the phase $\phi (\tau )$ to be continuous. If the dissipation is by discrete charge transfer, the path integral implies a sum over winding numbers so that

$$\begin{eqnarray}&&\mathop{\int}^{}\mathcal{D}\phi \to \underset{i}{\mathop \prod}\,\int_{0}^{2\pi}\text{d}{{\phi}_{i0}}\underset{{{m}_{i}}=0,\pm 1,\cdots}{\mathop \sum}\,\int_{{{\phi}_{i0}}}^{{{\phi}_{i0}}+2\pi {{m}_{i}}}\mathcal{D}{{\phi}_{i}}(\tau ),\end{eqnarray} \tag{ 214 }$$

because the charges on the grains are integer multiples of 2e [274].

To obtain the dual Coulomb gas action for a Josephson junction array, one expresses the array partition function in terms of a sum over charge and vortex configurations as $Z=\sum_{q,v}\text{exp}\left[-S(q,v)\right]$ by [274]

$$\begin{eqnarray}&&S[q,v]=\int_{0}^{\beta}\text{d}\tau \underset{ij}{\mathop \sum}\,\left[2{{e}^{2}}{{q}_{i}}(\tau )C_{ij}^{-1}{{q}_{j}}(\tau )+\pi {{E}_{\text{J}}}{{v}_{i}}(\tau ){{G}_{ij}}{{v}_{j}}(\tau )\right. \\ &&\left. +\frac{1}{4\pi {{E}_{\text{J}}}}\overset{\centerdot}{\mathop{q}}\,(\tau ){{G}_{ij}}{{\overset{\centerdot}{\mathop{q}}\,}_{j}}(\tau )+\text{i}{{q}_{i}}{{\Theta}_{ij}}{{\overset{\centerdot}{\mathop{v}}\,}_{j}}(\tau )\right],\\ &&{{\Theta}_{ij}}=\text{ta}{{\text{n}}^{-1}}\left(\frac{{{y}_{i}}-{{y}_{j}}}{{{x}_{i}}-{{x}_{j}}}\right),\\ &&{{G}_{ij}}=-\frac{1}{2}\text{ln}\left({{r}_{ij}}/a\right)+{{E}_{a}}\end{eqnarray} \tag{ 215 }$$

where the charge ${{q}_{i}}(\tau )$ and the vorticity ${{v}_{i}}(\tau )$ are integer valued functions of imaginary time τ. G_ij is the interaction between two vortices separated by r_ij/a lattice spacings with E_a the short distance contribution. Provided the capacitance to ground ${{C}_{0}}\ll C$, the capacitance between neighboring grains, the interaction between charges has the same form as between vortices

$$\begin{eqnarray}&&{{e}^{2}}C_{ij}^{-1}=\frac{{{E}_{\text{C}}}}{\pi}{{G}_{ij}},\end{eqnarray} \tag{ 216 }$$

so that, from (215), the charges and vortices are dual with a self dual point at ${{E}_{\text{J}}}/{{E}_{\text{C}}}=2/{{\pi}^{2}}$ which is exact when C₀  =  0 and in the absence of the spin wave term $\overset{\centerdot}{\mathop{q}}\,G\overset{\centerdot}{\mathop{q}}\,$ in (215). This term is irrelevant for the statics but not for dynamics. These considerations lead to a phase diagram as shown in figure 34, right, which should be compared with the experimental phase diagram of figure 34, center.

Zoom In Zoom Out Reset image size

There are some interesting predictions arising from these considerations related to the experimental observation that, at the border of the superconducting -insulator transition, the 2D system is metallic with a resistivity ${{\rho}^{\ast}}(T=0)\approx {{\rho}_{Q}}=h/{{(2e)}^{2}}$, the quantum of resistance. An oversimplified argument is in [326] based on vortex flow and duality. The Josephson relation for the voltage V across a junction and the current I are

$$\begin{eqnarray}&&V=\frac{\hslash}{2e}\overset{\centerdot}{\mathop{\theta}}\,=\frac{h}{2e}\overset{\centerdot}{\mathop{v}}\,,\\ &&I=2e\overset{\centerdot}{\mathop{q}}\,,\\ &&\Rightarrow \,\,\rho =\frac{V}{I}=\frac{h}{4{{e}^{2}}}\frac{\overset{\centerdot}{\mathop{v}}\,}{\overset{\centerdot}{\mathop{q}}\,}\end{eqnarray} \tag{ 217 }$$

At the self dual point of (215), we have $\langle q\rangle =\langle v\rangle$, which means that for every Cooper pair crossing a square, one vortex crosses perpendicular to the Cooper pair. Thus, at the T  =  0 SI transition the system is metallic with resistivity ${{\rho}^{\ast}}=h/\left(4{{e}^{2}}\right)$. Of course, the self duality is not exact but one might expect that ${{\rho}^{\ast}}=\mathcal{O}\left({{\rho}_{Q}}\right)$, especially as the self dual point is a universal fixed point in the RG sense. Some explicit calculations of ${{\rho}^{\ast}}$ have been performed, mostly as a 1/N expansion [326–328]. In the exactly soluble $N=\infty$ limit, the universal resistivity of a superconducting array ${{\rho}^{\ast}}(0)=8{{R}_{Q}}/\pi$ in zero applied field [327] and, in the fully frustrated case of f  =  1/2, ${{\rho}^{\ast}}(1/2)={{\rho}^{\ast}}(0)/2=4{{R}_{Q}}/\pi$ [328]. Measurements of the T  =  0 universal resistivity [329] are not easy.