Elasticity and melting of vortex crystals in anisotropic superconductors: Beyond the continuum regime

The elastic moduli of vortex crystals in anisotropic superconductors are frequently involved in the investigation of their phase diagram and transport properties. We provide a detailed analysis of the harmonic eigenvalues (normal modes) of the vortex lattice for general values of the magnetic field strength, going beyond the elastic continuum regime. The detailed behavior of these wavevector-dependent eigenvalues within the Brillouin zone (BZ), is compared with several frequently used approximations that we also recalculate. Throughout the BZ, transverse modes are less costly than their longitudinal counterparts, and there is an angular dependence which becomes more marked close to the zone boundary. Based on these results, we propose an analytic correction to the nonlocal continuum formulas which fits quite well the numerical behavior of the eigenvalues in the London regime. We use this approximate expression to calculate thermal fluctuations and the full melting line (according to Lindeman's criterion) for various values of the anisotropy parameter.

The elastic moduli of vortex crystals in anisotropic superconductors are frequently involved in the investigation of their phase diagram and transport properties. We provide a detailed analysis of the harmonic eigenvalues (normal modes) of the vortex lattice for general values of the magnetic field strength, going beyond the elastic continuum regime. The detailed behavior of these wavevectordependent eigenvalues within the Brillouin zone (BZ), is compared with several frequently used approximations that we also recalculate. Throughout the BZ, transverse modes are less costly than their longitudinal counterparts, and there is an angular dependence which becomes more marked close to the zone boundary. Based on these results, we propose an analytic correction to the nonlocal continuum formulas which fits quite well the numerical behavior of the eigenvalues in the London regime. We use this approximate expression to calculate thermal fluctuations and the full melting line (according to Lindeman's criterion) for various values of the anisotropy parameter.

I. INTRODUCTION
Among the fascinating aspects of the new hightemperature superconducting materials, there are many which are due to the rich and complex behavior of vortex lines 1 . These quanta of magnetic flux penetrate the superconductor above a certain threshold value of the external magnetic field H, the so-called lower critical field H c1 (roughly 10 −2 T , and their concentration increases with H up to the upper critical field H c2 (approximately 10 2 T ), above which the material is normal 2,3 (see Figure 1). In a classic work, Abrikosov showed that the minimum energy arrangement of flux lines in a conventional superconductor is a triangular lattice, with a lattice spacing a which varies with the magnetic field strength 4 .
Recently, with the discovery of high temperature cuprate superconductors, there has been a resurgence of interest in the properties of vortex matter. One novel aspect of these materials is the symmetry of the superconducting order parameter. A series of experiments 5 seem to provide evidence consistent with d-wave symmetry, in contrast to conventional superconductors which have s-wave symmetry. In addition to modifying the internal structure of the vortex lines, d-wave symmetry has implications for the global mean-field arrangement of the vortex lattice 6 . In particular, an oblique square lattice appears to replace the standard triangular arrangement, as the most stable configuration in part of the B − T phase diagram 7 . (The new phase appears for high values of the magnetic field, whereas the triangular arrangement is confined to low fields.) Furthermore, the relative flexibility of vortex lines in these materials makes them susceptible to distortions by thermal fluctuations, and other sources of disorder (oxygen impurities, grain boundaries, etc.). The regular lat-tices obtained in mean-field theory are thus distorted, giving rise to a rich variety of vortex-matter phases 8,9 . It is thus necessary to understand the elastic response of vortex lattices to distortions, a subject that has been intensely studied in the context of conventional superconductors. The stability of the triangular lattice against small distortions is guaranteed as long as the characteristic length of the field variations, the so-called penetration depth λ, remains smaller than the size of the system 10 .
(An infinite λ renders the lattice of vortex lines unstable against fluctuations or, more precisely, against shear deformations.) At low temperatures, fluctuations are well described by small corrections to the mean-field results. Although the phenomenology is similar for d-wave superconductors, less is known about the effects of thermal fluctuations and disorder in this case. We shall thus focus on a triangular lattice of vortex lines. The extension of the results to oblique configurations is possible, and should prove interesting.
The elastic properties of the triangular vortex lattice at long distances are characterized by its compressional, shear, and tilt moduli. These moduli are frequently involved in the theoretical and experimental determination of the properties of the material, as for example, the intricate vortex phase diagrams. As such, the derivation of the elastic energy and the elastic moduli has been undertaken in several publications [11][12][13][14][15][16][17][18][19][20] . While it is well known that the elastic moduli of the vortex lattice depend strongly on the magnetic field strength, the values currently presented in the literature are only strictly useful in certain limiting situations. Most frequently, results are obtained in the so-called continuum limit, in which one disregards the discrete nature of the underlying vortex lattice. Obviously, this description is not suitable for the whole range of possible flux-line densities in the mixed state, which extends from H c1 up to the upper critical field H c2 , at which the magnetic field penetrates uniformly into the material. The elastic properties of the vortex lattice, and hence its stability against thermal fluctuations, depend sensitively on the value of the magnetic field in this range. Furthermore, the important fluctuations sometimes occur at short wavelengths, where a simple elastic description may not be appropriate. It is thus worthwhile to obtain the general dependence of the energy cost of harmonic distortions of the vortex lattice. As temperature is increased, thermal fluctuations cause the vortex lattice to melt into a vortex liquid. Several experiments employing quite different techniques have provided firm evidence for such a transition [21][22][23][24][25][26][27][28][29] . In the new cuprate superconductors, the melting transition can occur at temperatures well below the mean field point, so that the Abrikosov lattice is melted over a substantial portion of the phase diagram. Furthermore, the vortex lattice can melt not only by increasing temperature, but also by decreasing the magnetic field to the vicinity of H c1 (T ). In this region of the phase diagram, the concentration of vortices is very dilute; the separation a between neighboring flux lines is larger than the penetration length λ, and the vortex-vortex interaction decays exponentially. As a consequence, the elastic moduli become exponentially small, and correspondingly, thermal fluctuations are greatly enhanced. This behavior gives rise to an interesting reentrant behavior of the melting line at low fields.
In this paper, we take up the computation of the elastic properties of a vortex lattice, in a systematic manner which is not restricted to the most commonly considered continuum limit. Our analysis allows the illustration of some of the subtleties involved in the calculation, and to successfully resolve some of the unclear points that we encountered in our review of the literature, and which may also have puzzled other investigators of the subject. The paper is organized as follows: In Sec. II, we introduce the Ginzburg-Landau Hamiltonian for an anisotropic superconductor. This model is commonly used in the literature as the starting point for studying the effects of fluctuations and disorder. Minimization of this Hamiltonian gives the optimal lattice configuration, around which we then study the cost of distortions in the harmonic approximation. This section also introduces the main parameters and notation used throughout the paper. In Sec. III, we provide an expression for the harmonic kernel in terms of a sum over Bravais lattice vectors or, equivalently, over reciprocal lattice vectors. The normal mode eigenvalues are explicitly written in terms of this kernel. Sec. IV is composed of subsections in which we introduce certain limiting situations, corresponding to a single line, local, local-continuum, and nonlocalcontinuum limits. These limiting cases are often quoted in the literature, and are widely used in studies of vortex matter in high temperature superconductors. Within our framework, we recalculate the analytic values of the elastic eigenvalues in these limits before discussing our more general results in Section V, where the harmonic eigenvalues are numerically evaluated. Our results are graphically presented in several plots, and compared to the limits introduced in Sec. IV, to easily visualize the accuracy of the approximations involved. As a practical illustration of the potential applications of our results, we use the harmonic eigenvalues to calculate the thermal distortions of the vortex lattice in Section VI. The leading contribution to flux-line fluctuations in real space comes from the transverse modes. In conjunction with the Lindeman's criterion, we can then find the full form of the melting line as a function of the magnetic field. This is one of many potential applications that our general analysis makes possible, without the need to extrapolate the elastic behavior of the lattice to regimes beyond their range of validity. Finally, in the last section we summarize our main conclusions.

II. GENERAL FORMULATION
Our starting point is the continuum Ginzburg-Landau free energy for an anisotropic superconductor in the London limit 2 . In this limit, the penetration length λ (∼ 10 3 A) is much larger than the coherence length of the superconductor ξ (∼ 10Å), and fluctuations in the magnitude of the order parameter Ψ o are neglected. The phase degree of freedom is then the only variable that needs to be considered. For the anisotropic superconductors under consideration, this approximation breaks down in a narrow band close to H c2 , where the separation between vortices becomes comparable to ξ. The Ginzburg-Landau Hamiltonian in this limit is given by Here θ(r ⊥ , z) is the phase of the complex order parameter field Ψ(r ⊥ , z) = Ψ o e iθ(r ⊥ ,z) , A(r ⊥ , z) is the magnetic vector potential related to the magnetic induction through b(r ⊥ , z) = ∇ × A(r ⊥ , z), α = (hΨ o ) 2 /2m, φ o = hc/2e is the flux quantum, and ǫ 2 = m/M is the usual anisotropy parameter, with m and M being the effective masses in the Cu-O planes of the material, and along the perpendicular c axis chosen to coincide with the z direction, respectively. The external magnetic field H is also oriented parallel to the z axis H = Hê z . Note that ∇ ⊥ , r ⊥ , and A ⊥ , denote the planar components of ∇, r, and A, respectively. Consider a configuration with N vortex lines, and decompose the corresponding phase field into two parts: , which represents the singular phase due to the N vortices passing through points R n (z) = R n,⊥ (z) + zẑ (n = 1, . . . , N ) on independent planes at each z; and θ r , a regular phase field accounting for the relaxation due to the couplings between the planes. By construction, each θ s n is the solution of a two-dimensional problem with the circulation constraint C dθ s n = 2π on any closed circuit C around the n-th vortex. We have chosen the coordinate z to parametrize the trajectories of the different flux lines. Nonetheless, one has to bear in mind that the results should be invariant under an arbitrary reparametrization.
Variations of Eq. (1) with respect to the phase θ r , and the vector potential A, provide the differential equations for these quantities, whose solutions minimize the energy in (1). After considering the Coulomb gauge ∇ · A = 0, these equations read where we have introduced the squared penetration depth in a plane perpendicular to the z axis, λ 2 = φ 2 o /(32π 3 α) = mc 2 /(16πe 2 Ψ 2 o ), and ∆ ⊥ stands for the in-plane component of the Laplacian operator. The coupled equations (2) and (3) are easily solved in Fourier space to give θ r (k) = ik z (λ 2 k 2 + 1) k 2 (λ 2 c q 2 + λ 2 k 2 z + 1) where k = q + k zẑ , λ c = λ/ǫ is the penetration length along the c-axis of the superconductor, and we have introduced the function which represent the Fourier transform of the vortex function ∂ z θ s n . Because of its singular nature, the Fourier transform of ∇ ⊥ θ s n , actually has a component which is perpendicular to q. We thus define the transversal and longitudinal components of the in-plane projection of the gauge field, as A T ⊥ = (I −qq) · A ⊥ , and A L ⊥ =qq · A ⊥ respectively. The former is obtained from Eq. (4) as while the latter is obtained from the Coulomb gauge condition as For a given distortion of the flux-lattice, the above solutions provide the form of the phase and gauge fields that minimize Eq. (1). The next step is to evaluate the energy cost of such distortions by substituting these solutions into the Ginzburg-Landau free energy, resulting in Let us now assume that the position vector R n,⊥ (z) is the sum of a perfect lattice vector, plus a small displacement field due to fluctuations, i.e. R n,⊥ (z) = R o n + u n (z). Up to second order in the displacements u n (z), the energy cost can then be expressed as The first term H o , gives the free energy of an array of N straight flux lines oriented parallel to the external field, and located at positions R o n , Here L is the sample thickness, H c1 = φ o /(4πλ 2 ) ln(κ) is the lower critical field, κ = λ/ξ is the Ginzburg-Landau parameter, H is the magnetic field strength, ǫ o = φ 2 o /(4πλ) 2 is an interaction energy per unit length, K o (x) is the modified Bessel function of zeroth order, and d nm = R o n − R o m are the relative position vectors of any pair of flux lines. On a perfect triangular lattice, these vectors are of the form d nm = a(nê 1 + mê 2 ), (n, m = 0, ±1, ±2, . . .), with a the lattice spacing,ê 1 oriented, for instance, along the x axis,ê 1 =x, andê 2 = cos(π/3)x+sin(π/3)ŷ (see Fig. 2). The modulus of one of these vectors is then given by d nm = a(n 2 + m 2 + nm) 1/2 .
The first term on the right hand side of Eq. (13) represents the energy cost of a single vortex line in a type II superconductor, times the number of lines N . The last term is due to the interactions among flux lines, which naturally depend on the interline separations. The penetration length λ, sets the extent of the interaction potential K o (x), which diverges logarithmically at short distances, and decays exponentially for x ≫ 1.
The quantity ∆H represents the harmonic contribution of fluctuations to the free energy. After writing the displacement fields u n (z) in terms of Fourier modes (14) and taking advantage of the translational symmetry of the lattice, the energy ∆H can be written as where the index BZ indicates that the integration is performed over the first Brillouin zone in reciprocal space. All the relevant information is therefore contained in the harmonic kernel M αβ (Q, k z ). Formally, we have to find the eigenvalues and eigenvectors of this matrix, and we can then calculate the extent of fluctuations in real space, the fluctuation corrections to the free energy, and other relevant quantities. In Fourier space, the eigenvectors of M αβ (Q, k z ) are N longitudinal modes, u L (Q, k z ) =Q · u(Q, k z ), and N transversal modes, u T (Q, k z ) = |Q × u(Q, k z )|, with corresponding eigenvalues Λ L (Q, k z ) and Λ T (Q, k z ). In practice, it turns out that the analytic expressions for these eigenvalues is rather complex. That is why in the literature only certain limiting regimes are usually treated. We shall describe these limits later on in Sec. IV, and then go on to discuss their accuracy in comparison to our more general results.

III. LONGITUDINAL & TRANSVERSAL MODES
The longitudinal elastic eigenvalue Λ L (Q, k z ) is typically expressed in the literature in terms of the socalled compressional, c 11 (Q, k z ), and tilt, c 44 (Q, k z ), elastic moduli as c 11 (Q, k z )Q 2 + c 44 (Q, k z )k 2 z . (In general, however, we shall see that this decomposition may be inaccurate.) In terms of the matrix M αβ (Q, k z ), the longitudinal eigenvalue is given by where, as usual, a repeated index is summed over. On the other hand, the transversal eigenvalue Λ T (Q, k z ), commonly written in the literature in terms of the shear, c 66 (Q, k z ), and tilt moduli as (17) There are two alternative ways of expressing the interaction kernel M αβ (Q, k z ): in terms of a sum over the Bravais lattice vectors, or as a sum over the reciprocal lattice vectors. The former yields where we have introduced the areal density of flux lines n = N/A. For compactness of notation, let us also introduce the dimensionless variables x = λ 2 k 2 z + 1, and x c = ǫ x. In Eq. (18), the quantity nǫ o E(k z )/λ 2 represents the line tension of each vortex, which, in general, contains contributions from both the Josephson coupling between the different Cu-O layers in the material, and the magnetic interlayer interactions, as The interaction kernel R αβ (d, k z ) is given by Note that we use the symbol d = |d| to indicate the distance between a pair of lattice points. The sum in Eq. (18)  The matrix R(d, k z ) depends on the different relative position vectors d = d nm of the perfect triangular lattice. When summing over all lattice vectors with a common origin at a point P , we can take advantage of the lattice symmetries. For example, the sum over products of an odd number of vector components vanishes, as in On the other hand, the sums involving products of an even number of components are nonvanishing and constrained by the symmetries; for example, (Higher order terms have a more complex structure.) The sums d =0 on the right hand side of Eqs. (22), are equivalent to d =0 g(d), with a degeneracy factor g(d), which counts the number of vectors whose modulus is d.
Equivalently, we can express the kernel as a sum over reciprocal lattice vectors G. The same symmetries (21)- (22) are, of course, valid for these vectors, which, according to our choice of d nm , are given by 3a is the reciprocal lattice spacing, and the unit vectors areĝ 1 = sin(π/3)x − cos(π/3)ŷ, andĝ 2 =ŷ. In terms of these vectors, the interaction kernel reads The symbol S indicates symmetrization with respect to The second term on the right-hand side of the expression is included to properly account for the particular case n = m and z = z ′ in Eq. (11).) In addition, from Eq. (11) we obtain Our goal is to provide accurate values of the harmonic eigenvalues as a function of Q and k z within the BZ, and for different concentrations of flux lines. The dimensionless parameter a/λ will be used as the indicator of the areal density of flux lines. A small value of a/λ, corresponds to a highly dense regime with strong and non-local interactions among the flux lines. On the other hand, for a/λ ≫ 1 the concentration of flux lines is very low, and interactions among them very weak. In this dilute limit, the elastic behavior of the lattice reflects to properties of a single flux line. We thus evaluate numerically the harmonic eigenvalues as a function a/λ throughout the BZ. The equivalence of expressions (18) and (23) renders the use of one or the other a matter of convenience. For instance, for very high areal densities of vortices n, it is common to disregard the discreteness of the underlying arrangement and approximate the sums over lattice positions by integrals. This so-called continuum limit is very often used in the literature [17][18][19][20] . The elastic moduli which follow from this approximation can be read directly from Eq. (24), when only the reciprocal lattice vector G = 0 is taken into account. We shall comment on the accuracy of this limit in the following sections.

IV. LIMITING REGIMES
In this section, we first present some special cases previously discussed in the literature. By comparing these limits to our numerical results, one can then see their range of validity and the accuracy of the approximations involved in their formulation. We shall also emphasize the roles played by the interline distance, the penetration length λ, and the symmetries of the triangular lattice.

B. Elastic moduli
The energy cost of long wavelength distortions is described by an elastic theory, whose form is constrained by symmetries. The triangular lattice is isotropic, and governed by the compressional, c 11 , shear, c 66 , and tilt, c 44 , moduli. In terms of the Fourier modes, the elastic energy is (25) and the harmonic eigenvalues have the simple forms Formally, Eq. (15) together with (18) or (23) provide the harmonic eigenvalues throughout the Brillouin zone. Naturally, the elastic limit is regained by expanding these results up to second order in Q and k z . Neglecting the higher order terms is referred to as the local limit, as it is usually obtained by including only short-range interactions. Sometimes, non-local elastic moduli are introduced which depend on Q. This is not always useful, as it constrains the form of the harmonic eigenvalues as in Eqs. (26), whereas symmetry allows higher order powers of Q to appear in other forms. After expanding the cosine, and using the symmetry properties of a triangular lattice in Eqs. (21)-(22), as well as certain relationship among the modified Bessel functions 30 , we arrive at Comparing this expression with Eq.(25) allows us to identify the compressional modulus with the shear modulus with and the tilt modulus with Note that the mean-field lattice spacing a, is related to the strength of the applied field H through the relationship which is obtained from ∂H o /∂a| a = 0. Then, by comparison, the local tilt modulus for an isotropic material (ǫ = 1) can be written as c 44 = nφ o H/4π, in agreement with the result expected for the local tilt modulus of a rotationally invariant superconductor 13 . Strictly speaking, in Eq. (32) there is an extra factor of 1/2, coming from the local line tension. However, as pointed out by Fisher 19 , one could chose a short distance cutoff ξ 31 inside ln(λ/ξ) to reproduce exactly the local isotropic limit. In Figs. 3a) and 3b), we depict the functions F c and F s , for different values of the dimensionless quantity d/λ, to illustrate some subtleties associated with the estimation of the shear modulus, often forgotten or misunderstood in the literature. As shown in Fig. 3b), for small values of d/λ the function F s is negative, whereas if d/λ > 2 it becomes positive and then decays exponentially to zero. The implications of this functional form for the shear modulus are as follows: As indicated before, the range of interactions among vortex lines is determined by the penetration length λ. For a dilute vortex lattice whose lattice spacing a is comparable or even greater than λ, almost all the terms in the sum in Eq. (30) are positive, and as they rapidly decay to zero, only the first few lattice vectors for neighboring sites are needed to calculate the shear modulus. On the other hand, in a dense system there are many neighboring vortices within the interaction range λ, contributing a negative amount to the sum in Eq. (30). It is then clearly not sufficient to consider only interactions among a few neighboring lines, and to account for the stability of a dense lattice against shear deformations, we need to sum over many lattice vectors (finally giving rise to a positive shear modulus). In Figs. 3a) and 3b), we also indicate with a dashed line the limit λ → ∞ (for a finite value of d). In this case there is a logarithmic interaction among the vortices, as in the two dimensional Coulomb gas 32 . Note that in this limit, all factors in the sum for the shear modulus are negative, and the lattice is unstable to shear deformation. (This is a manifestation of the more general instability of a system of point charges in the absence of external potentials.) Figure 4 depicts the compressional, shear, and tilt moduli as a function of a/λ, obtained by numerical evaluation of the sums in Eqs. (28), (30), and (32), respectively. These moduli have been normalized by their respective values in the local continuum limit, which for the long wavelength compressional modulus is defined as where we have replaced the sums by integrals, after first adding and susbtracting the d → 0 element explicitly, and introducing the appropriate unit area A pc = √ 3a 2 /2, i.e. the area of the primitive cell. We have also defined the average magnetic induction B = nφ o , with n = N/A = 1/A pc . Similarly, the long wavelength shear modulus is and, for the tilt modulus we obtaiñ Note that for the latter, the term d → 0 is already included in Eq. (32). Similar results for the continuum local limit have been previously reported in the literature 1,12,14 . In Ref. 14, Kogan and Campbell calculated the shear modulus of a flux lattice in different situations. Essentially, they provided an expression for the shear modulus c 66 in terms of a sum over components of the reciprocal lattice vectors. Because of the symmetry relations of Eq.(22) quoted in the previous section, all the terms included in their sums cancel each other out, and should yield a shear modulus equal to zero. Kogan and Campbell numerically evaluated this sum over reciprocal lattice vectors, and reported a value close to Eq. (36). They argued that there is a contribution due to the specific form of the cut-off at short distances 31 . However, in view of the above derivation of Eq.(36), it is clear that cutoffs are not relevant to the local continuum value of the shear modulus.
As soon as the lattice spacing is greater than the penetration length, a > λ, all the compressional, shear, and tilt moduli deviate strongly from the continuum results. The compressional and shear moduli then decay exponentially to zero, while the tilt modulus attains a constant value of nǫ o (ǫ 2 ln(κ/ǫ) + 1/2). (This is why when normalized by its continuum value of Eq.(37) the ratio grows as (a/λ) 2 .) As we shall demonstrate in Sec.(VI), under these circumstances the lattice again becomes rather soft and liable to melting, as its fluctuations are greatly enhanced for large interline separations.

