Monte Carlo study of the local pair superconductor

PHYSICA

Physica C 199 (1992) 403-413 North-Holland

Monte Carlo study of the local pair superconductor D. Reefman, S.I. M u k h i n i and L.J. de Jongh Kamerlingh Onnes Laboratory, Leiden University P.O. Box 9506, 2300 RA Leiden, The Netherlands

Received 27 May 1992 Revised manuscript received 25 June 1992

The local Pair (LP) superconductoris investigatedby means of a Monte Carlo (MC) method. The ground state phase-diagram in the density versus repulsion/transfer integral plane is calculated, and compared to analytical (RPA) results. The phase diagram features the superconducting phase, a charge-orderedphase and a region where superconductivity and charge ordering coexist. Also, the finite temperature phase diagram and the temperature dependences of the specific heat and the superconductingorder parameter are presented.

1. Introduction It is now generally accepted that the experimental behaviour of the high-temperature superconductors differs markedly from what is observed in "classical" superconductors which can be described by the BCS-theory. Therefore, several different models have been proposed to describe high-T~ superconductivity. A model which has gained interest in this respect is the local pair (LP) or large negative-U Hubbard approach [ 1 ]. The essential difference of this model with respect to BCS is that the charge carders are supposed to form pairs already above T¢. At To, the pairs Bose-condense and form a superconducting condensate. Owing to the possibility of direct mapping of the negative-U Hubbard model in the large I UI limit onto the anisotropic spin ½Heisenberg model (XXZmodel), many of the properties of the LP superconductor can be inferred from those of 2D and 3D antiferromagnets. On this basis, one of the authors has argued already [2 ] that many of the experimental observations might be explained at least qualitatively by this model. To provide a more rigorous basis for the comparison of the LP model with experOn leave of absence from: Moscow Institute of Steel and Alloys, Theoretical Physics Dept., Leninskii prospect 4, 117936 Moscow, Russia.

iment, detailed calculations on the tunneling characteristics and specific heat of the LP superconductor have been performed in refs. [ 3,4 ] within the Random-Phase-Mean-Field Approximation (MFRPA). This method, which was developed in ref. [ 1 ], differs from the Mean Field (MF) approach in that it allows for fluctuations in the order parameter which are not included in MF. The MF-RPA method, however, does not allow for short-range fluctuations above Tc because of the absence of an order parameter above Tc in the MF approach. Also, the effects of the occurrence of a charge-ordered state at higher LP concentrations were not discussed in refs. [ 3,4 ]. From MF-results [ 5 ], it is known that formation of a charge-ordered state suppresses the superconducting transition temperature. The MF approach, however, suffers from the fact that it predicts a superconducting transition temperature Tc which is much too high. In an attempt to meet these deficiencies, a Monte Carlo (MC) method has been developed for this problem. MC studies on quantum Heisenberg or XXZ models have been performed by various authors already [ 6-9 ], but always in constant external magnetic field. Such a procedure cannot be used in this case, because the mapping of the LP model on the XXZ model imposes the constraint that, in order to fix the density of pairs, the longitudinal magnetization rather than the magnetic field should be con-

0921-4534/92/$05.00 © 1992 ElsevierSciencePublishers B.V. All rights reserved.

D. Reefman et aL / Monte Carlostudy of LP superconductor

404

stant. Because of the large computational work involved in quantum MC calculations on the XXZmodel, a semi-classical MC method is developed which is much less expensive in terms of CPU-time. The paper is organized as follows. In section 2, the LP model and the MC method are discussed. In section 3, a comparison of our MC method with quantum MC ( Q M C ) results on the X X Z model is given to provide some indication of the errors introduced by the use of our simplified MC method. Also, the effect of the finite lattice used in the MC calculation is discussed, in connection with a comparison with experimental data on a 2D Heisenberg antiferromagnet. In section 4, the results for the LP model are discussed. Finally, in section 5 some concluding remarks are presented.

2. Model and methods

2. I. The model In the calculations, we base ourselves on the hamiltonian as developed in ref. [ 1 ]. In second quantized notation, the hamiltonian reads:

~LP=½ E Vob~b,b~bj-½ E tob~,b,-ltZ bt, b, • i~j

i~j

i

