Numerical simulations on two-dimensional vortex lattices

PHYSICA

Physica C 183 ( 1991 ) 212-220 North-Holland

Numerical simulations on two-dimensional vortex lattices D. Reefman and H.B. Brom Kamerlingh Onnes Laboratory, Leiden University. PO Box 9506, 2300 RA Leiden, Netherlands

Received 12 August 1991 Revised manuscript received 8 October 1991 We have studied the dynamics of a two-dimensional vortex lattice with the use of a molecular dynamics method for systems of up to 2000 vortices. The interaction potential is taken to be of the Pearl form and temperature is taken into account via a white Gaussian noise source. We find that the low temperature state is predominantly a triangular array of two-dimensional vortices. Long range order is absent. Defects are present at very low temperatures already, and possibly cause at higher temperatures Kosterlitz-Thouless (dislocation mediated) melting of the lattice. With a London penetration depth of 1250 ,~, we calculate the field distribution at zero temperature and the T-dependent magnetic field time-correlation function. The results compare well with NMR experiments on TI in T12Ba2CaCu2Oa in the temperature range below the melting temperature Tm.

1. Introduction In the past two decades much attention has been p a i d to the structure a n d d y n a m i c properties o f vortex lattices in type II superconductors. It has been well established since then that in general the vortices form a triangular ( A b r i k o s o v ) lattice. D y n a m i c properties can be dealt with by introducing an elasticity tensor ¢. Irreversible properties can be described well by e.g. thermally activated flux creep. Also fluxline melting has been the subject o f intensive research in the recent past, a n d m a n y theoretical a n d numerical studies on this subject have a p p e a r e d [1-8]. W i t h the discovery o f the highly anisotropic highTc superconductors TI2Ba2CaCu2Os and Bi2Sr2CaCu208 interest in the intermediate state (i.e. Bc~ < B < B c 2 ) properties o f purely two-dimensional systems has revived. In 1964 Pearl [9,10] d e r i v e d the current distribution o f a vortex in a superconducting slab with a thickness d smaller than the p e n e t r a t i o n depth 2. The intervortex interaction due to this current behaves as r - ~ for large distances. This Coulomb-like interaction gives rise to a n u m b e r o f effects very different from the 3D case. The compressibility o f the lattice becomes very large, a n d the energy difference be-

tween various types o f flux lattices (square or triangular) becomes negligible as long as the vortex density is kept constant. Another consequence o f twod i m e n s i o n a l i t y is lattice instability, in the sense that perfect long range o r d e r is absent, and, for elevated temperatures, the possibility o f a K o s t e r l i t z - T h o u less transition ( " d i s l o c a t i o n m e d i a t e d m e l t i n g " ) [ 11,12 ]. It is, however, not even clear that a vortex lattice should exist; the energy per vortex associated with the Coulomb-like interaction between vortices diverges for an infinite lattice. A good starting p o i n t in the study o f the flux phase o f the anisotropic high-To materials m a y be to study the properties o f a single superconducting stab. The properties o f the bulk superconductor are then a superposition o f the single slab properties, if we m a k e the a p p r o x i m a t i o n o f representing the bulk superconductor by a set o f uncoupled superconducting planes. As at quite low t e m p e r a t u r e s already the different layers in T12Ba2CaCu2Og are d e c o u p l e d [ 13 ], this m a y not be a too b a d a p p r o x i m a t i o n . It is then the m e n t i o n e d dislocation m e d i a t e d melting o f the vortex lattice that m a y be responsible for e.g. the line narrowing observed in N M R experiments on the layered superconductors. In the N M R experiments on Tl2Ba2CaCu2Os ( T o ~ 101 K ) b y Jol et al. [14] a n d M o o n e n et al. [ 15 ] a constant N M R linewidth down

0921-4534/91/$03.50 © 1991 Elsevier Science Publishers B.V. All fights reserved.