C. (Non-local) Continuum limit
The continuum limit refers to a situation in which the concentration of flux lines is high. In this limit, one disregards the underlying triangular lattice arrangement and treats the array of vortices as a continuous medium 17-20 . In the continuum limit, the sum of Eq. (18) over Bravais lattice vectors is approximated by an integral, or equivalently, in the sum over reciprocal lattice vectors in Eq. (23) only the vector G = 0 is considered. As a result, the elastic kernel reduces to The last term in this equation is a rather intricate integral, which, on the other hand, can be easily evaluated in the local approximation, i.e. when keeping terms up to Q 2 and k 2 z only. With this approximation, we can rewrite Eq. (38) as By comparing with Eq. (25), we can now directly read off the nonlocal continuum (again, up to the approximation carried out to evaluate the integral) values of the elastic moduli as Note that, up to second order in Q 2 and k 2 z , Eqs. (40) through (42) reproduce the local limits calculated in the previous section from the real space expression of the interaction kernel.

V. NUMERICAL RESULTS
In this section, we present and discuss our results from numerical evaluations of Eqs. (16)-(17) using, in particular, the real space expression of the interaction matrix in Eq. (18). In most cases, we selected a Ginzburg-Landau parameter κ equal to 100, and an anisotropy ratio of ǫ = 0.1. However, later in this section, we present some results for different values of ǫ, which may be more appropriate to describe materials with stronger anisotropy. For the sake of simplicity, we introduce the dimensionless eigenvalues Λ * L = Λ L /(nǫ o /λ 2 ) and Λ * T = Λ T /(nǫ o /λ 2 ), and discuss their behavior throughout this section.