(1) Here, t 0 and Vo are the LP hopping integral and repulsion for LPs on neighbouring sites i and j, respectively. # is the chemical potential, which determines the density of pairs. The operators b~ and bi, respectively, create and annihilate a local pair on site i. As the on-site local pairs are hard-core bosons, the creation and annihilation operators have to fulfill the following Pauli-like (anti)commutation relations: [b~,b,t]+=l;

[b,,btj]_=O(i~j),

(2)

thus excluding double occupancy of a single site. By introducing the following transformation [ 1,10 ]:

S:[=½(b~, +b,) ; i S~=~(b~,-b~);

z ~ ,, S,=~-bib,"

(3)

where the S~ are the usual spin ½operators, the hamiltonian ( 1 ) may be written up to a constant as:

~LP = ½ E V o S f S ] - ~ ~-, to(SxSy i~j

i#j

+ S f S ] ) - # ~. (1 - S ~ ) .

(4)

i

If we take to=t and Vo= V for nearest neighbours, and zero otherwise, eq. (4) describes the anisotropic Heisenberg ( X X Z ) model, where the chemical potential # plays the role of the external field. The density n of LPs is now related to the average magnetization in the z-direction by the following equation: 1

n= 2

1

~

(Sf>

(5)

where the brackets denote a thermodynamical average. To ensure a fixed density n of LPs, a calculations on eq. (4) have to be performed at constant magnetization instead of constant field. The onset of Bose condensation in this model corresponds to the onset of a non-zero thermodynamical average of the transverse magnetization Mxy:

/

N

l

Therefore, Mxy can be taken as the order parameter for superconductivity (Off-Diagonal-Long-RangeOrder; ODLRO) in the model. On the other hand, the existence of a charge-ordered state (DiagonalLong-Range-Order, D L R O ) is marked by the appearance of antiferromagnetic ordering along the zaxis in the pseudo-spin system. An extension of the model can be made by including possible spatial anisotropy in the values for t and V, which we will refer to as the anisotropic XXZ model. We include only axial anisotropy, thus introducing to= t± if i, j are neighbours in the z-direction, and to= t, if i, j are neighbours in the xy-plane. The same notation will be used for the repulsion parameter V along these directions.

2.2 The Monte Carlo method used As indicated in the precious section, the characteristics of the LP model can be calculated from the anisotropic XXZ model. Although efficient Quantum MC methods exist for dealing with XXZ models, like Handscomb's method, which is elaborated in refs. [ 7,8 ], and the more conventional Suzuki-Trot-

D. Reefman et aL / Monte Carlostudy of LP superconductor ter decomposition [6,9,11 ], the CPU time needed for calculations on a three-dimensional lattice still prohibits the use of those methods for the present problem. To circumvent this complication, we write the state I C ) for the lattice of spins as the product of single spin states: I C> =//~(ai l a > +b, I#> ) •

(7)

Here, the product runs over all lattice sites; [a > denotes the spin-up, Ifl> the spin-down state. The complex parameters ai and bi are only restricted by the normalization conditions: l ai 12+ I b~ 12= 1 .

(8)

We now make the approximation that the value of any operator A can be estimated in the ensemble as follows: Ec< CIAI C > e - a < c ) ~ l c >

Yce_B

=

= Z p(C)A(C).

(9)

C

The distribution p ( C ) is sampled using a standard Metropolis algorithm [ 13 ], giving a series of states I C>. From these, the expectation value of A can be estimated as:

1

~ ~ ~, A(C,)

(10)

i=1

Each sweep through the lattice (updating all a~ and bi) is defined as a single MC step. All results in this paper are obtained from averages of 50 bins of 50 data points each, resulting in a total of 2500 data points. To diminish the correlation between the data points, it was necessary to separate them by a number of MC steps. In the low-temperature regime, over 250 steps between different data points were needed. Owing to these correlation effects, it is impossible to calculate zero-temperature properties with the MC method. This is a feature the algorithm has in common with many other (numerical) methods, like the cluster variation method (CVM) [ 16 ], although the reasons for this shortcoming are different for the two methods. The calculations were performed assuming periodic boundary conditions for the lattice. All calculations up to the 16 X 16 × 8 lattices were performed on either a SUN SPARC workstation or

405

an IBM RS6000 workstation. Simulations on larger lattices were run on a CONVEX-220 and a CRAYYMP4 with veetorized code.

3. Comparison with other results

