# Absence of a Finite-Temperature Melting Transition in the Classical Two-Dimensional One-Component Plasma

\[

## Abstract

Vortices in thin-film superconductors are often modelled as a system of particles interacting via a repulsive logarithmic potential. Arguments are presented to show that the hypothetical (Abrikosov) crystalline state for such particles is unstable at any finite temperature against proliferation of screened disclinations. The correlation length of crystalline order is predicted to grow as $`\sqrt{1/T}`$ as the temperature $`T`$ is reduced to zero, in excellent agreement with our simulations of this two-dimensional system.

PACS numbers: 64.70.-p,74.60.-w,

\]

It has been commonly assumed for many years now that for the physically important case of particles moving in two dimensions interacting with each other via a repulsive logarithmic potential (a situation sometimes called the two-dimensional one-component plasma problem) one would have the usual phases expected on the KTHNY scenario. This scenario describes two-dimensional melting as a defect–mediated phenomenon (Halperin, Nelson  and Young ) and is based on ideas of Kosterlitz and Thouless . It is supposed that the crystalline phase – a triangular lattice – melts at a continuous transition into an hexatic liquid due to the proliferation of dislocations. The hexatic liquid becomes an ordinary liquid at temperatures which permit the creation of disclinations. This scenario is well-established for particles with short-range interactions , but we will argue that particles interacting with a logarithmic potential behave completely differently. For them we can show that the crystalline state is unstable at any temperature against the proliferation of (screened) disclinations and as a consequence the system stays in the liquid state down to arbitrarily low temperatures. The ground state of the system is of course crystalline; the correlation length of short-range crystalline order is predicted to grow as $`\sqrt{1/T}`$ as the temperature $`T`$ approaches zero. Our numerical simulations reported here confirm this behavior.

The one-component plasma problem is of considerable physical significance as it relates to the thermodynamics of vortices –“the particles” – in thin film superconductors. For thin enough films the screening length in the intervortex potential may be greater than the transverse dimensions of the film, which makes the logarithmic potential an accurate approximation for the potential. Most papers on thin film superconductors assume the vortices have a freezing transition at low enough temperatures (for a review see Ref. ), although clear experimental evidence for this is lacking. For a contrary view, however, see  and references cited therein.

For particles interacting with a repulsive potential a device is needed to stop them escaping to infinity. In numerical studies of two-dimensional melting the most commonly used device is periodic boundary conditions. Unfortunately the use of this boundary condition with either short-range interactions or with the logarithmic interaction  produces an apparently first order transition between the crystal and liquid states rather than the KTHNY scenario. (For a review of early work on short-range interactions see ; for some more recent work see Refs.  or ). This is probably a finite size effect: studies on systems with over 60,000 particles indicate that the van der Waals loops associated with the apparent first-order transition shrink in these very large systems , . We ourselves have found that placing the particles on the surface of a sphere is very effective for short-range interactions : no van der Waals loops occur with this topology and the results obtained even with modest numbers of particles are in excellent agreement with expectations based on KTHNY theory. As a consequence the numerical work which we are reporting in this paper has been carried out for the two-dimensional system represented by the surface of a sphere.

The ground state configuration of the particles on the sphere has to contain at least 12 disclinations (5-fold rings) by Euler’s theorem. We have made extensive studies of these ground states and discovered that for larger systems the disclinations are screened by lines of dislocations , . These defects within the crystalline state seem to overcome the problem of the spurious first order transition induced by finite size effects when periodic boundary conditions are employed and so enable one to get results closer to those obtaining in the thermodynamic limit. It is noteworthy that an early simulation of the one-component plasma on the surface of a sphere  did not find a finite temperature phase transition either, in agreement with our results.

The Hamiltonian for particles moving on the surface of the sphere interacting via a logarithmic potential is

$$H=J\underset{i<j}{}ln(|𝐫_i𝐫_j|/R),$$
(1)   
where $`𝐫_i`$ is the position of the $`i`$th particle on the surface of the sphere, $`R`$ is the radius of the sphere and $`J`$ is a measure of the magnitude of the repulsive forces between the particles. The key feature which distinguishes the logarithmic potential from other potentials is that all the stationary states of $`H`$ have zero dipole moment . This was proved by noting that the force on the $`i`$th particle due to all the others must be directed radially for any equilibrium configuration since otherwise the particle would move along the sphere, so