A. Angular dependence
Because of the point group symmetries of the triangular lattice, it is sufficient to consider only the irreducible Brillouin zone (IBZ), indicated in black in Fig. 5. In particular, we have chosen five paths corresponding to values of the angle α (defined within the figure) equal to 0, π/24, π/12, π/8, and π/6. Figure 6 shows the numerically evaluated eigenvalues along these paths within the IBZ. The horizontal axis is the dimensionless quantity (aQ) 2 . The plots also correspond to the choice of λk z = 1.0. The dotted lines represent the longitudinal eigenvalue, Λ * L , while we use dashed lines for the transversal eigenvalue, Λ * T .
The different plots in Fig. 6a) clearly show that there is an angular dependence, which becomes more marked as we approach the edge of the BZ. This fact was already pointed out by Brandt in his early work on the elastic properties of vortex lattices 12 . The anisotropy is more pronounced in the dense regime, and diminishes as we increase a/λ. In the concentrated case with a/λ = 0.2, the relative differences |Λ(α) − Λ(π/12)|/Λ(π/12) are up to 20% for the longitudinal eigenvalue, and 60% for the transversal one. These ratios decrease in the intermediate region with a/λ = 1.0, and are of the order of 10% and 15%, respectively, for a/λ = 5.0, within the dilute regime. We also observe that the longitudinal eigenvalue is largest along the α = 0-direction, while the transversal eigenvalue attains its smallest value for α = 0.
Note that both eigenvalues, for all angles α, have the same value at Q = 0, determined by λk z and a/λ. (Since this is not immediately apparent from Fig. 6, the intersection point is marked with a black dot on the vertical axis.) Away from Q = 0, the transversal eigenvalue drops sharply, more so in the dense regime. Transversal modes are therefore less costly than longitudinal modes. Figure 6 also indicates that at a finite value of λk z (equal to 1.0 in this case), the minimum cost is obtained for a finite value of the rescaled wavevector aQ, which depends on the density of flux lines. The different lines are obtained for different paths along the IBZ, at angles α = 0, π/24, π/12, π/8, and π/6. The two solid lines in a) and b) depict eigenvalues from the nonlocal continuum limit in Sec. IV.
The solid lines in Figs. 6a) and 6b) are the results of the nonlocal continuum approximation introduced in subsection IV C, i.e., Λ L ≃ (c nl 11 Q 2 + c nl 44 k 2 z )/2 and Λ T ≃ (c 66 Q 2 + c nl 44 k 2 z )/2, with the appropriate parameters. From these figures, one can conclude that (at least for λk z = 1.0) the continuum approximation reproduces extremely well the behavior of both eigenvalues along the α = π/12 direction. This is certainly the case at the high-est concentration of flux lines (a/λ = 0.2), and even at the intermediate concentrations of a/λ = 1.0. However, the approximation fails at lower densities, and is not even included in Fig. 6c), as the predicted eigenvalues fall well outside the range of our plot.