3.1. QMC results To have an indication of the reliability of the resuits obtained with the MC algorithm described above, we have calculated two cases of the 2-dimensional (2D) XXZ model for wich Quantum Monte Carlo (QMC) simulations have been performed [6,7]. The general hamiltonian used for the XXZ model is: x x z z ~ = - - g E ( SiS} + s i.Vs ) y dI-z~SiS~)

(11)

For J < 0 and A= 1 this hamiltonian describes the isotropic antiferromagnetic Heisenberg (AFH) model, for which the energy and specific heat have been calculated in ref. [7]. For J > 0 and A= - 1, it describes the XXZ model with ferromagnetic and anti-ferromagnetic interactions between, respectively, the planar (x, y) and the z-components of the spins [ 6 ]. In figs. 1 and 2 we show the results on the specific heat obtained with our method and the QMC resuits; both calculations were performed on a 16 × 16 lattice with periodic boundary conditions. Although there is no one-to-one correspondence, the overall agreement is quite good. When the error bars in the QMC calculation on the AFH model are taken into account, no substantial difference is found. For the X X Z model with J > 0 some statistically significant discrepancies are found, although the position of the maximum in the specific heat is in good agreement, and the overall trends seem to conform more to the QMC result than is the case for the AFH model. The correct prediction of the maximum in the specific heat is in agreement with the common lore that, even for the XY model, quantum fluctuations are not too important in the temperature region where the magnetic correlations become strong [ 11 ]. For temperatures lower than that of the maximum (especially in fig. 1 ), our values of CH (T) dearly deviate from the QMC result. The same holds for the energy

D. Reefman et al. /Monte Carlo study of LP superconductor

406

I

I

I ' ' ' ' l

L

/ ~-(t

i

,

I

,5

I

i

,

I

l

I

'

~

1.5

'

I

I

I

I

I

~

2

I

,

L

I

2.5

I

I

3

kBT/J'

Fig. 1. MC results for the specific heat in zero field for the 2D A F H model on a 16× 16 lattice. Units are in J ' to adhere to the notation in ref. [ 7] (J'---J/2). Triangles are the Q M C results from Lee et al., squares results from the present simulation.The drawn curve is a guide to the eye.

6'

•2

.4

.6

.8 kBT/J

1

1.2

1.4

Fig. 2. M C results for the specific heat in zero field for the 2 D X X Z m o d e l ( J > 0, zf = - 1.0) on a 16 × 16 lattice. Triangles are the Q M C results f r o m O k a b e et al. [6 ] ( e r r o r bars are s m a l l e r t h a n s y m b o l size), squares results from the present s i m u l a t i o n . The d r a w n curve is a g u i d e to the eye.

D. Reefman et al. / MonteCarlostudyofLP superconductor

407

U (not shown). This phenomenon is probably linked with the fact that our method does not include quantum fluctuations exactly. We must note, however, that the maximum in the specific heat in these calculations does not correspond to a phase transition, since long-range order cannot occur in this model. However, the specific heat should decay as T -2 for higher temperatures (in the bosonic picture due to the finite width of the boson band), and approach zero for temperatures close to zero. This implies that a maximum in the specific heat should be present anyway, which marks the onset of short-range order and thus occurs at a temperature of about J/kB. Therefore, it may safely be concluded that the present method is able to describe the properties of the anisotropic XXZ model quite well. For the purely isotropic 2D AF Heisenberg model, the method does not seem to work as well. This may be understood from the fact that quantum fluctuations in the 2D isotropic Heisenberg hamiltonian are more severe than for the XXZ model. However, below we shall be interested in the 3D LP model with rather large Vvalues, which system is much closer to the classical limit, since the 3D lattice and, in addition, strong Ising-type anisotropy is considered. Therefore, the results of the simulation for this model should be quite trustworthy, except of course at the lowest temperatures.

peratures higher than T~ 1.2J/kB the HTS expansions are reliable and in good accord with the experimental data; the value of J (the only free parameter) resulting from the fit being in good agreement with other independent determinations [ 12 ]. The discrepancy between the MC data and the experimental results is likely to be a result of finite size effects. If the entropy corresponding to the specific heat is calculated according to:

3.2. Experimental results and HTS expansion

4. Results