$$\underset{ji}{}\frac{𝐫_i𝐫_j}{|𝐫_i𝐫_j|^2}=f_i𝐫_i.$$
(2)   
By multiplying both sides by $`𝐫_i`$ one can show that $`f_i=(N1)/2R^2`$ where $`N`$ is the total number of particles. By summing Eq. (2) over all $`i`$ and using the fact that $`𝐫_i𝐫_j`$ is antisymmetric in $`i`$ and $`j`$, it follows that the dipole moment, $`_i𝐫_i`$, is zero. No other potential has this feature and it has important consequences.

Our basic contention is that thermally excited screened disclinations will destroy crystalline order for particles interacting with each other via a logarithmic potential. The energy cost of an unscreened disclination is $`O(N)`$  and such a disclination will not be thermally created in a crystalline state. However, disclinations can be “screened” by a cloud of dislocations and it turns out that the energy cost of such screened disclinations can be much smaller, of $`O(ln(N))`$ for non-logarithmic potentials and $`O(1/N)`$ for the logarithmic potential. The phenomenon of screening of a disclination by dislocations is well-known , , and is illustrated in Fig. 1. The figure shows a five-fold coordinated ring – a disclination – screened by five lines of dislocations where the dislocations are spaced a distance $`cl_0`$ apart; $`l_0`$ is the lattice spacing. The strain field of the central disclination can be largely cancelled by that arising from the lines of dislocations. The resulting strain field from the dislocations along a line may be approximated at large distances by that which arises from a positive and negative disclination at each end of the line with a “disclination” charge of size $`q_s=1/c`$ (where $`q_s=+1`$ for a fivefold disclination). If we consider a line of dislocations with $`c=5`$ and five lines as in Fig. 1 then the strain field of the central disclination is exactly screened away. As shown in Ref.  the contribution of the disclinations at the other end of the lines can be made arbitrarily small by allowing the spacing of the dislocations to increase with distance $`r`$ from the disclination as $`c(r)=5+S/l_0g(r/S)`$ with the condition that $`g(0)=0`$; $`S`$ is the size or scale of the screened disclination. Then the residual charge associated with the screened disclination can be made as small as $`O(l_0/S)`$ but, as we shall show, must be as small as $`O((l_0/S)^2)`$ for the special case of the logarithmic potential.

First let us review some features of two-dimensional continuum elasticity theory . Small strains $`u_{ij}(𝐫)`$ are related to the stress field by Hooke’s law, $`\sigma _{ij}=B\delta _{ij}u_{kk}+2\mu (u_{ij}\delta _{ij}u_{kk}/2)`$, where $`B`$ is the bulk modulus and $`\mu `$ is the shear modulus. In the presence of topological defects it is convenient to introduce the Airy stress function, $`\chi `$, defined by $`\sigma _{ij}=ϵ_{ik}ϵ_{jl}_k_l\chi `$. A fivefold disclination is defined by a change in bond angle $`2\pi /6`$ when a path encircles the defect. Dislocations are defined by their Burgers vector density field $`𝐛(𝐫)`$ which for the dislocations in Fig. 1 points perpendicular to the line upon which they lie. The stress field is related to the densities of disclination $`s(𝐫)`$ and dislocations via

$$\frac{1}{Y_2}^4\chi =s(𝐫)ϵ_{ik}_kb_i(𝐫)\stackrel{~}{s}(𝐫),$$
(3)   
where $`Y_2=4B\mu /(B+\mu )`$. For a single disclination at the origin as in Fig. 1, $`s(𝐫)=(2\pi /6)\delta (𝐫)`$. $`\stackrel{~}{s}(𝐫)`$ can be regarded as a total disclination density made up of a “free” disclination density $`s(𝐫)`$ and a “polarization” contribution $`ϵ_{ik}_kb_i`$ from dislocations. The energy of the screened disclination expressed in Fourier space is

$$E=\frac{1}{2}Y_2\frac{d^2q}{(2\pi )^2}\frac{1}{q^4}\stackrel{~}{s}(𝐪)\stackrel{~}{s}(𝐪).$$
(4)   
The expectation  for non-logarithmic interaction potentials is that the screening can at best make the Fourier transform of $`\stackrel{~}{s}`$, $`\stackrel{~}{s}(𝐪)=ql_0f(qS)`$ when the amount of “disclination charge” within a region of radius $`S`$ around the disclination is of $`O(l_0/S)`$. Substituting this form for $`\stackrel{~}{s}(𝐪)`$ into Eq. (4), one finds that in a system of $`N`$ particles the energy of the screened disclination would be of $`O(lnN)`$ – which explains why screened disclinations would be unlikely to modify the KTHNY scenario for non-logarithmic interactions.