B. Variations with kz
To explore the k z -dependence of the results, in Fig. 7 we plot the longitudinal and transversal eigenvalues as a function of (aQ) 2 for a fixed angle α = π/12, and for a number of values of a/λ and λk z . From these figures, we see that, not surprisingly, both eigenvalues increase as a function of λk z , and consequently, modes with higher k z are in general costly. For the smallest value of a/λ = 0.2 in Fig. 7a), the λk z dependence is most pronounced in a small region close to the BZ center. This region becomes larger as we increase either λk z or a/λ. We also note that both eigenvalues show a linear dependence on (aQ) 2 for a large fraction of the IBZ. The slope of this line is almost constant within each plot, i.e. for a given value of a/λ. The transversal eigenvalue for λk z = 0 is a straight line throughout the whole range of values of (aQ) 2 , as is the longitudinal eigenvalue for high enough values of λk z . These features are qualitatively well accounted for by the continuum limit results. According to Eqs. (40), (41), and (42), the transversal eigenvalue reduces toc 66 Q 2 /2 for λk z = 0, i.e., grows linearly with (aQ) 2 , with a slope nǫ o /8a 2 . In addition, for suficiently large values of λk z ≫ λ c Q, c nl 44 (Q, k z )k 2 z reaches a saturation value, rendering the termc 66 Q 2 /2 as the leading contribution depending on Q. Similarly, the longitudinal eigenvalue reaches a saturation value for λk z ≫ λQ, so that the term −nǫ o Q 2 /8 controls essentially the Q dependence of the eigenvalue at large values of λk z . Eqs. (43)-(44) fit very well the numerical resuls obtained for a/λ = 0.2 and λk z = 0.0, 1.0, 10.0 (we are not showing these results in Fig. 7a) for the sake of clearity of the figure). For the larger value of λk z = 50.0, although the qualitative form of the functions is quite well approximated by these equations, there are major quantitative differences with respect to the numerical results mainly due to a smaller value of both eigenvalues at Q = 0. In Fig. 7b) with a/λ = 1.0, there are important differences in the value at Q = 0 already for λk z = 10.0, as well as small deviations from the value of the slope of the linear part of the curves. These deficiencies of the continuum approximation are even more evident for a/λ = 5.0 in Fig. 7c).
In the latter case, we also note the emergence of singleline characteristics: The eigenvalues are clearly mostly a function of k z , with only small variations with (aQ) 2 .