In fig. 3 we have depicted experimental results on the 2D Heisenberg antiferromagnet Cu(C6H5NO)6(BF4)2 [12]. Included in the same figure are the results from high-temperature series (HTS) expansion and the MC data on the AF Heisenberg model presented above. Although the present MC results and the QMC results appear to be mutually in reasonable agreement, it should be noted that both MC results for the F H model fall substantially below the experimental results, as well as the HTS expansion prediction. The HTS expansions are in the parameter kBT/J and are in principle exact, but since only a limited number of terms in the expansion can be calculated, the predictions based on HTS expansions become uncertain for low temperatures (below T ~ J / k a ) . Nevertheless, for tem-

4.1. Ground state phase diagram

oo

0

then for the experimental results the correct value of kB In 2 is obtained, whereas the entropy from the MC results only amounts to 0.84ks In 2. This difference can be accounted for by the finite size of the MC system and the ensuing (periodic) boundary conditions, by which a gap is introduced in the excitation (spin-wave) spectrum. The resulting reduction in the magnetic entropy is of the order of the number of spins at two edges divided by the total number of spins. For the 16 × 16 lattice used, this amounts to a decrease of the total entropy by a factor of 0.875, in good agreement with the observed discrepancy. This illustrates that the finite size of the lattice may have important effects on the thermodynamic properties.

The LP model exhibits a rather rich phase diagram. From mean field and MF-RPA calculations it was shown [ l ] that the T= 0 phase diagram features a region for low values of the density n of LPs, where for all values of the repulsion V superconductivity is observed. For larger densities, the system simultaneously shows superconductivity and a charge-ordered state (CO). This part of the phase diagram disappears for very low values of V. Only for exact half filling (n = 0.5 ) does superconductivity disappear for all values of V/t> 1 and only the CO remain. In fig. 4 we depict our MC results for the isotropic LP model. Simulations were performed on a 16× 16× 16 lattice, at a temperature T=O.OO5t/ks,

408

D. Reefman et al. /Monte Carlo study of LP superconductor

' ' ' ' l ' ' ' ' l ' ' ' '

I ' ' ' ' l ' ' ' ' l

....

p./

.2

' ' " ' I l l

0