The change in the particle density due to the presence of the topological defects is, when Fourier transformed, given by

$$\mathrm{\Delta }\rho (𝐪)=nu_{ii}(𝐪)=\frac{n}{2B}q^2\chi (𝐪),$$
(5)   
where $`n`$ is the number density of the particles. A feature of the logarithmic potential is that for it the bulk modulus $`B`$ is not constant  but diverges at small wavevector: $`B(q)=2\pi Jn^2/q^2`$. The shear modulus is well-behaved: $`\mu =Jn/8`$. Equations (4) and (5) can still be used if the replacements $`Y_24\mu `$ and $`BB(q)`$ are employed. At small wave-vector it follows that for the logarithmic potential

$$\mathrm{\Delta }\rho (𝐪)=\frac{1}{8\pi }\stackrel{~}{s}(𝐪).$$
(6)   
Only for the logarithmic case does a finite small $`q`$ limit exist for the density change $`\mathrm{\Delta }\rho `$ associated with a screened disclination.

We now exploit the fact that all stationary states of $`H`$ have vanishing dipole moment to show that for the case of logarithmic interactions between the particles the screening of the disclination is more efficient than for non-logarithmic interactions. The Fourier transform of the particle density is defined by

$$\rho (𝐪)=\frac{1}{N}\underset{i}{}e^{i𝐪.𝐫_i}.$$
(7)   
(Formally, as our system is the surface of a sphere rather than a plane we should use spherical harmonics rather than plane waves, as was done in Ref., but the distinction is unimportant for our argument). Consider now a state which differs from the groundstate by the presence of a screened disclination of size $`S`$. The density difference $`\mathrm{\Delta }\rho (𝐪)`$ of the two states must differ as $`q0`$ by terms of $`O(q^2)`$ ; (if one of the states had had a dipole moment then one can see from expanding the exponential in Eq. (7) that $`\mathrm{\Delta }\rho `$ would have been of $`O(q)`$). Eq. (6) implies that $`\stackrel{~}{s}(𝐪)=q^2l_0^2f(qS)`$. By Fourier transforming $`\stackrel{~}{s}(𝐪)`$ one can then show that the “disclination charge” within a distance S of the center of the disclination is of $`O((l_0/S)^2)`$.

Substituting this form for $`\stackrel{~}{s}`$ into Eq. (4) one finds that the energy of the screened disclination is of order $`J(l_0/S)^2`$. By increasing the scale S it can be made arbitrarily small. At a temperature $`T`$ a region of linear extent $`\xi `$, where $`J(l_0/\xi )^2=T`$, will be unlikely to contain a disclination and so $`\xi `$ will be a measure of the short-range crystalline order present in the system at temperature $`T`$. $`\xi `$ diverges as $`\sqrt{1/T}`$ as $`T0`$. This means that by investigating numerically the structure factor one can find from the widths of its peaks the correlation length $`\xi `$ and its temperature dependence will tell us whether the arguments above are valid.

We studied using molecular dynamics, specifically a velocity Verlet algorithm , a system of $`N`$ particles confined to move on the surface of a sphere and interacting with the logarithmic potential. Reduced units were used, i.e. $`m=k_\mathrm{B}=R=J=1`$, where $`m`$ is the mass of the particle and $`k_\mathrm{B}`$ is the Boltzmann constant. The acceleration $`𝐚_i`$ of the $`i`$th particle equals $`𝐟_i/m`$, where $`𝐟_i`$ is the force produced by the other particles on the $`i`$th particle. After a small time interval $`\delta t`$ the position of the particle will be $`𝐱_i=𝐫_i(t)+𝐯_i(t)\delta t+\frac{1}{2}𝐚_i(t)\delta t^2`$, where $`𝐯_i(t)`$ is the velocity of the particle. In general $`𝐱_i`$ will not lie on the surface of the sphere. The $`i`$th particle is brought back to the surface by acting on it with a fictitious force $`2\lambda _i𝐫_i(t)`$ where