C. Comparison to analytic approximations
We observe in Fig. 8 that Λ * L (Q = 0, k z ) = Λ * T (Q = 0, k z ) increase monotonously with λk z , with the singleline behavior of E(k z ) in Eq. (19) eventually taking over for values of λk z greater than a characteristic λk c z , proportional to (ǫa/λ) −1 . The contribution of the nonzero reciprocal lattice vectors in Eq. (23) becomes relevant for λk z > λk c z . In the nonlocal continuum approximation, i.e., considering only G = 0, Λ * L (Q = 0, k z ) = Λ * T (Q = 0, k z ) reaches a saturation value of 2πn 2 ǫ o at large values of λk z . Indeed, for high areal densities of a/λ = 0.2, our numerical outcome is very close to the nonlocal continuum approximation only up to λk c z . On the other hand, at low areal densities, such as for a/λ = 5.0, the numerical data are obviously much closer to the single-line limit, than to any of the other forms, over the whole range of values of λk z . As shown in Figs. 6c) and 7c), the variations of Λ L and Λ T with (aQ) 2 in the latter case are very weak. The eigenvalues throughout the IBZ essentially coincide with their values on the Q = 0-axis, which is set by the form of E(k z ). In the interval between these two values of a/λ, Λ L (Q = 0, k z ) = Λ T (Q = 0, k z ) gradually cross over from saturation type behavior to the parabolic-log dependence expressed in Eq. (19). The results of the local approximation naturally fit well the behavior of both eigenvalues for small values of Q and k z , as expected.
In discussing Fig. 7, we noted that while the non-local continuum expressions capture the qualitative form of the numerical results, there are also important quantitative differences. Using the insights gained from the numerics, we shall now present an analytic form that corrects some of these discrepancies. One difference from the continuum results arises from the limiting slope in Fig. 7: In the continuum limit, the termc 66 Q 2 /2 eventually determines the slope of the linear regime of the transversal eigenvalues. However, as shown in Fig. 4a), the actual shear modulus c 66 can be quite different from the limiting value ofc 66 used in the continuum approximation. The differences in slope thus merely reflect the differences between c 66 andc 66 as a function of a/λ, as already discussed in Sec. IV B. This deficiency of the continuum eigenvalues is thus removed by using the exact value of c 66 from Fig. 4a).
A second difference is an overall shift of the eigenvalues from the continuum limit prediction, which becomes more pronounced at larger λk z . This clearly originates in the differences appearing already in Figs. 8a) and b) for Λ L (Q = 0, k z ) = Λ T (Q = 0, k z ) at Q = 0 and sufficiently large k z . Since the continuum approximation uses only the term with G = 0 in Eq. (23), we introduce a correction by replacing the sum over all the remaining reciprocal lattice vectors with an integral. These corrections result in the following expression for the transversal eigenvalue: In this equation, the correct local shear modulus, which depends on the field strength or lattice spacing a through Eq.(30), is used in place of its continuum limit. Also the non-linear modulus c nl 44 is corrected by the additional integral over the nonzero reciprocal lattice vectors, properly normalized by the BZ area A BZ = 8π 2 /( √ 3a 2 ). To exclude the point at G = 0, the radial component of G in the above integral goes from C, defined by A BZ = πC 2 , to ∞ (or, if necessary, to the short distance cutoff ξ −1 ). After evaluating the integral, we obtain We have ascertained that this expression provides an excellent fit to the numerically obtained results for Λ T (Q, k z ) along the α = π/12 direction. In the next section we provide explicit comparisons of the analytic expression and numerics for different values of the anisotropy. The longitudinal eigenvalue along the α = π/12 direction can also be reasonably well fitted by a similar analytic expression, which replaces the term [c 66 Q 2 + c nl 44 k 2 z ]/2 in Eq. (46) with the expression in Eq. (44), with a further substitution of c 66 Q 2 /2 for nǫ o Q 2 /8 =c 66 Q 2 /2, as

