30 January1995 PHYSICS LETTERS A
Physics Letters A 197 ( 1995) 345-349
EIs
Quantum charge fluctuations and global superconductivity of Josephson media: path integral Monte Carlo simulations Yu.E. Lozovik a,1, L.M. Pomirchy b a Institute of Spectroscopy, Russian Academy of Sciences, 142092 Troitsk, MoscowRegion, Russian Federation b Institute for High PressurePhysics, 142092 Troitsk, MoscowRegion, Russian Federation
Received 30 October 1994; accepted for publication 16 November 1994 Communicatedby V.M. Agranovich
Abstract
The effect of quantum charge fluctuations on the global superconducting state of Josephson media (arrays of microscopic Josephson junctions, films consisting of tiny granules, etc. ) is considered. Path integral Monte Carlo simulations are used. The system is described by the 2 + 1 XY-model. The dependences of potential, kinetic and total energy, and mean square fluctuations of phases versus the temperature T and the quantum parameter q= h / x / ~ are studied (J is the Josephson constant and C is the granule electrical capacity), The phase transition line Tc= Tc(q) was determined by using the dependences of the helicity modulus versus the temperature and the quantum parameter q. The helicity modulus drops to zero at the phase transition line T= Tc(q). The quasi-3D behaviour in the quantum region of the phase diagram T-q is discussed.
1. Introduction
It is known that in thin metal films the transition to superconducting state occurs according to the Kosterlitz-Thouless scenario [ 1-5 ]. It occurs if the film thickness d is less than the coherent length ~(T) (for uniform film according to the BCS theory ~ ( T ) = h v F / A ( T ) , where vF is the Fermi velocity, A(T) is the superconductive gap). The KosterlitzThouless phase transition relates with the dissociation of pairs of 2D vortices with opposite signs. The low temperature phase has quasi-long-range order of the phases of the order parameter (i.e. power-low decay of correlation fluctuation). At the temperatures T> Tc short-range exponential decay of the correlation function replaces the low temperature behaviour. Granular structure of films may result in some E-mail:
[email protected].
interesting phenomena [5-8 ]. In particular, quantum charge fluctuations in granules in 2D granular film (with thickness d ~ D, D is the granular diameter) result in quasi-three-dimensional (3D) behaviour [ 5,8 ] of the correlation function at distances up to the region crossover R ~ Rc ~ aq/ t ( q = h / @ is the dimensionless quantum parameter characterizing quantum charge fluctuations, t = T / J is the reduced temperature, J is the Josephson constant, and C is the granule capacity). At T < To, in the crossover region R ~ Rc mentioned above, the quasi-3D behaviour is substituted by quasi-long-range order [ 1,3 ] inherent to two-dimensional systems at great distances. Another important effect which has its origin in quantum charge fluctuations is breaking down of superconductivity, even at zero temperatures, at sufficiently large values of the quantum parameter q (i.e. at small size of the granules for granular systems), see, e.g. Refs. [5,6] for the 2D case and Ref. [7] for
0375-9601/95/$09.50 9 1995 ElsevierScience B.V. All rights reserved SSDI0375-9601 (94)00937-6
346
Yu.E. Lozovik, L.M. Pomirchy / PhysicsLetters A 197 (1995) 345- 349
the 3D case). This effect is analogous to "cold" quantum melting of crystals due to zero-point quantum oscillations (see, e.g. Refs. [9,10] ). It is interesting also to describe a possible modification of the mechanism and type of the phase transition for sizes of the system greater than some crossover radius (see, e.g., Ref. [ 8 ] ). Indeed, the partition function for a two-dimensional quantum system is equivalent to the partition function for a classical three-dimensional anisotropic system, namely, for a thick classical anisotropic film with thickness proportional to q/T. The critical behaviour for films with large q may be analogous to the one for a massive anisotropic superconductor. It is very interesting, how and in which region (at real finite systems) the mechanism of Kosterlitz-Thouless phase transition is replaced by behaviour inherent to 3D anisotropic superconductors [ 8,11 ]. It is also interesting to analyze topological excitations responsible for the phase transition in the region of q near qr (in particular, instantons in the 2 + 1 XY-model [ 8 ] ). The phase diagram of this system in the self-consistent harmonical approximation (SCHA) was considered in Ref. [5]. Mean field approximation (MFA) was also used for analogous calculations - see, e.g., Ref. [ 7 ] for the 3D case (MFA does not take into consideration correlations of phases in neighbouring granules). Thus it is interesting to carry on calculations of phase diagram of the system from the first principles (i.e. not using approximations like SCHA, MFA, etc.). This is the purpose of the present article. By path integral Monte Carlo simulations we will calculate the dependences of the physical values which characterize the system (energy, mean square fluctuations, helicity modulus) versus the quantum parameter q and the temperature. From the results of these calculations we will determine the line of the phase transition T~= T~(q).
2. Path integral simulation of the quantum XY-model
The property of granular films in the region of globa! superconductive state is determined mainly by fluctuations of the phases (fluctuations of the gap at small Josephson constants between granules may be neglected in the region near To(q))'. The Hamilto-
nian of the system taking into account quantum charge fluctuations is h2 32 H = - -~-~ ~ ~ { + J ~ ' [ 1 - c o s ( Oi-Oj) ] ,
(1)
where the summation is taken over all sites i (j) and their nearest neighbours (here we take into account only the electrical capacity of the granule, a more complex case will be discussed elsewhere). The canonical ensemble for the quantum system may be described by the Boltzman classical partition function,
Wp(rl, r2, ..., rp,fl)=exp[-flVeff(rl, rE, ..., rp) ] . (2) For the system with Hamiltonian ( 1 ) the effective potential (3) takes the form P
Vefr= E ~ { ( J / P ) [ 1 - c o s ( O ~ - f b y ) ] s= I i,j
+ [P/2J(qfl) z ] (03-0}+1)2},
(3)
where the summation was calculated with respect to i,L s (s is the imaginary time index). In path integral Monte Carlo simulations the average potential energy is equal to
Up(q, f l ) = ( ( J / P )
~ [1-cos(01-0~)])
.
(4)
The average kinetic energy has the form (see, e.g., Ref. [12])
Kp( q, fl) = ( 2NP/2fl -- [P/2Y(qfl) 2] E (q~;--0J+x) 2) ,
(5)
where N is the number of granules. Analysis of the convergence of the results of Monte Carlo calculations of Eqs. ( 5 ), (6) was carried out in Ref. [ 13 ]. The total average energy of the system is
Ep( q, fl) = Up(q, fl) + Kp( q, fl) .
(6)
The important question for computation of the path integral in a discrete point representation is how large should be the value of the parameter P for convergence of thermodynamic averages. The convergence is determined by the discretization parameter 3= flJq/ P. For considered values of the quantum parameter q, temperature T, and for values of P = 10-20 the discretization parameter is z=0.1-0.5. As it was seen from our simulations at such values of the discreti-
Yu.E. Lozovik, L.M. Pomirchy / Physics Letters A 197 (1995) 345-349 r a t i o n p a r a m e t e r the relative variance o f the average potential energy versus P is a p p r o x i m a t e l y equal to 0.01. The analogous value for the average kinetic energy is 0.1, as the kinetic energy ( 6 ) is calculated as the difference between two large values, which are a p p r o x i m a t e l y p r o p o r t i o n a l to P [ 12,13 ]. W h e n the averages are calculated according to Eqs. ( 5 ) and (6), the m e a n square error is t 7 2 ~ O ( P / M ) , where 3 / / i s the total n u m b e r o f M o n t e Carlo steps [ 12 ]. The simulated square lattice c o n t a i n e d N = 5 • 5 sites with phases 0i. We choose the initial state for the Monte Carlo calculation as an equilibrium one for the classical 2D XY-model at the corresponding temperature. The lattice was duplicated P times along the z direction. Initial phase distribution was taken in the interval 0~ = Vt~, where ~ is a uniform r a n d o m value in the interval (0, 1 ), V= ( J / T P ) I/2, q is the maxim u m variance o f phases on one step. Then the quant u m p a r a m e t e r q was increased by a small step A q = 0 . 1 - 0 . 2 . On every step the relaxation o f the system was carried out by the Metropolis technique. The relaxation time was 106P M o n t e Carlo steps. The average values o f kinetic and potential energy were calculated by Eqs. ( 5 ) , ( 6 ) at every step o f the relaxation, that is at every fixed value o f q.
347
1.5
,o
i
0.5
i
t
0
t
t
0.5
i
1.0
Fig. 2. The helicity modulus ? as a function of the quantum parameter q for the temperature T/J= 0.35.
0,8
0,4
1 0
3. Discussion
of results
The dependences o f the average energy versus q are presented in Fig. 1. The dependence o f the potential
'
'
"
1.0'
"
2.0
'
Fig. 3. The temperature dependence of the system energy per site for the quantum parameter q=0.3: ( 1) potential energy U, (2) kinetic energy K, (3) total energy E= U+K.
0.8
1.0
I I I l I
0.4
0.5
I ,
0
i
0
i
1
I
I
i
2 q
Fig. 1. The energy of the system per site as a function of the quantum parameter q for the temperature T/J=0.35: ( 1) potential energy U, (2) kinetic energy K, ( 3 ) total energy E = U+ K.
*
i
0.5
i
T/J
l
l
t
J
1.0
Fig. 4. The temperature dependence of the helicity modulus y for the quantum parameter q = 0.2.
Yu.E. Lozovik, L.M. Pomirchy /Physics Letters A 197 (1995) 345-349
348
8~
6
4
2
0.'5 . . . . V~
).b
Fig. 5. The m e a n square fluctuation of the phases, ~qxy, as function of the q u a n t u m parameter q for the temperature T/J= 0.2.
9(~)
t
0.8
04
o.0
~
lo
1's
r Fig. 6. The correlation function of the phases, g(r), for the system with N= 30• 30 sites at the temperature T/J=0.2 for two values of the quantum parameter q: ( 1) q= 0.3, (2) q = 2.1.
-/3 5,0
O.5
. . 0
1.0
q
. \ ~ ~l~ , 2.0
Fig. 7. The Tc/J-q phase diagram for 2D granular system. The solid line is the calculation in SCHA [ 5]. The crosses represent results of our path integral Monte Carlo simulations.
energy versus q was found to be approximately a linear function. To determine the phase transition point in the classical XY-model the behaviour of the helicity modulus 7 = O2F/Oh2 versus temperature may be used [ 3 ], h is the imposed gradient of the phase and F is the free energy of the system. In Re/'. [ 14 ] the behaviour y (T) was used to determine the phase diagram for the classical XY-model with percolation. In the present article the helicity modulus 7 and its dependence versus the temperature T a n d the quantum parameter were calculated. The characteristic plot y(q) at fixed temperature is presented in Fig. 2.~The helicity modulus 7 as function of q decreases and has a sharp jump to zero at the critical value q = qc (Fig. 2). This point determines the quantum breakdown of phase correlations in the system due to quantum phase fluctuations. The value of the quantum parameter q at which sharp falloff of 7(q) has occurred at fixed temperature T is taken as the point of phase transition to plot the q-T diagram. In some cases another procedure, in comparison with the one described in Path 2 for "melting" of the model system, was used. The initial state was prepared, as described in Path 2 by duplicating the relaxed XY-model. Then the value q increased up to some q less than qc (T) with relaxation at every step dq. To determine the temperature dependence of the thermodynamic averages (Figs. 3, 4) the temperature of the system increased with relaxation for every step A T = 0.1. When the parameter q increases y(T) decreases. At increasing q the jump in the dependence of the helicity modulus 7(T) at the point T = Tc becomes different from the universal j u m p 7(T;- ) / T~- = 2 / r t characteristic of the Kosterlitz-Thouless phase transition for classical two-dimensional systems. Besides the dependence 7(T) at fixed q becomes smoother when q increases. This may be considered as the manifestation of 3D nature of the system when q becomes larger (crossover from 3D to 2D behaviour is discussed in Refs. [5,8]. The temperature of the phase transition was determined from the location of the j u m p of 7(T) to zero. The points of boundary Tc (q) on the phase diagram determined by such a way from the dependences 7(T) are consistent with the points determined by jumps on the curves y(q). The mean square fluctuations of the phases in the
XYplane,
Yu.E. Lozovik, L.M. Pomirchy / Physics Letters A l 97 (1995) 345-349
349
References 1
~r
4 N P ~ [ ~ ( i ' s ) - O ( j , s) 2 ] ~/2,
(7)
were also determined. The dependence of this value versus q at low temperatures is presented in Fig. 5. The value ~Oxycharacterizes the relative diffusion of phases of neighbouring sites during Monte Carlo simulations. The sharp increase of 8~xy is the evidence of phase transition - "quantum melting" which occurs increasing the parameter q at low temperatures (compare Refs. [ 9,10 ] ). Besides the correlation function of phases, g ( r ) = ( c o s [ ~ ( r ) - ~ ( 0 ) ] ) , was computed. In this case the system of larger sizes with N = 30 • 30 was used. The correlation function has a weak dependence on distance as in the case of quasi-3D order and decreases when T or q increase (Fig. 6). The phase diagram resulting from our simulations is shown in Fig. 7 together with the phase boundary To(q) calculated by SCHA in Ref. [ 5 ]. We point out that the results obtained from calculations of the values ~,(T) and 7(q) are in agreement with each other and do not have a too large difference from the resuits obtained by SCHA [6]. Note that the phase transition according to the simulation results occurs slightly before the SCHA stability line. This situation is analogous to Wigner crystal melting (see, Refs. [9,15]).
[ 1 ] J.M. Kosterlitz and D.J. Thouless, J. Phys. C 6 ( 1973 ) 1181. [2] B.I. Halperin and D.R. Nelson, J. Low Temp. Phys. 36 (1979) 599. [3] P. Minnhagen, Rev. Mod. Phys. 59 (1987) 1001; P. Minnhagen and H. Weber, Physica B 152 ( 1988 ) 50. [41 A.T. Fiory, A.F. Hebard and W.I. Glaberson, Phys. Rev. B 28 (1983) 5075; A.F. Hebard and A.T. Fiory, Phys. Rev. Lett. 58 No. 11 (1987) 1131. [5] Yu.E. Lozovik and S.G. Akopov, Solid State Commun. 35 (1980) 693; J. Phys. C/14 No. 2 (1981) L31; 15 (1982) 4403. [6] E. Simanek, Solid State Commun. 31 (1979) 419. [7] K.B. Efetov, Zh. Eksp. Teor. Fiz. 78 (1980) 2017. [8] N.K. Kultanov and Yu.E. Lozovik, Solid State Commun. 88 (1993) 645; to be published. [ 9 ] Yu.E. Lozovik and V.M. Farztdinov, Solid State Commun. 54 (1985) 725; B. Abdulaev, Yu.E. Lozovik and V.M. Farztdinov, Fiz. Tverd. Tela 27 ( 1985 ) 3124. [ 10] Yu.E. Lozovik and V.A. Mandelshtam, Phys. Left. A 165 (1992) 469. [ 11 ] D.R. Nelson and J. Toner, Phys. Rev. B 24 ( 1981 ) 353; V. Kotsuvo and G.A. Williams, Phys. Rev. B 33 (1986) 6106; G.A. Williams, Phys. Rev. Lett. 59 (1987) 1926; S.R. Shenoy, Phys. Rev. B 40 (1989) 5056; B. Chattopadhyay and S.R. Shenoy, Phys. Rev. Lett. 72 (1994) 400. [ 12] M.F. Herman, E3. Bruskin and B.J. Berne, J. Chem. Phys. 76 (1982) 5150. [ 13] Yu.E. Lozovik and V.A. Mandelshtam, Nuovo Cimento B 109 (1994) 575. [14] Yu.E. Lozovik and L.M. Pomirchy, Fiz. Tverd. Tela 35 (1993) 2519; Solid State Commun. 89 No. 2 (1994) 145. [ 15 ] Yu.E. Lozovik, to be published.