D. Reefman, H.B. Brom / Numericalsimulations on 2D vortexlattices to about 60 K is observed. Below this temperature the line width starts to increase linearly and amounts to more than 300 kHz at 4 K. This phenomenon cannot be explained by the temperature dependence of the penetration depth 2, as is the case in the more three-dimensional YBCO [ 16 ], which stimulated the interest in the magnetic properties of a Single superconducting layer. In the many papers that have appeared on the twodimensional vortex crystals, only for very limited temperature ranges analytical results on the dynamic properties have been derived. We therefore have developed a molecular dynamics method to simulate a large number of vortices in a 2D slab. Though Monte Carlo based methods are more efficient in finding various temperature dependent properties, our method enables us to watch the system evolving in real time, which in turn offers the possibility to calculate various time-correlation functions needed for comparison with e.g. the N M R experiments. The paper is organized as follows. In section 2 we describe the model developed by Pearl, and the numerical methods we have used. In section 3 we discuss the results of the calculations, and compare them with experiment and with analytical expressions obtained by Fisher [ 17 ]. In section 4 we present our conclusions.

2. Model and methods

213

is a Bessel function of the second kind of order n. From the limiting expressions we see that, although the screening length is given by A, the magnitude of the current density falls off only slowly with distance, in contrast with the 3D case. This is due to the fact that in a bulk superconductor the screening current can get very large due to the infinite extension in the z-direction, whereas in the 2D case the total screening current must come from the large radial extent of the current density. The interaction between vortices in a two-dimensional plane, therefore drastically differs from the three-dimensional interaction. For two dimensions the potential energy due to the interaction between two vortices at sites r~ and rj is given by: ~ 2o

r,j) = 2uo- T) U(ro) ~