D. Anisotropy
We conclude this section by discussing the dependence of our results on the anisotropy of the superconductor. The results presented so far correspond to a particular value of the anisotropy parameter, namely ǫ = 0.1, which falls within the range (from 1/10 − 1/5) reported for YBa 2 Cu 3 O 7 (YBCO) in the literature 1 . However, smaller values of ǫ, in the interval 1/100−1/50, characterize highly anisotropic materials such as Bi 2 Sr 2 CaCu 2 O 8 (BiSCCO). For such highly anisotropic materials, the discreteness of their layered structure becomes important, and one may well question the validity of a three dimensional Landau-Ginzburg description. An alternative model is a set of weakly (Josephson) coupled superconducting layers. The applicability of the three dimensional description is typically assessed by comparing the coherence length along the c axis ξ c = ǫξ, with the distance d between the Cu-O layers in the material 1 . A coherence length ξ c larger than the layer spacing usually justifies a continuous description along the c direction 33 . This is certainly the case for YBCO, but not necessarily for BiSCCO. Nevertheless, extrapolating the results of the continuous description may also provide insights into the elastic properties of highly anisotropic superconductors such as BiSCCO. To this end, we compare results obtained for three different values of the anisotropy parameter, namely ǫ = 0.02, 0.1, and 0.5, for the usual areal densities of a/λ = 0.2, 1.0, 5.0.
In Figs. 9 and 10, we again plot the elastic eigenvalues along the Q = 0-axis, and the transversal eigenvalue along the α = π/12 direction inside the IBZ, respectively, for different choices of ǫ. There is clearly a strong dependence on ǫ, mainly originating from the ǫ-dependence of the line tension E(k z ) in Eq. (19). Naturally, the most pronounced tendency is the decrease in eigenvalues with ǫ, reflecting the general softening of the vortex lattice. It is worth pointing out however, that the slopes of the linear part of the curves in Fig. 10 are independent of ǫ, reflecting merely the fact that the shear modulus c 66 is independent of the anisotropy parameter. The dashed lines in the three graphs of Fig. 9 correspond to the analytic expression of Eq. (46) (for the appropriate parameters a/λ and ǫ), and clearly provide an excellent fit to the numerical data along the Q = 0 axis. Similarly in