$$\lambda _i=\frac{𝐫_i(t)𝐱_i\sqrt{\left[𝐫_i\left(t\right)𝐱_i\right]^2R^2\left[\left|𝐱_i\right|^2R^2\right]}}{R^2\delta t^2}.$$
(8)   
Then, the velocity Verlet algorithm updates particle positions using the equation:

$$𝐫_i(t+\delta t)=𝐱_i\lambda _i𝐫_i(t)\delta t^2.$$
(9)   
We chose $`\delta t=0.005(mR^2/J)^{1/2}`$. The velocities of the particles were chosen from a Boltzmann distribution appropriate to the temperature $`T`$ and were re-selected at equally spaced time intervals . The system was equilibrated at high temperatures, then the temperature was slowly reduced. For each temperature we determined the structure factor which is related to the Fourier transform of the pair correlation function $`h(r)`$ by :

$$S(q)=1+2\pi \rho R^2_0^\pi h(R\theta )\mathrm{sin}\theta J_0(qR\theta )𝑑\theta ,$$
(10)   
where $`J_0`$ is the Bessel function of zeroth order. This adaption of the conventional relation to particles moving on the surface of the sphere is valid provided $`q`$ is of $`O(1)`$ and not of $`O(1/R)`$. Peaks in the structure factor grow as the temperature is reduced at wavevectors $`q`$ corresponding to the reciprocal lattice vectors $`|𝐆|`$ of the triangular lattice expected for the groundstate. The correlation length, $`\xi `$, is the inverse of the width of the first peak of the structure factor. To determine it, we fitted the first peak to a Lorentzian curve.

We studied systems of 1442 and 2252 particles. The simulation time for each temperature was 100,000$`\delta t`$. In Fig. 2 $`\xi `$ is plotted against $`\sqrt{1/T}`$. The vertical line is drawn where other authors found a first-order melting transition using periodic boundary conditions . The predicted behavior $`\xi \sqrt{1/T}`$ is clearly seen in Fig. 2. When the temperature was reduced to $`T=0.01`$ the correlation length for the system of 1442 particles reached the system size and stopped growing as the temperature was reduced further. However, the simulation with 2252 particles indicates that this levelling off is just a finite size effect. True equilibration in numerical studies of two-dimensional melting phenomena is always problematic  and may be the cause of the scatter in Fig. 2.

The apparent absence of a finite temperature phase transition for particles interacting with a logarithmic potential cannot be attributed to the fact that we have done the simulation for the two-dimensional geometry represented by the surface of a sphere. We can demonstrate this by simulating particles interacting with a $`1/r^{12}`$ potential also moving on the surface of a sphere. A finite temperature melting transition of KTHNY character is seen. We had already found indications of such a melting transition  but the numbers of particles used in that reference were rather small. Using the Verlet algorithm described above we are able to simulate much larger systems eg. 5882 particles. The calculations for the short-range potentials run faster than with the logarithmic potential as it is possible to use look-up tables of nearest neighbors. In the KTHNY picture, the correlation length has the following density dependence along an isotherm for a $`1/r^{12}`$ potential :

$$\xi (\rho )\mathrm{exp}\left(\frac{b}{\left((\rho _c/\rho )^61\right)^\nu }\right),$$
(11)   
where $`\rho `$ is the density and $`\nu =0.36963\mathrm{}`$. In Fig. 3, ln$`\xi `$ is plotted versus $`(1/\rho ^61)^{0.36963}`$. We have assumed that $`\rho _c=1`$, a value obtained by other authors  working at the temperature which we used, $`T=1`$. The slope of the straight line is $`1`$ which corresponds well with KTHNY predictions. Note again that finite size effects cut off the growth of $`\xi `$ when it is of $`O(R)`$. Thus simulations on the sphere do produce for a non-logarithmic potential the expected crystalline phase.

In summary, we have shown that thermal excitation of screened disclinations removes at non-zero temperature the crystalline phase of the vortex system. Numerical simulations have confirmed our prediction that as the temperature is lowered the correlation length of short-range crystalline order should grow as $`\xi \sqrt{1/T}`$.

APG would like to acknowledge a grant and financial support from the Dirección General de Investigación Científica y Técnica, project number PB 96/1118 and EPSRC under grant GR/K53208. We have had useful discussions with H. Bokil.