[ ro \ ¢'~ ln~A-~) nltoA (T)

ro<
U(rij)~ *~ ( 1 ~

ro>>A(T) .

Ito X \riy /

(2)

Here to= Ir,-rjl is the distance between the vortices. Because we use the London approach of nonoverlapping vortex cores, i.e. the density of vortices is low, the potential energy Uv for a vortex at site r/ can be written as the sum of the contributions of interactions with all vortices: N

Uv(r~)=½

2. I. The model As is shown by Pearl [9,10 ], the current distribution due to a single vortex in a superconducting slab of thickness d with d<<2, ;t is the bulk London penetration depth, is given by: q~o

r

J(r)= 2 - ~ [ H I ( A - ~ ) - Y I ( A - - - ~ ) ~o [ 1 ~ j(r).~ -n---~~-~)lo

Oo (l)r lo

j(r)~ n~A(r)

r

2]lo ,

~

j=l,j#i

U(r o)

(3)

where N is the total number of vortices and U(r)defined in eq. (2). From this equation we see that the total potential at position r diverges, as the energy is approximately linear in the size of the system under consideration. This infinite energy however is cancelled up to a constant by the field energy, Uf, given by: Uf(r) = - ½Hogo J d2r ' [r' × j ( r - r ' ) ]z,

r>>A(T), r<
A

(l)

Here A = 222/d is an effective penetration depth for the slab, H . is a Struve function of order n, and Y.

l/_/ R2 -4' 2 x/~-R'~ Uf(r) ~ - ~ o/~o~--~r L~ -~-~-r )

R>>r>>A(T) ,

(4)

which is also linear in the radius R of the system. In

214

D. Reefman. ll.B. Brom / Numericalsimulationson 2D vortexlattices

eq. (4) A is the area of the superconducting slab; the field is taken along the z-axis. E is a complete elliptic integral of the second kind. If we replace the sum in eq. (3) by an integral over the total area A and take the density n of vortices as a constant, then the expression for the total magnetic energy of a superconducting slab in a field becomes:

4A2(T) Ho r

2

r

bz(r),~ 2nA2(T )

n _t- d2r[ Uv(r) + Uf(r) ] ~ -R2Ho B

q~o A ( T)

A

b~(r)~, 2n R> > >A(T) ,

r3

r>>A(T).

(8)

(5)

which compensates the loss in condensation energy in the vortex cores. In the calculation the field energy has to be evaluated numerically for each position, which has to be done only once for each field and system size R. The size of the system is set of course by the number of vortices and the applied field. In the present study we have taken the field to be 5 T, which, in a triangular lattice, corresponds to a lattice constant of 220 A, and have taken sets of 250, 500, 1000 and 2000 vortices. With the parameters 2(0) and d for TlzBa2CaCu2Os, set 1250 A and 10 A respectively, A(0) is chosen as 30 ~tm. Tc is taken 100 K. The temperature dependence of 2 ( T ) is taken from ref. [ 18 ]. Also we have not approximated the potential of eq. (2) by any more convenient form.

2.2. NMR linewidth The lineshape measured in N M R experiments is a reflection of the field distribution in the sample, because the resonance frequency to of a nucleus is proportional to the local field: to(r)=y/roB(r) (7 is the gyromagnetic ratio). Therefore knowledge of the "field density" function X ( B ) is sufficient for calculating the static lineshape. The field density X ( B ) is defined as:

The total field B is the superposition of the fields of all vortices. When the local field, however, becomes time dependent (motion of the vortices) the line shape will be changed due to motional narrowing. To correct for this effect, one can approximate the static line profile by a Gaussian shaped curve (which turns out to be a good approximation ). The line profile F(to) for a time dependent field is then given by [ 19 ]

F(to) = J dtG(t)e -i'°' ,

(9)

where G(t) is given by l

G(t)=ei°~°texp[-to2f (t-z)go~(z)dz],

(10)

0

where top is the linewidth for a time independent field, and 090 is the position of the absorption. The field correlation function go(t) occuring in eq. (10) is defined as

go(t) = ( t o ( 0 ) t o ( t ) ) / ( o 9 ( 0 ) 2) .

( 11 )

Here the brackets denote an average over all space, which is calculated with a Monte Carlo method. For most temperatures except the very low ones go(t) could be approximated well by an exponential:

g,o(t)~,e -t/~, dr~(B'-B(r)),

Y(B')= ~

(6)

In this case eq. (10) goes over in

A

where A is the total area of the specimen. The field b(r) due to a single vortex is found from the London equation: b ( r ) -- - ~ A ( T ) V x j ( r )

resulting in

(7)

G( t ) = ei~t( -to2z2 [ e x p ( - t / % ) - 1 + t/% ] ) . (12)

From eq. (12) it is clear that for small % the lineshape changes from broad, Gaussian-like, to narrowed Lorentzian-like.

D. Reefman, H.B. Brom /Numerical simulations on 2D vortex lattices 2. 3. Molecular dynamics algorithm The equations of motion for the set of N vortices are given by r/v, ( t ) = - g r a d [ Uv + Uf],, +R~ (t) i=1, ..., N ,

(13)

where r/ is a friction coefficient, given within the framework of the Bardeen-Stephen model by [ 20 ] h B ?/~ 4~2 Bc2 "

(14)

h is Planck's constant. R ( t ) is a random force corresponding to a temperature T as defined by the fluctuation dissipation theorem [ 21 ]:

( R~( t)gB( t' ) ) =2rl~c,p~( t - t ' ) k a T ,

(15)

where ka is Boltzman's constant. In the frequency range we are interested in (i.e. frequencies below 100 MHz) we can assume the vortices to be massless. The set of eqs. ( 13 ) is solved by a molecular dynamics (MD) method. In this simulation we used a modified "leapfrog" algorithm:

x(t)=x(t-At)+

l[v(t-At)+v(t-2At)]At;

(16)

x and v are both 2N dimensional vectors representing the positions and velocities respectively of the N vortices. The algorithm is more stable for relatively large timesteps than the original method [ 22 ], which does not involve the velocity at t - 2 A t . This is because oscillations around minima tend to be reduced due to the "memory effect" of the second term in eq. (16). In the discrete integration of eq. (13) one must check whether thermal equilibrium is reached at the correct temperature T [ 6 ]. In thermal equilibrium the following relation holds: ( ( A x ) 2 ) = 2kB____T, r/

(17)

where ( ( A x ) 2) is the mean squared displacement per vortex in a time r >> At. In practice this means that, after application of the Langevin force Ri, some equilibration steps may be necessary to reach ther-

215

mal equilibration. When At is very smail, special care is needed to fulfill eq. ( 17 ). In all simulations the vortices were placed at random positions at t = 0 and T=0.9Tc. We then allowed the vortices to relax, with very small time steps, until the total energy of the system did not change steeply anymore. After this step, we gradually decreased the temperature and increased the time steps, until a new equilibrium was reached which we will call the ground state of the system. This procedure of starting at high temperature and gradually lowering the temperature ensures that the system explores configuration space completely and will not easily settle down in a local minimum. With this ground state we then ran simulations at different temperatures, adjusting the time steps such that the displacement in a single time step was less then 0.5% of the average vortex-vortex distance. For each temperature run, consisting of at least 100,000 time steps for good statistics, we calculated the following correlation functions:

g(t) = ( ( r ( t ) - r ( O ) )2) , )p(r)) , g(r) = ( p ( O (p)2

(18)

where p denotes the density of vortices and the brackets indicate an average over all space. At tern, peratures close to temperatures where the lattice became unstable runs were repeated with 500,000 time steps.

2.4. Computational details All simulations were run on either a SUN sparc workstation or an IBM RS6000 workstation.

3. Results and discussion

3.1. Low temperature We have calculated the ground state configurations for the circular systems of 250, 500, 1000 and 2000 vortices respectively, all in a field of 5 T. These vortex numbers Nv correspond to sample diameters of 0.28, 0.39, 0.56 and 0.78 ~tm respectively. When looking at fig. 1, where the results for N v = 2 5 0 and

2 !6

D. Reefman, H.B. Brom /Numerical simulations on 2D vortex lattices 0.2

(a)

Q

g

0.1

>,

-0.1

-0.2 -0.2

0.6

I -0.1

b

'

( )

~

I 0 × (pro)

I 0.1

'

'

.... .v.'.'.':

0.2

....

0"30

( 2 ) ) however, the precise nature o f the vortex lattice is o f little influence to the total energy, as long as the density is kept constant. In fig. 2 we show the energy difference for a square a n d triangular lattice for both the usual 3D type v o r t e x - v o r t e x interaction [24] and for the 2D type interaction, both for an applied field o f 5T. As is well-known, for the 3D case the energy difference between these two lattices is quite large; o f the o r d e r o f 2 percent. In fig. 2 however, no significant energy difference between the two vortex arrays can be seen. Also, the long range interactions become clear in fig. 2. F o r the 3D type interaction the energy per vortex stays constant if more vortices are a d d e d at distances large c o m p a r e d to 2. On the contrary, for the 2D type interaction the energy per vortex increases almost linearly when m o r e vortices are added. Still, the hexagonal structure is the more stable one from out m o d e l calculations• The very small energy difference between the different types o f ordering, however, makes the lattice somewhat m o r e a p t to disorder than is the case for a 3D VL. This corresponds well with the models p r o p o s e d by Kosterlitz a n d Thouless [ 1 1 ] for the melting o f a two-dimensional lattice. In their m o d e l lattice melting takes place in a few steps. Already at low t e m p e r a t u r e lattice dislocations exist. Because the energy o f a single dislocation increases logarithmically with the size o f the dislocation, at low t e m p e r a t u r e s only dislocation pairs occur. At elevated t e m p e r a t u r e s these pairs can

- 0 " 36 m

-0.

i

-0.6

-0.3

l

0

0.3

~•6

x (~)

Fig. 1. Ground state configurations for the sets of 250 (a) and 2000 (b) vortices.

s g .o v

c uJ

2000 are shown, the most p r o m i n e n t feature is the almost perfect triangular arrangement o f the vortices. This is in agreement with the h y d r o d y n a m i c calculations o f Fetter a n d H o h e n b e r g [23] on the basis o f which they concluded that a triangular vortex lattice ( V L ) in a superconducting slab is stable. Due to the extremely long range interaction (see eq.

-200

0

5~00

10~00 Number

ZOO

15JO0

600

I000

20~00

25 O0

of v o r Z Icoe

Fig. 2. Energy per vortex in a 2D slab for a square and triangular lattice (curves coincide on this scale) and (insert) a 3D system (upper curve: square lanice, lower curve: triangular lattice). 2 is taken to be 250 A.

D. R e e f m a n , H.B. B r o m / N u m e r i c a l simulations on 2D vortex lattices

dissociate into single dislocations, which can drift through the lattice, thereby causing melting of the lattice. Indeed in fig. 1 (b) only dislocation pairs can be observed. The absence of long range order is also clearly demonstrated in fig. 3 where we have depicted the vortex-vortex radial distribution function g(r) for the set of 2000 vortices. The distance r is measured in units of ao, the lattice parameter for an ideal triangular VL. The peak structure due to the neighbours disappears gradually for larger distances. Because the longe range interaction tries to keep the density constant, the effect of dislocations on g(r) will be significant only for larger distances, though there is some smearing of the r=~l ao peak visible already. Interestingly, the structure in g ( r ) disappears faster for the sets of 1000 and 2000 vortices than for g(r) of the smaller sets. Probably this is indicative for size effects. As the vortex-vortex correlation functions g(r) of the sets of 1000 and 2000 vortices are identical, for these systems the finite size seems not to be important. In fig. 4 we depict the positions (in units ofao also) a single vortex typically visits around its equilibrium position at low temperature (1 K). The maximum displacement from the equilibrium position is roughly 1% of the intervortex distance ao.

3.2. Higher temperatures When the temperature is raised, but is still fairly low, i.e. there is no melting of the lattice, the results

7

i

i

I

4

8

12

6

5 q t. v

3

1 0

i 0

r

16

(a.)

Fig. 3. Vortex-vortex radial distribution function g ( r ) for the set of 2000 vortices.

217

0.010

'

0.005 •..: .'.: ~.:,".,:~.~.:... %.,,.?. "~

• :-,'r.-~

.',.".,

..

. . . . . ..r% ::/~
O v >.

.

t., : •. -a

t .,s1.~.~.j ~..-. ..... -: > .

• .

i

-0.005

-0.010 -0.015

i -0.0[0

i -0.005 ×

i 0

0.005

(ao)

Fig. 4. Positions typically visited by a vortex at low temperature in a sequence of about 500 time steps•

for the simulations for small and large numbers of vortices differ. For the simulations with a small number of vortices, we observe that for non-zero temperatures sometimes a cluster of three or four vortices makes a jump. Probably this is a transition from local triangular order to square array order, or discrete motion of a lattice defect (dislocation pair). This effect is observed for the simulations on the set of 250 vortices. For simulations with 150 vortices, no motion up to temperatures close to Tc is observed. This is in contrast with simulations on larger numbers of vortices. For large scale simulations it is not possible anymore to observe distinct movements of the vortices; instead, there is some diffusion, i.e. the mean square displacement of a vortex is approximately linear in time. This behaviour is observed for the simulations on the set of 2000 vortices. Simulations on the intermediate set of 500 vortices show a behavior that is somewhat in between the discrete "hopping" and the continuous diffusion. Therefore, we believe that the "hopping" phenomenon is merely an artefact of the small size of the system. In fig. 5 we have plotted the diffusion coefficients for the system of 2000 vortices for various temperatures. The diffusion constant D is defined by [ 21 ]: g(t)=

4Dt,

(19)

218

D. Reefman, H.B. Brom / Numerical simulations on 2D vortex lattices

i

vortex systems, the latter explanation may be the most plausible. Limiting formulas for various magnetic field ranges for the temperature Tm where dislocation mediated melting should occur have been derived by Nelson and Halperin [25], and Fisher [ 17]. If we take the regime for Be1 <
i

0.030

~. o. o2o

0.010

1 ¢~ Tm = 16x2x/~ ksltoA(Tm)"

0 ZS

50

(20)

75

r (K)

Fig. 5. Diffusion constant D vs. temperature for the set of 2000 vortices.

where g(t) is defined in eq. (18). Within the scatter, the diffusion coefficients are increasing almost linearly with temperature up to roughly 50 K. Above this temperature, the diffusion coefficient D starts to increase rapidly. We attribute this steep rise to the melting of the vortex lattice. It is worth noting that a similar phenomenon occurs for the smaller set of 250 vortices. At a temperature of 60 K the number of "hops" per unit time increases suddenly. At still higher temperatures the structure of the hops disappears, and diffusive behavior appears. The linear increase in D observed at low temperatures is in contrast with the exponential increase found by Jensen and Brass [ 5 ]. In their simulations, the interaction between the vortices was of the 3D type, i.e. exponentially decaying for distances larger than 2, and randomly distributed pinning centers were included. Also, at very low temperatures, diffusion was absent. This difference in low temperature behavior can have several causes. Firstly, Jensen and Brass used periodic boundary conditions in their simulations. In our model simulation simple use of periodic boundary conditions was impossible, due to the longe range of the interaction. Secondly, due to the use of a short ranged force between the vortices, melting is in their model a very local process, in contrast with our model where melting is a global process. As we do not see a large difference between the diffusion coefficients for various sizes of

If we take A in this expression equal to 30 pm, the same value we have used in our simulations, a melting temperature of approximately 30 K is obtained. This value is substantially lower than the numerical result. Although our simulations do not show a significant size effect in the melting temperature, except for the set of 250 vortices, this discrepancy still may be due to the limited number of particles in the simulations. For the largest set we have studied, the dimensions of the substrate are still much less than A, and the effect of the very long range interactions is ignored. However, due to the positional disorder, the effective potential due to the vortices at a large distance is, though large (infinite), structureless. In our view it will thus probably not be of much influence to the melting temperature.

3.3. Field distribution and NMR line shapes In fig. 6 we show a cross section of the field profile of the slab corresponding to the set of 2000 vortices. Due to the long range current distribution of a single vortex, the "vortex-field" in the center of the slab is higher than at the edge of the slab. Only for slabs with a diameter much larger than A, the field in the slab is constant except for a region of order A near the edge of the sample. Because the applied field can also penetrate the sample over a distance A, this edge-effect will not be very pronounced. Still, one has to be aware of field variations due to this effect. To make a direct comparison to the N M R data on 2°STl in T12Ba2CaCu2Oa by Jol et al. and Moonen et al. [14,15], we assume that we can image the T12BaECaCu2Os superconductor as a set of uncoupled planes of CuO2 double layers. The validity of

D. Reefman, H.B. Brom /Numerical simulations on 2D vortex lattices

0.0025 0.0020 0.0015 0.0010

--L_,

j 0.0005 0 -t6

I

-12

I

I

I

I

-8

-'4

0

q

,

I

I

8

12

16

(a,) Fig. 6. Cross section of the vortex-induced field in a slab with 2000 vortices in a field of 5 T. Because the sample size is m u c h less than A, reaches very low values only. If the penetration of the external field is taken into account, is somewhat lower than 5 T. The "spikes" due to the vortices remain, however.

1.|

i

i

i

0.9 ~o

~- 0.7 8 .o

0.5

v

0.3 0.1 -0.1 t~.980

i

i

q. 985

tt. 990

0,. 995

5. 000

219

appeared [26,27] and the shape of the field distribution is more or less Gaussian. For T = 0 K, this shape represents the line shape that would be obtained in an N M R experiment. For higher temperatures, the motion of the vortices leads to motional narrowing of the absorption line. Also, the line shape changes from Gaussian-like to more Lorentzian-like, which is in agreement with the careful lineshape analysis by Moonen et al. [ 15 ]. In fig. 8 we have plotted both the linewidth for the 2°5T1 resonance in T I 2 B a 2 C a C u E O s from our calculations and from the most recent experiment (insert; B perpendicular to the layers). We note that the normal state linewidth (about 80 kHz) is not included in our calculations. For the conversion of field distribution to distribution in resonance frequencies we have used the gyromagnetic ratio 7 for 2°ST1 (0.1977('H)). The temperature where according to the calculations line broadening occurs does not correspond with the temperature where the diffusion constant starts to increase faster than linear. The vortex mobility is large enough already below this temperature to reduce the linewidth. Although qualitatively simulation and experiment agree rather well, the temperature where drastic line broadening sets in is lower for the experiment than for the calculation. This may be either due to a possible size effect in the simulations, or to the effect of pinning centers, which is completely neglected in our calculation. For weak pinning it is known [5 ] that

B (T)

Fig. 7. Magnetic field density function Y ( B ) as function of internal field for a sample large compared to A (effects due to the edge negligible). Except for the tail at the high field side the field density function is almost Gaussian. For low temperatures this profile is identical to the N M R lineshape.

300

'

250 -,,"

350 300

'

' ,

,'

,

'

250 -'~

\

o

o o - %,

x

this assumption is discussed by Clem [ 13 ]; at very low temperatures the layers are decoupled already• The in-plane inter-vortex interaction however, falls off more rapidly as compared to the purely two-dimensional system. In fig. 7 we plot the field density function sff(B) for a slab with a radius of 1 ram, where we have used g(r) as obtained from the simulations with Nv = 2000 to generate the distribution of vortices. Due to the disorder, peaks arising from saddle points have dis-

,50 ~

150-

='=="= o .

tOO ~

S

.

" ~

i

o

.

O

50

0

.

I

20

i

I

tto

,:

.

oo

L

,

,

50

tO0

150

i

I

6o T (K)

i

o

o

200

251

~-"?'-=----,--=~

i

8o

~oo

Fig. 8. N M R linewidth for 2°STl in Tl2Ba2CaCu2Os as calculated. The experimental data (B perpendicular to the layers) are shown in the insert for comparison [ 15 ].

220

D. Reefman, H.B. Brom /Numerical simulations on 2D vortex lattices

the p i n n i n g centers may reduce the temperature at which melting occurs. Also, the linewidth starts to increase already significantly just below To. Most likely, this is cased by the effects of stronger p i n n i n g centers.

4. Conclusions We have shown that for a two-dimensional slab the diffusion constant of the vortices deviates from zero already at low temperatures, a n d increases roughly linear with temperature. For still higher temperatures, the diffusion constant increases much more rapidly. We associate this behavior with the dislocation-mediated melting. For a slab characterised by parameters that are representative for a single CuO2 double layer in TI2Ba2CaCu208, the calculated melting temperature is about 50 K, which is still substantially higher than the analytical result (30 K ) obtained by Fisher [ 17 ]. By assuming that the coupling between these double layers is negligible, comparison of N M R lineshapes of 2°5T1 in Tl2Ba2CaCu2Os with line shapes calculated in our model is possible. We find, that the calculated and experimental line shapes are in qualitative agreement. Differences can be attributed to the neglect of p i n n i n g in our model, a n d possibly, to size effects in the calculation. The calculation presented here may also be relevant for N M R experiments where the external field is oriented parallel to the layers, because, due to the high anisotropy, a slight misorientation of the layers results in a situation very similar to B perpendicular to the layers [ 28,29 ]. Although we focussed on the N M R lineshape only, we will extend the calculations to include the effect of vortex m o t i o n on the spin-lattice relaxation time Tl, a n d the s p i n - s p i n relaxation time 7"2.

Acknowledgements We would like to thank Peter H. Kes, Cees v.d. Beck and Jaco T. M o o n e n for enlightening discus-

sions, a n d M. Elout for the use of the IBM RS6000 workstation.

References [ 1] H.-R. Ma and S.T. Chui, Phys. Rev. Lett. 67 ( 1991 ) 505. [2] L.I. Glazman and A.E. Koshelev, Phys. Rev. B 43 (1991) 2835. [ 3 ] Daemen and Gubernatis, Phys. Rev. B 43 ( 1991 ) 2625. [4] A. Houghton, R.A. Pelcovitsand A. Sodbe, Phys. Rev. B 40 (1989) 6763. [5] H.J. Jensen, A. Brass, A.-C. Shi and A.J. Berlinsky, Phys. Rev. B 41 (1990) 6394. [6] A. Brass and H.J. Jensen, Phys. Rev. B 39 (1989) 9587. [ 7 ] Y. Enomoto, R. Kato, K. Katsumi and S. Maekawa, to appear in Physica C. [8 ] E.I. Kornlov and A.A. Litvin, Z. Phys. B 84 ( 1991 ) 3. [9] J. Pearl, Low Temperature Physics LT19, eds. J.G. Daunt, D.O. Edwards, F.J. Milford and M. Yaqub (Plenum, NY, 1965), Part A, p. 556. [ 10] J. Pearl, Appl. Phys. Lett. 5 (1964) 65. [ 11 ] J.M. Kosterlitzand D.J. Thouless,J. Phys. C 6 (1973) 1181. [ 12] W.F. Brinkman,D.S. Fisher and D.E. Moncton, Science217 (1982) 693. [ 13] J.R. Clem, Phys. Rev. B 43 ( 1991 ) 7837. [ 14] J.C. Jol, D. Reefman,H.B. Brom,T. Zetterer, D. Hahn, H.H. Otto and K.F. Renk, Physica C 175 ( 1991 ) 12. [ 15] J.T. Moonen, D. Reefman and H.B. Brnm, in preparation. [ 16] H.B. Brom and H. Alloul, Physica C 177 ( 1991 ) 297. [ 17] D.S. Fisher, Phys. Rev. B 22 (1980) 1190. [18] B. Rakvin, T.A. Mahl, A.S. Bhalla, Z.Z. Sheng and N.S. Dalai, Phys. Rev. B 41 (1990) 769. [ 19] A. Abragam, Principles of Nuclear Magnetism (Clarendon, Oxford, 1986). [20] J. Bardeenand M.J. Stephen, Phys. Rev. 140 (1965) A1197. [21 ] L.E. Reichl,A ModernCourse in StatisticalPhysics (Edward Arnold, Great Britain, 1987). [22] M.J.L. Sangster and M. Dixon, Adv. in Phys. 25 (1976) 247. [23] A.L. Fener and P.C. Hohenberg,Phys. Rev. 159 (1967) 159. [24 ] M. Tinkham, Introduction to Superconductivity(McGrawHill, New York, 1975). [25] D.R. Nelson and B.I. Halperin, Phys. Rev. B 19 (1979) 2457. [26] P. Pincus, A.C. Gossard, V. Jaccarino and J.H. Wernick, Phys. Lett. 13 (1964) 21. [27] E.H. Brandt, J. Low. Temp. Phys. 73 (1988) 355. [28] S.L. Thiemann, Z. Radovi~ and V.G. Kogan, Phys. Rev. B 39 (1989) 11406. [29 ] P.H. Kes, V.M. Vinokur and C.J.v.d. Beck, Phys. Rev. Lett. 64 (1990) 1063.