VI. THE MELTING LINE
The melting of a vortex lattice by thermal fluctuations has attracted considerable attention in the context of high temperature superconductors. This transition has been observed by means of several experimental techniques such as bulk magnetization, local induction, and latent heat measurements [21][22][23][24][25][26][27][28][29] . The line-like nature of the constituent elements provides intriguing challenges to theoretical analysis. Important features of the melting transition are the negative slope of the melting curve T m (B) at high fields, its reentrant behavior at low fields, and its marked dependence on anisotropy. Some of these features can be extracted from simple models of the vortex lattice, as in the so-called XY 32 , Bose 34 , and cage 35 models. These models agree in the prediction of certain universal features, such as the scaling of the melting temperature with the anisotropy parameter, or the magnetic field, in the high field region of the phase diagram 36 .
In the absence of a rigurous theory for threedimensional melting, the position and shape of the vortex lattice melting line is usually estimated using the socalled Lindeman criterion. According to this criterion, the lattice melts when thermally induced fluctuations of a lattice point become comparable to the lattice spacing. This condition can be written as where c L ∼ 0.1 − 0.2 is the (empirically chosen) Lindeman parameter. The extent of fluctuations is measured by the autocorrelation function, u 2 (r 0 ) . This quantity was calculated by Nelson and Seung 13 in the local limit, and by Brandt 15 and Houghton et al. 16 for the nonlocal continuum limit, and in more general circumstances including variations of the amplitude of the order parameter in the Ginzburg-Landau hamiltonian.
In the preceeding sections, we calculated the harmonic energy cost of fluctuations of the vortex lattice. In a clas-sical equilibrium state, each independent harmonic mode acquires a thermal energy of k B T /2. By adding the corresponding squared amplitudes of the normal modes, the autocorrelation function is obtained as The Lindeman criterion now provides a rough estimate of the melting temperature T m , through the relationship The leading contribution to the fluctuations comes from the transversal modes, which as discussed in Sec. IV are always smaller than the longitudinal modes. We also ignore the angular dependence observed in Figs. 6, which is rather weak, and should not give rise to qualitatively different results. Furthermore, we assume that the transversal eigenvalue is isotropic, and given by Eq. (46). This expression can be integrated analytically over the two-dimensional BZ, resulting in a rather long expression, that we have finally integrated numerically over λk z . The numerical integral over λk z is, in all cases, between a long wavelength cutoff equal to 10 −6 and the ǫ-dependent short scale cutoff of λ/ξ c = κ/ǫ 33 .
Provided that we are sufficiently far away from the critical temperature T c , one can further assume that the penetration length is independent of temperature, i.e. λ(T ) ∼ λ(0). As we get closer to T c , however, one must consider the temperature dependence of λ in order to obtain sensible results. For the sake of consistency of our analysis, within the Ginzburg-Landau critical regime, we then use the mean-field temperature dependences of both λ and ξ, i.e. λ(T ) = λ(0)(1 − T /T c ) −1/2 and ξ(T ) = ξ(0)(1 − T /T c ) −1/2 . We expect this approximation to fail in close vicinity of T c due to critical fluctuations. In Figs. 11, 12, and 13, we depict the resulting melting lines for different values of the anisotropy parameter ǫ = 0.02 (BiSCCO), and ǫ = 0.1 (YBCO). Note that some of the data are presented in a log-log or in a loglinear plot in order to emphasize the scaling with the magnetic field in the high field region, as well as to better visualize the lower part of the melting line, i.e. the reentrant portion of the phase diagram. The dotted line in the last two figures represents the upper critical field H MF c2 (T ) = φ o (1 − T /T c )/2πξ 2 (0), with T c = 93 • K and ξ(0) = 14Å. Close to H MF c2 (T ) the amplitude of the order parameter is strongly reduced due to the overlap of the vortex cores and these results are no longer valid. However, well below this line, the London limit provides a good description of the system. In Fig. 11, we represent the results obtained assuming λ(T ) = λ(0) = 1400Å, for ǫ = 0.02, 0.1, 0.5. At high fields a/λ ≪ 1, the melting temperature decreases as B −0.53 in all cases, consistent with the prediction of T m ∼ B −1/2 , common to the above mentioned models 32,34,35 . We have also checked that the melting temperature decreases linearly with ǫ, for various values of a/λ within the high field region, i.e. T m ∼ ǫB −0.53 . This is another feature which is in agreement with previous predictions, and with available experimental data. In this figure we also encounter some unreasonably high melting temperatures which are the consequence of the assumption λ(T ) = λ(0).  Figure 12 shows the results obtained for the melting line T m assuming λ(T ) = λ(0)(1 − T /T c ) −1/2 . We have used the values of T c = 90 for BiSCCO (ǫ = 0.02), and T c = 93 for YBCO (ǫ = 0.1). As in the previous figure, within the high field region of the phase diagram, the melting temperature decreases as the external field is increased. Quantitative and qualitative differences between the new curves and the corresponding results for λ(T ) = λ(0) (depicted as dashed lines in the figure) can be observed even at low temperatures. Near the critical temperature, the low-field portion of the melting curves has a negative slope (as for their high field counterparts), but at low temperatures the slope is positive, so that the melting temperature decreases with the external field at very low densities of flux lines.
At the scale of Fig. 12, the high and low field portions of the melting curve T m (B) appear almost discontinuous close to T c . However, as indicated in the blow-up of this region in Fig. 13, this is in reality a sharp but continuous turn around. Following the literature 26-28 , we fit the melting fields B m (T ) to power laws in the reduced temperature t ≡ 1 − T /T c , with an exponent of α. The dashed lines in Figure 13 indicate the extent of this power law regime for each value of ǫ. We note that even though we are using the input value of T c , the power law regime extends at most over two decades. We naturally obtain different exponents for the upper and lower branches of the melting curve. The apparent exponent of the upper branch actually depends on the parameter ǫ, taking the value of α + (ǫ = 0.1) = 1.91 for 'YBCO', and a slightly smaller value of α + (ǫ = 0.02) = 1.67 for 'BiSCCO'. Both curves overlap in the reentrant region at very dilute con-centrations where we can fit to an effective exponent of α − = 0.75.
The phase diagram of Fig. 12 is qualitatively similar to the predictions of Ref. 9, where the power law forms for the phase boundary close to T c were also first discussed. In particular, it is reasonably straightforward to estimate the shape of the phase boundary in the low field region, where the leading contributions to the transversal modes Λ T come from the shear modulus c 66 (which decays exponentially with a/λ, and does not depend on ǫ), and from the last term in Eq. (46). As expected for very low densities of flux lines, the last two terms in this equation tend to the single line tension E(k z ) (see Eq. (19)) which, for small values of k z , is in turn dominated by the ǫ-independent magnetic contribution. This leads to B − m ∼ λ(T ) −2 ln −2 (λ(T )), and consequently to an exponent of α − = 1 with logarithmic corrections 9 . The effective exponent of α − ≈ 0.75 obtained in Fig. 13 is different from the α − = 1 expected on the basis of our mean-field form for λ(T ) indicating the difficulty of determining the true exponent from such fits. Close to the critical temperature Tc, the melting induction grows as a power of the reduced temperature 1 − T /Tc, as indicated by the asymptotic dashed lines. As in Fig. 11, the dotted-dashed lines show the melting lines resulting from the local elastic moduli.
For the upper portion of the mean-field line, the melting field is estimated 9 to scale as B m ∼ λ(T ) −4 , leading to α + = 2 for our choice of a mean-field λ(T ). While the observed value of α + (ǫ = 0.1) = 1.91 is not far from this expectation, the smaller value of α + (ǫ = 0.02) = 1.67 indicates the difficulty of finding a good asymptotic regime. Naturally, in close vicinity of the critical point we can no longer use the mean field forms for the divergence of λ ∝ ξ ∝ |t| −1/2 . Using the scaling forms, ξ ∝ |t| −ν ∝ λ 2 , with the XY critical exponent of ν ≈ 2/3, leads to α + = 2α − = 2ν 9 . Indeed, the values of the exponent of the upper branch (α + ∼ 1.33 − 1.36) reported from measurements of the melting transition in single crystals of YBCO [26][27][28] , are consistent with this prediction. However, given the difficulties of extracting exactly known exponents from the data of Fig. 13, we may well question the robustness of this procedure.
We would like to emphasize that most of these results cannot be obtained from the local values of the elastic moduli given in Sec. IV B. The results obtained from such a local elastic approximation are plotted as dotdashed lines in Figs. 11 and 13. The latter curves overlap in the high field region, where neither the local tilt, nor the shear modulus depend on ǫ. Rather than decreasing as T m ∼ B −1/2 , the resulting melting temperature reaches a constant value for high B, i.e. T m ∼ B 0 ǫ 0 . On the other hand, the local melting temperature provides a very good approximation to the melting line at very dilute concentrations, in the reentrant low field portion of the curve.