I'''

... t l l l l , , I , ,

.5

1

i l l

l * , , I , , , , l l l l ,

1.5

2

2.5

3

3.5

kRT/J'

Fig. 3. Specificheat for Cu (CeHsNO)6(BF4)2 (circles,from ref. [ 12 ] ) and high-temperatureseriesexpansion (drawn line). The squares are the present MC results, with a guide to the eye.

CO

il il il

SS + CO

Z'*..

*'\

",,

,, ......

c'

I ' '

~" "~ ~(RPA)

T

I ' ' ' I ' ' , I , I I I I I I I I , I I ,

2

4

6

8

10

12

, , I * i

14

i I , , i

16

I , , , I i i i-

18

20

22

V/t

Fig. 4. Low-temperaturephase diagram of 3D LP system in the n - V/t plane at T= 0.005t. SS:superconductingstate; CO:charge-ordered state; SSCO:mixedphase (both superconductivity and charge ordering). The results from the present MC calculation are indicated by squares (CO/SSCO boundary) and triangles (SSCO/SS boundary). Dashed and dotted lines: RPA and MF results, respectively [ 1]. which is the lowest temperature that could be reached with our method. To check for finite size effects, a few points for a 32 × 32 X 24 lattice were calculated, a n d f o u n d to reproduce the results for the

16 × 16 × 16 lattice within the errorbar. The existence of a charge-ordered state was tested by calculating the Fourier transform of the z-component o f the pseudo-spins for each layer:

D. Reefman et al. I Monte Carlo study of LP superconductor

1

z

S~(q) = N I~ St.,, exp i( lqx + mqy) .

(13)

Here n labels the index of the layer in the z-direction, and I and m run over the sites in the plane. The existence of a Fourier component for q # (0, 0) signals the presence of a charge-ordered state. The staggered magnetization Ms:

Ms" -- ~, ( - 1 )t+ ms~mn ' I,",

409

tion occurs [ 16 ]. In our calculation no phase separation is observed, although we did not restrict the system to a homogeneous magnetization. To make a more definite statement possible, we repeated an MC calculation on a 64× 64X 32 lattice, for which also no phase separation was observed.

4.2. Finite temperature phase diagram

(14)

where l, m, n label the lattice position of the spin in the three-dimensional lattice, is in fact the Fourier component of Sz for q= (~, ~). By calculating only this Fourier component, however, one might overlook possible charge-ordering with wave vectors which are different from (~, n), which was discussed in ref. [14]. From fig. 4, it can be seen that, grosso modo, the MF-RPA results and the MC results coincide. The MF-RPA result for the phase boundary between superconducting (SC) and the non-ordered (NO) phases for zero temperature is shown by a dashed line. A few differences however are evident. First, from our calculation, for a V/t ratio of 1.0 a mixed state is already present. Secondly, in our phase diagram there exists a phase for V/t> 2 for which a pure charge-ordered state is observed, not only for n = 0.5, but also for values slightly lower than 0.5. This, however, may be an artifact, because Tc may be lower in this part of the phase diagram than 0.005t/kB. In the next section, we will indeed show that Tc drops fast for higher LP concentrations. We should note that charge ordering below half-filling cannot be characterised by a single wave-vector of (~, n), but should be described by several wave-vectors (close to (~, n ) ) . Lastly, the boundary between the superconducting region and mixed phase region for high values of V/t settles at a limiting density n of roughly 0.075, which is somewhat higher than the 0.0675 predicted by MF-RPA [ 1 ]. This fact is in line with the observation that MF-RPA in turn predicts for this boundary a higher limiting value of n than mean-field theory (which in fact predicts the limiting value of n to be zero). The nature of the mixed CO/SC phase has gained a lot of interest. For example, it is not clear whether the mixed phase really exists or that phase-separa-

In fig. 5 the finite temperature phase-diagram is shown, for V=4t and (in part) for V= 15t, for the 3D isotropic LP superconductor. Also shown is the MF-RPA result obtained in refs. [ 4,17 ] for the phase ,,,i,l,,i,,,,i,,,,I,,~

+ +

+

I ,

~0

,

,

I

.1

,

,

~ ~ I

.2

~ ,

,

~ libel,

I

.3

.4

,

i

i

! .5

Fig. 5. Finite temperature phase diagram. Squares represent the data for V=4.0t, triangles represent the data for V= 15.0t. The phase boundary between the superconducting phase and non-ordered (metallic) phase, as well as the boundary between mixed phase and charge-ordered phase, is indicated by filled symbols. Open symbols represent the boundary between charge-ordered phase and non-ordered phase. The drawn line is the MF-RPA result (eq. (15) ) for the phase boundary between superconducting phase and metallic phase.

410

D. Reefman et al. /Monte Carlo study of LP superconductor

boundary separating the superconducting phase from the metallic phase: Tc ~ tn 2/3 .

(15)

a bit concave, in disagreement with the MF expression for this phase-boundary which predicts it to be convex: Tco=2n(l-n)V.

(17)

It can be seen from fig. 5 that, for low values of the density n where the system is far from charge order, the agreement with MF-RPA is good; values of T¢ obtained in this paper are only 10%-15% below the MF-RPA values. However, the discrepancy between MC and mean-field in the low-density regime is tremendous. This can be seen from the ratio of the corresponding expressions for the critical temperature [1,51:

Also, the maximum value of Tco, reached for n = 0.5, differs by almost 30% from the value obtained from eq. (17). It should be noted, however, that this is in agreement with results found for the Ising model: here the discrepancy between the mean-field result ( T ¢ = J ) and HTS expansion (T¢~0.75J [15]) is 25%.

T~ F

4.3. Temperature dependence o f the order parameter and specific heat o f the L P system

T eRPA

~

- 1 n2/31nn

(16)

which diverges for low values of n. This is in clear contrast to the usual difference of MF from more exact results for magnetic systems at constant field, where the MF Tc is at most a factor of 2 offfrom the exact result [ 15 ]. As the MC method used is not the correct quantum mechanical approach, the result for Tc presented here is also only an upper bound for the exact To. Also, for low values of n, T¢ seems not to depend heavily on the repulsion parameter V, in concordance with eq. ( 15 ). For high enough n, T~ starts to decrease with increasing pair concentration. In fact, this transition just marks the entrance in the region where the superconductivity coexists with the charge-ordered state, thus decreasing T¢. The large discrepancy with the MF-RPA results in this part of the phase-diagram originates from the fact that RPA is done only on the order parameter for superconductivity. Therefore, the mixed phase cannot be described by the simple MF-RPA results presented in refs. [4,17 ]. For higher densities, stronger repulsion between the local pairs results in the formation of a chargeordered state at lower concentrations so that, after going through a maximum, T~ moves to lower temperatures. This feature in fact remarkably resembles [2] what is observed for the high-T¢ compounds, where a maximum in T~ is also observed as a function of hole concentration. The phase boundary between the charge-ordered and non-ordered phase is indicated by open squares in fig. 5, for a repulsion parameter V = 4 t . For values of n slightly higher than 0.12, this line seems to be

To calculate the temperature dependences of various properties, the chemical potential has to be adjusted for each temperature in order to keep the density constant. In fig. 6 we show the dependence of the pair concentration as a function of the chemical potential, for a fixed temperature of O.Olt/kB, and V = 5 t . The simulations were performed for a 16 × 16 × 16 lattice. From very low densities up to a density of approximately 0.08, n is not strongly dependent on g. In the region 0.08 < n < 0.42, however, a small decrease of I/tl already results in a drastic increase of n. In terms of the magnetic system, this corresponds to the spin-flop phenomenon. Therefore, especially in the region of larger values of n but still not too close to 0.5, special care is needed in adjusting # properly. In fig. 7 we show the temperature dependence of the x y component of the magnetization for a fixed density n=0.071. The results were obtained on a 32 × 32 × 16 lattice. At low temperatures T, the magnetization is slowly decreasing with T, while at Tc the magnetization suddenly vanishes, marking the disappearance of superconductivity. This behaviour is also reflected in the behaviour of the chemical potential # as a function of temperature. Below Tc chemical potential is still varying with temperature (albeit very slow), owing to the fact that for a Bosecondensate of interacting bosons # is not constant. Above T¢ it is decreasing nearly linearly with temperature, conforming to the expressions obtained in ref. [4]. The temperature dependence of the specific heat

D. Reefman et al. / Monte Carlo study o f L P superconductor

I

'

'

'

I

'

'

'

I

'

'

'

411

I

'

| | Ol

T

I

I

I

I

-16

I

I

,

,

I

,

-14

i

i

i

I

-12

l

,

-10

,

I

-8

~2t Fig. 6. Density n o f local pairs as a function o f chemical potential/z, for T= O.02t, V= 5.0t.

,,,,i,,,,i,,,,i,,,,i,,,,i,,,,I,,,,i,,~

•

R

J,,,l,,,,I,,,,I,,,,I,,~,l,,,,I,,,,l~,, 0

.05

.1

.15

.2

I

.25

3

.35

.4

T/t Fig. 7. Transverse magnetization M ~ as a function o f temperature for the 3D LP superconductor. Squares: isotropic 3D system, V= 5.0t. Triangles: anisotropic 3D system, V j . / V I = t ± / t I = 0.2. For both cases n = 0.071.

is shown in fig. 8, for a fixed density of 0.071. Clearly, an ordering anomaly is present in the specific heat at To~ 0.2t/kB, coincident with the appearance of an xy component of the average magnetization. This marks the onset of ODLRO.

Evidently, in the temperature regime below To, the specific heat is not decreasing exponentially, at least not down to the temperatures that were accessible by Monte Carlo simulation. Also, ti can be seen that the peak in the specific heat is quite symmetric around

412

D. Reefman et al. /Monte Carlo study o f LP superconductor

'

'

'

I

'

~

'

I

'

I

I

I

.I

I -

z

]'.]. "

0

.2

.4

I

I

.6

the critical temperature ("2-like"), which is also observed in many experiments on high-T¢ oxides [ 18,19]. This feature was not present in the results of ref. [4], because the short-range order fluctuations were not included properly above T~. Also included in fig. 8 is the specific heat for an anisotropic LP superconductor. At a ratio of t±/ t n= V± / V, = 0.2, T~ drops by only 30%. This is in line with the RPA results [4,17], in which approximation anisotropy only leads to a logarithmic correction to the isotropic case. In fig. 9 the specific heat for a system with a high bipolaron density of approximately 0.3 is shown. Two maxima in the specific heat can now be observed; the strongly peaked maximum at high temperature corresponding to the onset of charge-order, the average of Ms becoming non-zero. At still lower temperatures a second maximum occurs, corresponding

i

I

I

I

.4

.2

keT/t Fig. 8. Temperature dependence o f the specific heat Cv. Squares: isotropic 3D system, V= 5.0t. Triangles: anisotropic 3D system, V L/ VI = t ± / t I = 0.2. For both cases n = 0.071. Lines are guides to the eye.

I

.6

kBT/t Fig. 9. Specific heat for n = 0 . 3 + 0 . 0 2 3 as function of temperature. Isotropic case with V= 5t. The drawn lines are guides to the eye.

to the appearance of superconductivity. Due to the difficulties we had in adjusting # properly for this particular density, large error bars had to be included. For densities larger than 0.4, only a single phase transition (disordered-charge ordered) could be observed in the accessible temperature range ( T> 0.005t).

5. Conclusions It is shown that the critical temperature for LP-superconductivity is strongly dependent on carrier concentration n. For low values of n, Tc increases with n, for higher values it decreases. This feature is in qualitative agreement with experiment and MF predictions [ 5 ]. The reason for this behaviour is the competition between the superconductivity and the charge-order, that suppresses To. For special sets of

D. Reefman et al. / Monte Carlo study of LP superconductor

parameters, a fairly high density and not too small repulsion, the local pair system even has two transitions: from the disordered phase at high temperatures the system first enters the charge-ordered state. On lowering the temperature, the purely charge-ordered state becomes unstable with respect to the mixed phase where charge-order and superconductivity coexist. Also, the specific heat has been calculated and shown to feature a lambda-like anomaly at a temperature corresponding to the onset of superconductivity. This is in agreement with the experimental observation that the specific heat anomaly of the high-To superconductors shows strong short-range order contributions above Tc [ 18,19 ].

Acknowledgements This work is part of the research program of the "Stichting voor Fundamenteel Onderzoek der Materie" (FOM), which is financially supported by the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek" (NWO). It was also sponsored by the "Stichting Nationale Computeffaciliteiten" (National Computing Facilities Foundation, NCF) for the use of supercomputer facilities, with financial support from the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek" (NWO). We would like to thank M. Elout for the use of the IBM RS6000 workstation, and D.J. Bukman for valuable discussions and a critical reading of the manuscript.

413

References [ 1] R. Micnas, J. Ranninger and S. Robaszkiewicz, Rev. Mod. Phys. 62 (1990) 113, and references therein. [ 2 ] L.J. de Jongh, Physica C 161 (1989) 631; idem, J. Less Common Metals 164&165 (1990) 1449. [ 3 ] S.I. Mukhin, D. Reefman and L.J. de Jongh, Physica C 171 (1990) 42. [4] S.I. Mukhin, D. Reefman and L.J. de Jongh, Physica C 174 (1991) 455. [5] S. Robaszkiewicz, R. Micnas and K.A. Chao, Phys. Rev. B 23 (1981) 1447. [6]Y. Okabe and M. Kilmchi, J. Phys. Soc. Jpn. 57 (1988) 4351. [7] D.H. Lee, J.D. Joannopoulos and J.W. Negele, Phys. Rev. B30 (1984) 1599. [8]E. Manousakis, Rev. Mod. Phys. 63 ( 1 9 9 1 ) 1, and references therein. [9] S. Miyashita, J. Phys. Soc. Jpn. 57 (1988) 1934. [10] T. Matsubara and H. Matsuda, Prog. Theor. Phys. bf 16 (1956) 569. [ 11 ] H. de Raedt and A. Lagendijk, Physics Reports 127 (1985) 233. [12] H.A. Algra, L.J. de Jongh and R.L. Carlin, Physica B 93 (1978) 24. [ 13 ] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller and E. Teller, J. Chem. Phys. 21 (1953) 1087. [ 14] K. Kubo and S. Takada, J. Phys. Soc. Jpn. 52 (1983) 2108. [15] M.E. Fisher, Rep. Prog. Phys. 30 (1967) 615; L.J. de Jongh and A.R. Miedema, Adv. Phys. 23 (1974) 1. [ 16] D.J. Bukman, G. An and J.M.J.v. Leeuwen, Phys. Rev. B 43 (1991) 13352. [ 17 ] A.S. Alexandrov, J. Ranninger and S. Robaszkiewicz, Phys. Rev. B 33 (1986) 4526. [ 18] J.W. Loram, ICA. Mirza, W.Y. Liang and J. Osborne, Physica C 162-164 (1989) 498. [ 19] M.B. Salomon, S.E. Inderhees, J.P. Rice and D.M. Ginsberg, Physica A 168 (1990) 283.