Phase Diagram of the Two Dimensional Lattice Coulomb Gas

ar Xiv:c on d -m a t/9609031v 13 S e p1996Phase Diagram of the Two Dimensional Lattice Coulomb GasPramod Gupta and S.TeitelDepartment of Physics and Astronomy,University of Rochester,Rochester,New York 14627(February 1,2008)We use Monte Carlo simulations to map out the phase diagram of the two dimensional Coulomb gas on a square lattice,as a function of density ρand temperature T.We find that the Kosterlitz-Thouless transition remains up to higher charge densities than has been suggested by recent theo-retical estimates.The nature of phase transitions in the two dimensionalneutral Coulomb gas (CG)has remained a topic of con-siderable interest.The CG can be related via dualitytransformation to the 2D XY model,and thus to super-fluid and superconducting films.1The pioneering workof Kosterlitz and Thouless 2(KT)showed that,at low charge density,there is a second order transition from aninsulating to a conducting phase,due to the unbinding ofneutral charge pairs.As the charge density is increased,several authors 3–6have predicted that this KT transitionshould become first order.Recently,Levin et al.5,usinga modified Debye-H¨u ckel approach,have estimated that,for a continuum CG,the KT transition ends in a tricriti-cal point at the surprisingly low density of ρc ≃0.004/a 2,where a is the hard core diameter of the charges.To investigate this issue,we report here on new MonteCarlo simulations of the 2D CG on a square lattice,com-puting the phase diagram as a function of density andtemperature.We make sensitive tests of the nature ofthe transition,and conclude that it remains second order up to densities much higher than estimated by Levin etal.Our work follows that of Lee and Teitel,7with a fewmodifications.We take for the Hamiltonian,H =1L +V 2πP y2ln(8e 2γ),where γ≃0.5772is Euler’s constant.2,8n i =0,±1,±2...are the integer charges,and neutralityi n i =0is imposed.The third termin H tends to suppress charges with |n i |>1,and isneeded to stabilize the system in the very dense limit.Assuming therefore that all charges satisfy |n i |≤1,so that |n i |=n 2i ,the second term is just −uρL 2withρ=L −2i |n i |the charge density.Thus u is the chem-ical potential.The last two terms in H are effectiveboundary terms,which arise in the duality mapping tothe CG fromthe 2D XY model with periodic boundary conditions.9–11Here P ≡i n i r i is the net dipole mo-ment of the charges,and V (φ)is the Villain function,12e −V(φ)/T ≡ ∞m=−∞e −1L δµν−2πL V ′ 2πP ν0.20.40.60.81.01.2T r i c r i t i c a l P o i n t I s i n g (M e l t i n g )K T c h a r g e d e n s i t y ρI n s u l a t i n g S o l i d C o n d u c t i n g S o l i d I n s u l a t i n g G a s C o n d u c t i n g L i q u i d P h a s e S e p a r a t e d R e g i o n FIG.2.Phase diagram of the lattice CG in the ρ−T plane.Open symbols are from simulations with system size L =16;solid symbols are from L =32.0.0000.0020.0040.0060.0080.0100.0120.0140.30.40.50.60.70.80.91.0P r o b a b i l i t y P (ρ)T = 0.142T = 0.139T = 0.137L = 32 u = 0.4 > u 0L i q u i d C o e x i s t e n c e S o l i d FIG.3.Histograms of charge density ρfor u =0.4just above u 0=π/8.At the coexistence boundary the histogram is bimodalobtained by simulating with u just above and just below u 0,measuring the average density ρas T increases.To determine the boundary closer to the tricritical point,we simulate with fixed u >u 0,increasing T ,and measuring the histogram of the values of ρfound at each value of T .When one is in either the solid or the liquid phase,this histogram has a single peak.However,when one crosses the first order line,this histogram develops double peaks.15The locations of the two peaks determine the densities of the two coexisting phases.In Fig.3we show an example of such histograms.Varying u >u 0then maps out the rest of the coexistence boundary.In Fig.2we show the coexistence boundary found in this way from simulations with L =16and L =32.As is seen,and as is expected,our results near the tricritical point are limited by finite size effects.However it is clear that the KT line at small u joins the first order line at a density ρ≃0.1,much larger than the estimate of Levin et al.Next we verify that what we have called the “KT”line at small u actually does remain a second order KT transition all the way up to the first order line separating the insulating gas and insulating solid.To locate the transition temperature T ∗for u =0.39,just below u 0,we compute ǫ−1(T )for various L ,using Eq.(2).Our results are shown in Fig.4.The intersection of these curves with the line 4T determines the upper bound T ∗≃0.125.No hysteresis,or other suggestion of a first order transition 0.20.30.40.50.60.70.80.90.1180.1220.1260.130L =64L =32L =16L =8in v e r s e d i e l e c t r i c c o n s t a n t u = 0.39 < u 0 = 0.392699ε−14T FIG.4.Inverse dielectric function ǫ−1(T )for various lattice sizes L ,for u =0.39just below u 0=π/8.The intersection with the line 4T gives the KT upper bound on the transition temperature T ∗.was observed in ǫ−1.As a more precise test that the transition is indeed KT–like,we use the finite size scaling procedure of We-ber and Minnhagen 16,17.Precisely at the KT transition temperature,the finite size dependence of ǫ−1is given by,ǫ−1(L )=ǫ−1(∞) 1+10102030400.1200.1210.1220.1230.1240.1250.126χ2offit toE q.(3)u=0.39 < u 0(a)3.54.04.55.00.1200.1210.1220.1230.1240.1250.126L=8,12,16,24,32,48,64L=12,16,24,32,48,64L=16,24,32,48,64fit te dv alueofε-1(∞)/T Temperature Tuniversal KT criterion(b)T *FIG.5.Fitting of ǫ−1(T,L )to Eq.(3).(a )shows the χ2of the fit vs.T ,(b )shows the value of the fitted ǫ−1(∞)/T vs.T .0.0000.0050.0100.0150.0200.0250.0300.000.050.100.150.200.25P r o b a b i l i t y P (ρ)L = 64 u = 0.39 < u 0T =0.110T =0.115T =0.120T =0.125T =0.127T =0.130FIG.6.Histograms of charge density ρfor u =0.39just below u 0=π/8.Histograms remain single peaked as one passes though the transition at T ∗≃0.1235of excited isolated dipoles and quadrupoles.For the lower branch we get ρ−(T )=1L 2(2N ′d e −βE ′d +4N ′q e −βE ′q )where E ′d and E ′q are the energies for removing an isolated dipole and quadrupole,and N ′d =2L 2and N ′q =L 2are the number of ways this may be done.Clearly these expressions will change with the geometry of the discretizing lattice,or if a continuum is used.Nevertheless,at the very low den-sities ρ∼0.004and high temperatures T ≃0.25where Levin et al.estimate a tricritical point,we would be very surprised if the lattice is qualitatively different from the ing a discrete lattice also has the effect that it tends to stabilize the charge solid phase above u 0to high tem-peratures.Indeed,our charge solid only melts after it has already become conducting via a KT transition arising from the excitation and diffusion of vacancies through-out the solid.The present model does not possess any first order transition from an insulating to a conduct-ing phase.In contrast,recent simulations 18,19of the CG in a flat continuum with periodic boundary conditions find that the charge solid phase is melted at any finite temperature.Here,a first order line separates the in-sulating gas from a dense conducting liquid,ending at a critical point at relatively low temperature and high density:(T c ,ρc )=(0.056,0.21)according to Ref.18,and T c =0.032according to Ref.19.Similar results were found earlier for the continuum CG on the surface of a sphere:20(T c ,ρc )=(0.087,0.11).In these models,the KT line ends either at,or near this critical point.Although the geometry of the the CG system clearly affects the location of the end of the KT transition line,it is interesting to compare our results for the square lat-tice with the predictions of the continuum self-consistent screening theory of Minnhagen and Wallin.3To do so,it is necessary to note 2,8that if the interaction G ( r )on the lattice is chosen so as to asymptotically match the continuum −ln r as r →∞,then the lattice CG with u =0acts like a continuum model with a chemical potential µ0=−1。