VII. CONCLUSIONS
In conclusion, we have obtained the elastic moduli of the flux line lattice in a systematic and detailed way which is not restricted to the most commonly considered continuum limit. The transversal and longitudinal harmonic eigenvalues have been computed numerically for different areal densities of flux lines, and as a function of Q (within the irreducible Brillouin zone) and k z . Several features emerge from the analysis of our results: (i) There is a weak angular dependence of the harmonic eigenvalues which becomes more pronounced on approaching the BZ boundary. (ii) Throughout the BZ, transversal modes are less costly than longitudinal modes, and are the main cause of lattice fluctuations. (iii) Not surprisingly, both eigenvalues increase with λk z , and modes with higher k z are more costly. (iv) Rather surprisingly, due to a rapid decrease of c nl 44 , the energy of a transversal mode with a nonzero k z actually goes down with Q. The minimum cost occurs for a finite Q which depends on λk z , and the density of flux lines. (iv) For a large portion of the IBZ, both eigenvalues exhibit a linear dependence on (aQ) 2 , whose slope depends on the local shear modulus c 66 which can be calculated as a function of a/λ. Some of the above features are qualitatively well accounted for by the nonlocal continuum limit results recalculated in Sec. IV. Nevertheless, beyond a characteristic value of λk z , there are major differences between this analytic form and those calculated numerically. Guided by our results, we propose analytic corrections to the nonlocal continuum limit, which fit quite well the behavior of the elastic eigenvalues throughout the London regime. Other limiting forms for the elastic energy often quoted in the literature, such as the single line, local, and local continuum approximations, are reconsidered within our framework, and their range of validity is examined. We have plotted some of these limits, together with our numerical results, for better ease of comparison.
We also consider different values of the anisotropy parameter ǫ, in particular corresponding to values quoted in the literature for high temperature superconductors such as YBCO or BiSCCO. The latter is the paradigm of a highly anisotropic superconductor with softer elastic eigenvalues. This fact can be corroborated in our analysis and is ultimately responsible for strengthening the magnitude of thermal fluctuations. We have made use of our proposed analytical expression to calculate the extent of thermal fluctuations for different ǫ. The full form of the melting curve is then obtained using the Lindeman criterion, capturing the following salient features of T m (B): (i) It decreases at high temperatures with the magnetic field, approximately as B −0.53 . (ii) It decreases with the anisotropy parameter, and consequently the molten vortex liquid covers a large fraction of the equilibrium phase diagram at small ǫ. (iii) There is reentrant melting at low fields due to the weak interactions between widely separated flux lines. The reentrant phase boundary is itself a non-monotonic function of temperature. (iv) Close to the critical temperature T c , both branches of the melting line B m can be fitted to power laws in the reduced temperature 1 − T /T c . However, the power law is observed at most over two decades, and the resulting effective exponents are different from the expectations from mean field theory, casting doubt upon the effectiveness of this method [26][27][28] .
Another interesting potential extension of this work is to include entropic contributions to the free energy, in order to explore the fluctuation induced effects which have been proposed to be important in the low field region of the phase diagram 37 . Our analysis has been applied to a conventional s-wave superconductor with a triangular lattice of flux lines. A square lattice of flux lines is also possible in d-wave superconductors, where similar approach and phenomenology can be carried out, although the details of the harmonic energy, as dictated by symmetry, will be different. Again, calculation of the entropic contributions to the free energy may provide a better estimate of the transition between triangular and square lattices 38 .

ACKNOWLEDGMENTS
We are grateful to R. Pastor-Satorras for critical reading of the manuscript. This research was supported by grants from the Direcció General de Recerca (Generalitat de Catalunya), and the National Science Foundation (Grant No. DMR-98-05833).