{'fraction_non_alphanumeric': 0.04038885288399222, 'fraction_numerical': 0.016176279974076474, 'mean_word_length': 4.427124366910523, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Density functional theory of vortex lattice melting in layered superconductors: a mean-field--substrate approach

We study the melting of the pancake vortex lattice in a layered superconductor in the limit of vanishing Josephson coupling. Our approach combines the methodology of a recently proposed mean-field substrate model for such systems with the classical density functional theory of freezing. We derive a free-energy functional in terms of a scalar order-parameter profile and use it to derive a simple formula describing the temperature dependence of the melting field. Our theoretical predictions are in good agreement with simulation data. The theoretical framework proposed is thermodynamically consistent and thus capable of describing the negative magnetization jump obtained in experiments. Such consistency is demonstrated by showing the equivalence of our expression for the density discontinuity at the transition with the corresponding Clausius-Clapeyron relation.


I. INTRODUCTION
The melting of the vortex lattice is one of the most striking aspects of the phenomenology of hightemperature superconductors 1,2,3 .
It is now both theoretically 4,5,6,7,8 and experimentally 9,10 established that enhanced thermal fluctuations render the Abrikosov flux-line lattice unstable to a vortex liquid over a large part of the B-T phase diagram in these materials. Much interest has been devoted to the melting transition in extremely anisotropic materials 7,11,12,13 such as Bi 2 Sr 2 CaCu 2 O 8 . The large value of the anisotropy 14 in these materials motivates their modeling in terms of a one dimensional array of magnetically-coupled twodimensional superconducting layers with no inter-layer Josephson coupling. Vortex lines are then described by linear arrangements (stacks) of pancake vortices 15 , whose highly anisotropic interactions are mediated only by the magnetic field.
The first theories for the melting of the vortex lattice 5 were based on the Lindemann criterion 16 . The Lindemann melting criterion yields qualitatively accurate estimates for the location of the first-order melting transition. However, since it is based on the crystal properties only, it describes an instability rather than the melting itself. More detailed approaches, in particular analyses encompassing both solid and liquid phases, are thus required for a quantitative description of the transition and of the discontinuities in thermodynamic quantities upon melting.
Apart from numerical simulations 17 , the extreme limit of zero Josephson coupling between the layers has been analysed through classical Density Functional Theories 11,18,19 (DFT) and using a "mean-field" substrate approach 13 . Even if only the electromagnetic interaction is retained, the analysis of the melting transition is a challenging task due to the long-range character of the interaction. At large magnetic fields the strong in-plane vortex repulsion dominates over the out-of-plane interactions and the melting line approaches the 2D melting temperature T 2D m of the individual planes. On the other hand, at low fields, ignoring the possibility of low-field reentrance 4 , only few pancake vortex stacks are present in the system. Thermal fluctuations trigger the evaporation of isolated vortex stacks 15 at T BKT , in correspondence to the zero-field two-dimensional Berezinskii-Kosterlitz-Thouless (BKT) transition. For intermediate values of the external magnetic field the melting line interpolates between the BKT transition temperature T BKT T c close to the material critical temperature and the twodimensional lattice melting temperature T 2D m ≪ T c . In this regime, the full three-dimensional character of the system must be accounted for. The weak long-range outof-plane interaction motivates a mean-field treatment 13 , in which the 2D pancake vortex system in each layer is considered separately, with the effects of the other layers included via a self-consistent substrate potential. Extensive Monte Carlo and molecular dynamics simulations have confirmed the validity of this approximation and the accuracy of its results 17 .
If the results of previous DFT approaches to the numerical melting line are compared with the simulation data, consistent results are obtained at magnetic fields larger than about 0.5B λ (B λ = Φ 0 /λ 2 , where Φ 0 = hc/2e is the unit flux and λ the planar penetration depth). However, at lower fields, the disagreement between the two theories is substantial, with the DFT melting line shifted to much higher temperatures in comparison with the simulation data. Another difficulty with previous DFT studies is associated with the prediction of the sign of the magnetization jump across the transition. As it is well known, the vortex lattice exhibits a negative jump in magnetization upon freezing, leading to a solid phase which is less dense than the liquid, as in the ice-water transition. Such anomalous behaviour was first obtained in direct measurements of the magnetization discontinuity across the transition 9 . These observations are consistent with the negative slope of the melting line in the T -H phase diagram taken together with the Clausius-Clapeyron relation. However, earlier DFT analyses of the freezing transition in the pancake vortex system obtained thermodynamically inconsistent results 18 , reporting a positive density change and a negative slope of the melting line.
In this paper, we adapt the DFT analysis to incorporate ideas from the substrate approach; our analysis is simpler and more accurate in comparison with simulation data. Moreover, it addresses and solves the problems with the thermodynamic consistency faced in previous studies. The approach of this paper was first described in a recent letter where the effects of a surface on the melting transition have been studied 20 . In the present work we concentrate on the implications of our approach for the bulk transition in the infinite system.
The main methodological difference between the work described here and previous DFT approaches lies in the determination of the direct correlation function. While in previous work, this correlation function was derived ab initio from the microscopic vortex interaction using the hypernetted chain approach or more elaborate extensions, here it is obtained by combining results from Monte Carlo simulations of 2D logarithmically interacting particles, i.e., the One Component Plasma (OCP), with the substrate potential approach. The correlations of the 2D OCP are used to describe the effects of the strong inplane vortex interaction, while the substrate potential accounts for the out-of-plane contributions. Within this new approach we obtain a simple expression for the freeenergy which can be extended to the inhomogeneous case.
We also prove the thermodynamic consistency of our approach by showing how to obtain the negative density jump across the transition, fully consistent with the Clausius-Clapeyron equation. We do this by including a constraint which enforces an integer number of particles per unit cell 21 . Contrary to earlier claims 18 , no higher order correlation functions (three point or more) are needed to obtain the anomalous sign of the density discontinuity.

II. MODEL
Strongly anisotropic layered superconductors, such as the Bi-based compounds, are conveniently described in terms of the Lorentz-Doniach model 23 . The basic topological objects are two dimensional vortices (pancake vortices) with a core limited to a single superconducting layer. The interaction between pancake vortices is strongly anisotropic due to the underlying layered structure. In Fourier space the potential reads 1 here d is the layer spacing and λ the bulk penetration depth in the plane. When placed on the same layer, pancake vortices feel a strong repulsive interaction, logarithmically dependent on their separation, where ε 0 = (Φ 0 /4πλ) 2 is the vortex line energy and ξ the (in-plane) correlation length. The interaction between vortices residing on different layers is attractive. Fourier transforming (1) back along the z coordinate we obtain with K + = 1/λ 2 + K 2 . Finally, transforming to planar real coordinates as well one obtains the potential in real space For large in-plane separations R ≫ λ a logarithmic attractive interaction is obtained V z =0 (R) ∼ −(d/λ)e −z/λ V 0 (R), suppressed by a factor d/λ when compared with the in-plane one. This out-of-plane interaction decays exponentially in the z direction over a distance λ, extending over a large number (λ/d) of layers. The strong anisotropy in the vortex interaction allows us to separate the strong in-plane repulsion from the weak out-of-plane interaction. The overall effect of the latter can then be accounted for via an effective substrate term 13 .

III. CLASSICAL DENSITY FUNCTIONAL THEORY
We follow the classical Density Functional Theory of freezing of Ramakrishnan and Yussouf described in Refs. 24,25 by choosing the uniform liquid as the reference state and estimating the difference in free energy due to the appearance of finite density modulations. For simplicity, we consider here a generic three-dimensional system of interacting particles; the modifications needed to describe the strongly anisotropic pancake vortex system are presented in the next section.
The spatial arrangement of particles is described through the density field where r i is the position of the i-th particle (the index µ emphasizes that ρ µ (r) describes the non-averaged microscopic density). To analyze the finite temperature behavior of the system we consider the density ρ(r), averaged over thermal fluctuations where the brackets . . . denote the thermal average. The liquid and solid phases are characterized by qualitatively different density fields ρ(r): in the liquid phase the particles are delocalized across the system and the averaged density ρ(r) =ρ 3D is constant; on the other hand the solid phase is characterized by a modulated density ρ(r), with peaks at the lattice points. The appearance of finite density modulations is a consequence of particle-particle correlations arising from microscopic interactions. The classical DFT is based on the assumption that the free energy can be written as a functional of the averaged density ρ(r). One starts from the ideal gas free energy describing a non-interacting liquid and includes the correlations via an effective quadratic term in the density modulations δρ(r) = ρ(r) −ρ 3D . Within this approximation the grand canonical free energy difference relative to the uniform liquid reads where the temperature T is measured in unit of energy with k B = 1. The first two terms generalize 26 the standard free energy of an ideal gas to the case of nonhomogeneous systems. The double integral term incorporates the effects of interactions up to second order in the density difference δρ. This is the term which is responsible for the appearance of finite density modulations which are absent in a non-interacting system. Therefore, the key input in this theory is the function c(r), the socalled direct pair correlation function, which accounts for correlations in the reference liquid. The direct pair correlation function can be related to more transparent physical quantities such as the static structure factor 25,27 . Here, we give a brief derivation of the main relations which will be needed in our following discussion and, in particular, in the derivation of the Clausius-Clapeyron equation. We start from the microscopic density-density correlator The Fourier transform of the density-density correlator defines the structure factor 28 Next, we calculate the structure factor from the free energy (7). The second functional derivative of the free energy δΩ with respect to ρ(r) evaluated at ρ(r) =ρ 3D is the (functional) inverse of the density-density correlator (see also Ref. 25) The functional inverse can be easily calculated in Fourier space. In this way, we find a relation between the structure factor and the direct correlation function 29 Therefore, apart from additive constants (or delta functions in real space), the direct correlation function c(q) is given by the inverse of the structure factor S(q). The q = 0 mode of the structure factor is related through the fluctuation-dissipation theorem to the isothermal compressibility of the system, i.e., S(q = 0) = ( N 2 − N 2 )/ N 2 =ρ 3D T κ T (the isothermal compressibility of the ideal gas is κ 0 The study of the free energy functional (7) requires knowledge of the direct correlation function c(r), a quantity which is usually obtained from the liquid state theory. It is possible to write a diagrammatic expansion for c(r) in terms of the microscopic two-body potential V (r) (or more precisely of the Mayer function f (r) = exp(−V (r)/T ) − 1). For weak potentials, high-order correlations can be neglected and c(r) is given by the first (unperturbed) term 27 In order to obtain a better estimate for c(r), higher order terms in the perturbation expansion must be included. This is usually done by selecting specific classes of diagrams out of the complete perturbative series. The approach that is most widely used is the hypernetted chain (HNC) closure, an approximation scheme which, however, is known to underestimate liquid-state correlations. The strategy we pursue here is different: we exploit the specific properties of the pancake vortex lattice by considering the system as a collection of two-dimensional systems of log-interacting particles subject to a periodic modulated substrate potential due to the other layers, as done in Ref. 13. The function c(r) then combines an in-plane correlator, arising from the strong in-plane logarithmic interaction, and the weak out-of-plane potential.

IV. DFT-SUBSTRATE APPROACH
In contrast to standard liquids, the pancake vortex system exhibits a strong uniaxial anisotropy. Hence, we consider separately the in-plane (R) and out-of-plane (z) dependencies: in particular ρ(r) → ρ z (R) becomes a sequence (in z) of two-dimensional densities. Similarly, we define the density variations δρ z (R) = ρ z (R) −ρ and the direct pair correlation function c z (R) (note thatρ is a 2D density). The DFT free energy of Eq. (7) can be adapted to the anisotropic vortex liquid The only input needed in the DFT free energy is the direct correlation function c z (R) which we obtain by implementing the substrate model for the pair correlation function, where V z (R) is the out-of-plane interaction of Eq. (4). Within the planes, vortices are strongly correlated due to the repulsive logarithmic interactions (Eq. (2)). Hence, we can approximate c 0 (R) with the direct correlation function c 2D (R) of the two dimensional logarithmically interacting particles (also known as one-component plasma, OCP). The in-plane component also contains contributions from the out-of-plane interactions V z =0 (R), which are small however, since V z =0 (R) appears in the perturbative expansion 27 of c 0 (R) at least to quadratic order, Hence, the overall contribution due to out-of-plane interactions of all planes adds We use results of Monte Carlo simulations of the twodimensional OCP at various coupling constants Γ = 2ε 0 d/T to extract c 2D (R). We have performed simulations on a system of 256 particles, using an alternative to the traditional Ewald summation method proposed recently by Tyagi 30 . Thermodynamic data and correlations were averaged over ∼ 3 × 10 3 independent measurements following equilibration. The program was benchmarked using available numerical results for correlations and thermodynamic functions. In the simulations, the structure factor S(K) was calculated to yield the direct correlation function via the relation c(K) = 1 − 1/S(K), see Eq. (10).
In the determination of the out-of-plane direct correlation function c z =0 (R) we neglect the higher order terms in the potential expansion, approximating it with the leading unperturbed value −V z =0 (R)/T (cf. Eq. (12)). Higher orders in c z =0 (R) involve at least terms of order ∼ (d/λ) 2 . In the following, we will see that the relevant quantity in our analysis is the total out-of-plane correlator, defined as (dz/d)c z =0 (R). Whereas the leading term is of order (λ/d)(d/λ) ∼ 1 and hence comparable to the in-plane component, the subleading term is of order (λ/d)(d/λ) 2 ∼ (d/λ). We neglect this contribution to be consistent with our approximation for the in-plane component of the correlator. At a mean-field level the thermodynamically stable state corresponds to the minimal free energy configuration of the functional (13). Then, the density functions ρ z (R) must obey the saddle point equation A key quantity in our discussion is the molecular field 24,31,32 ξ z (R) defined through At the minimum of the free energy, combining the saddle point equation (15) with (16), the molecular field be- . (17) Hence, the molecular field represents the screening potential, produced by the modulated density which is acting back on the density itself. However, while (16) defines the molecular field and hence applies always, the interpretation in terms of a screening potential following from Eq. (17) is valid only at the minimum.

V. FREE ENERGY IN FOURIER SPACE
In thermodynamic equilibrium all superconducting planes are equivalent and the averaged vortex density ρ z (R) becomes independent of the layer position z; we write ρ z (R) = ρ(R). Next, instead of seeking the exact form ρ(R) solving the non-linear integral equations (15), we restrict our analysis to a simple family of periodic functions which model the modulations of the density in the triangular crystalized phase. In the following, we concentrate on the simplest case, retaining only the first Fourier components of the density in a triangular lattice where the vectors K 1 are the first reciprocal lattice vectors of the frozen structure and depend on the area a of the unit cell, µ = δρ(K 1 )/ρ is the Fourier component of the density with wave length K 1 , and η = δρ(K = 0)/ρ is the relative density change upon freezing. Within this approach the density is characterized by the three variables µ, η, and K 1 . Note that, however, in a scheme where the number of particles in a unit cell is constrained to be an integer, the size of the unit cell a in the crystalized structure and the value of the density jump are related and, thus, K 1 and η are not independent variables; we return to this issue in the next section. The function g K1 (R) = K1 e iK1·R = 2 cos(2x) + 4 cos(x) cos(ỹ) (19) includes the sum over the six first reciprocal vectors in the triangular lattice; in the last equality we have defined the dimensionless variablesx = xK 1 /2 andỹ = √ 3 yK 1 /2. We also write a similar Ansatz for the molecular field retaining only the zeroth and first Fourier components, ζ and ξ respectively, consistent with (17) and the rapid decay of the correlator c(K) in the liquid phase, cf.  25)) for the full three-dimensional pancake vortex system at the melting field Bm (Bm/B λ ≈ 0.099 for T /ε0d = 0.1). We define our wave-length unit as G ≡ (8π 2ρ / √ 3) 1/2 . The substrate potentialṼstack(K) modifies the full correlator c 2D sub (K) in two different ways as compared to the in-plane correlation function c 2D (K): i) it enhances the correlations of the liquid, pushing the melting to temperatures larger than T 2D m and ii) it shifts the peak of the correlation function to a value of K which is smaller than G. In the incompressible limit K1 = G and at melting the full three-dimensional correlator assumes the critical valuē c 2D sub ≡ c 2D sub (G) =cc ≈ 0.856, cf. (42). Considering a finite compressibility, the system is allowed to gain correlation energy by crystalizing at K1 < G, hence, at a density smaller than that of the liquid. However, a large density change is prevented by the finite compressibility of the system and the crystalized structure is characterized by a first reciprocal lattice vector K1 G (see text).
help of the relations we project the zeroth and first Fourier components of ρ(R) =ρ exp(ζ + ξg K1 (R)) and obtain where we have defined the function This function is unaffected by rescaling of the unit cell area a (equivalently, it does not depend on K 1 ). Substituting the Ansätze (18) and (20) in the DFT free energy (13), we obtain the two-dimensional free energy density δω 2D sub ≡ (d/ρV ) δΩ as a function of the order parameters η and µ and the length K 1 of the first reciprocal lattice vectors, where ξ has to be understood as a function of η and µ through (22). In (24) we have defined the correlator whereṼ z (K) =ρV z (K) is the dimensionless Fourier transform (with an additionalρ factor) of the out-ofplane pancake vortex interaction (Eq. (3)). The out-ofplane interactions contribute to the correlator through the total stack potential which enhances the nearest-neighbor correlations, cf. Fig. 1.
Within this approach only the K = 0 and the K = K 1 components of the correlator are present in the expression for the free energy (24). The K = K 1 component measures the correlations which are responsible for the solidification of the liquid; its dependence on temperature and field yields the position of the melting line and its local dependence on K determines the discontinuities in the first order transition. On the other hand, the K = 0 component c 2D sub (0) of the correlator is related to the compressibility of the vortex system via (11). Using the relation κ T = 1/c 11 (0) between the compressibility and the compression modulus and using the expression 33 c 11 (0) = B 2 /4π, we find that consistent with the result of the present substrate-based approach: here, the K → 0 divergence in the correlator c 2D (K → 0) ∼ −4πρε 0 d/T K 2 of the incompressible Coulomb gas 34 is cancelled by the corresponding divergence in the out-of-plane component of the correlator, see (26), and the remaining term reproduces (27). In the Coulomb gas problem, this cancellation has to be achieved by the introduction of a compensating background.
Our expression for the free energy (24) in terms of the density ρ has to be compared with the original formula in Ref. 24, where the free energy has been given in terms of the molecular field ξ. The two approaches are related via a Legendre transformation which adds the energy of an external periodic potential u(r), δW (u) = min ρ [δΩ(ρ) − d 3 r u(r)ρ(r)]. The relation (17) then includes the external potential u, ξ z (R) = u z (R) + (1/d) dz ′ d 2 R ′ c z−z ′ (|R − R ′ |)δρ z ′ (R ′ ). Setting u = 0, both approaches provide identical results. However, the formulation given here is more appropriate when describing non-uniform configurations as they occur in a solid-liquid interface or near a surface, cf. Ref. 20.

VI. CONSTRAINED THEORY
Within our approximation for ρ(R), the state of the system is characterized by three parameters: the density modulation µ at the first reciprocal lattice vector, the density change η across the transition, and the wave number K 1 . However, a theory based on the Ansatz (18) and the functional (13), in which η, µ, and K 1 are independent variables, is not fully consistent. The problem, as was pointed out in Ref. 21, is that in (18) the density jump η and the wave number K 1 appear as independent variables. This theory therefore may lead to the appearance of states with a non-integer occupancy per unit cell and hence to an incorrect description of the solidification of the liquid 22 . To solve this inconsistency only states with a fixed integer total number of particles per unit cell should be considered. One must then proceed with a constrained minimization of the free energy, introducing a Lagrange multiplier χ which enforces the so-called 'perfect crystal' condition 21 i.e., each unit cell contains exactly one vortex. With this additional term the free energy density becomes where δω 2D sub (η, µ, K 1 ) is given by (24). Substituting the Fourier Ansatz in the Lagrange multiplier term, we obtain the expression for the constrained free energy where G = (8π 2ρ / √ 3) 1/2 is the length of the first reciprocal lattice vector associated with a solid with the same density as the liquid.
To obtain the constrained minimum of the free energy, one needs to solve the saddle point equations for the variables χ, µ, η, and K 1 . Taking the derivative with respect to the Lagrange multiplier χ, we recover the 'perfect crystal' constraint (28) in the form which yields a relation between η and K 1 . The case K 1 = G describes an incompressible system, since the solid and the liquid have the same density. However, an ordinary first order phase transition is characterized by a finite jump of the density, and hence by a non-zero η. Consistent with our constrained theory, a finite η corresponds to the crystallization into a solid with a vortex density n v =ρ and, hence, with a first reciprocal lattice vector K 1 ≡ (8π 2 n v / √ 3) 1/2 which is different from G, K 1 = G. When K 1 > G the solid is denser than the liquid, as in conventional materials, whereas K 1 < G leads to an anomalous density jump with a solid which is less dense than the liquid, as is the case in the water-ice transition.
Next, the extremum condition of the free energy (30) with respect to µ gives the relation which is essentially Eq. (17) in Fourier space. The equation for η is modified as compared with the standard unconstrained theory by the presence of the Lagrange multiplier, which can be found from the minimization of (30) as a function of K 1 Combining (32)- (34) and the relation between µ and ξ as given by (22), we can eliminate ξ and χ. The saddle point equations for η and µ can be written as , (35) where η and K 1 are related by (31). The homogeneous (µ = 0) and uncompressed (η = 0) liquid always solves these equations. However, at low temperatures, other non-uniform solutions (µ = 0) may appear. If their corresponding free energy is smaller than that of the liquid, δω 2D sub < 0, the system freezes into the periodic (crystal) structure (note that at the minimum the perfect crystal constrained is fulfilled and δω 2D sub = δω 2D sub ). The equations (35) and (36) can also be obtained directly from the free energy (24), where K 1 is written as a function of η via (31). The minimization of δω 2D sub (η, µ)/T with respect to η and µ then yields the saddle point equations (35)

VII. INCOMPRESSIBLE LIMIT AND MELTING LINE
The pancake vortex system is essentially incompressible (|c 2D sub (0)| ≫ 1) in a wide portion of the phase diagram and the corresponding value of the density jump is small, i.e., η ≪ 1. Hence, the determination of the melting line can be done within the incompressible limit, with η = 0 and K 1 = G. This approximation is sufficiently accurate also for small magnetic fields (see later for estimates of the value of η). In the incompressible limit, the free energy becomes a function of µ alone, wherec 2D sub ≡ c 2D sub (G) and ξ(µ) is implicitly defined by The effective three-dimensional correlatorc 2D sub entering (38) is given by the sum of two contributions: the OCP correlation functionc 2D ≡ c 2D (G) and the stack potential V stack (G)/T of Eq. (26), both evaluated at K = G. While the first depends on temperature only, the latter depends also on the vortex density and thus on the magnetic field, where we useρ/G 2 = √ 3/(8π 2 ) and At large values of B, the inter-plane interaction is negligible and the full correlator reduces to the 2D-OCP componentc 2D . The temperature enters via the T -dependence of the direct correlation functionc 2D (T ), changing the coefficient of the quadratic term (as in the φ 4 Landau theory). We obtainc 2D (T ) directly from MC simulations; the results are shown in the inset of Fig.  3. Increasing the temperature, the liquid correlations weaken, S(G) decreases, and so doesc 2D , cf. Eq. (10).
As a function of µ, the free energy exhibits the shape of a Landau theory describing a first-order phase transition. In Fig. 2 we plot the free energy as a function of µ for different values ofc 2D . At large temperature, the correlatorc 2D is small and δω 2D (µ) exhibits only one minimum at µ = 0 with a value δω 2D (0)/T = 0, in correspondence with the (homogeneous) liquid phase. Decreasing the temperature (which corresponds to increasingc 2D ), a second local minimum µ s (metastable solid) with energy appears in addition to the liquid minimum at µ l = 0. Freezing occurs when the liquid and solid minima assume the same value of the free energy, i.e., when δω 2D (µ s ) = 0. Within our single order parameter theory, this condition is equivalent to a simple equation for the correlator 31 Going to even lower temperatures,c 2D further increases, the solid minimum decreases in value, δω 2D (µ s )/T < 0, and the crystal becomes the only thermodynamically stable phase. Monte Carlo simulations 34 show that the 2D-OCP freezes at T 2D m ≈ ε 0 d/70 where, however, the correlator assumes the valuec 2D ≈ 0.77 <c c (Γ 2D m = 2ε 0 d/T 2D m = 140). This disagreement is due to the approximations we have adopted in our analysis. In particular, at low temperatures, the higher-order peaks become important and more terms in the Fourier expansion have to be retained 26 .
At lower magnetic fields, the inter-plane correlation becomes important and the 2D correlationsc 2D are augmented by the stack potential V stack (G). The critical conditionc 2D sub (T, B) =c c can be solved together with (40) and yields a simple expression for the melting line B m (T ) This melting line is plotted in Fig. 3 (lower) together with the numerical results of the MC/MD simulations 17 . Furthermore, we find improved agreement in comparison with previous DFT studies where the direct correlation function was derived ab initio through approximative closure schemes such as the hypernetted chain or the more elaborate Rogers-Young approach 11,19 . In particular our novel approach approximates well the numerical results for small values of the magnetic fields well (B 0.5 B λ ), in a regime where results from previous work exhibit a substantial disagreement.

VIII. DENSITY JUMP AND CLAUSIUS-CLAPEYRON RELATION
To quantify the density jump across the transition, we need to consider the saddle point equation (36) for η where we have used K 1 = G √ 1 + η ≈ G(1+η/2) to linear order in η from (31). Previous analyses 18 of the freezing of the pancake vortex system were based on the unconstrained free energy (24). In the unconstrained theory the first reciprocal lattice vector is fixed at K 1 = G and the equation for η contains only the first term in (44), since the second term is due to the Lagrange multiplier, cf. (34). Hence, within the unconstrained theory, one obtains ('nc' stands for not-constrained) where we used (27) for c 2D sub (0). Given that Φ(c 2D sub µ s ) > 0, the sign of η is always positive. This result contradicts the Clausius-Clapeyron relation and experimental evidence that vortices, like water, freeze into a solid which is less dense than the liquid. In a magnetic system the Clausius-Clapeyron relation reads 3 where ∆B = B l − B s = −ηB l = −ηΦ 0ρ is the jump in magnetic induction and ∆s = ∆S/V = (S l − S s )/V the (positive) jump in entropy density on heating. Ignoring the small difference between H and B, equation (46) can be rewritten as ('CC' stands for Clausius-Clapeyron) Combining (47) with the negative slope of the melting line B m (T ) of Eq. (43), we obtain that the density jump on heating is positive, which corresponds to a negative η, thus η CC < 0. Therefore, a theory with a positive η and a melting line with a negative slope is not thermodynamically consistent. The second term in (44) resolves this inconsistency. In 2D systems the first maximum of the direct correlation function shows up at K max ≈ G, and hence ∂ K c 2D sub (K)| G (and the second term in (44)) is zero and η > 0. However, for the 3D pancake vortex system the substrate potential shifts the first peak of c 2D sub (K) to a value of K which is smaller than G, K max < G (cf. Fig. 1). The system gains correlation energy by crystalizing at a K 1 < G, hence, at a density which is smaller than that of the liquid. However, a large density change is prevented by the finite compressibility of the system and the crystalized structure is characterized by a first reciprocal lattice vector K 1 below but still close to G, Fig. 1. The second term of (44) is therefore negative and a negative solution for η becomes possible.
Next, we confirm the thermodynamic consistency of the constrained theory by showing that (44) and the Clausius-Clapeyron equation (47) are equivalent. To compare (47) with (44), we need to calculate the jump ∆s in entropy density in (47). The latter is given by the temperature derivative of the free energy difference where we have used that δω 2D sub = 0 along the melting line. The second term is of order η 2 ≪ 1 and can be neglected when compared to the first one. In order to calculate the entropy jump and η CC , we need to evaluate the partial derivative ∂c 2D sub (K 1 )/∂T at melting. The standard way 18,31 is to estimate ∂c 2D sub (K 1 )/∂T from the temperature dependence of the solid structure factor. Here we proceed differently. Comparing (44) with (47) (combined together with (48)), we see that in (44) the partial derivative of c 2D sub with respect to K 1 appears whereas (47) contains the partial derivative with respect to T , once we use (48) for the entropy jump. To compare the two different expressions for η we need to find a way to connect these two partial derivatives. This relation can be found from the critical condition which determines the melting line as we show in the following.
The system freezes when the free energy at the solid minimum vanishes. Substituting the values of µ s and of η s at the minimum into the expression for the free energy (24) (or, equivalently, into (30)), the freezing condition reads where K 1 is related to η s through (31) and we have neglected in (24) the term (1−c 2D sub (0))η 2 s , which is quadratic in η s . At the minimum, the molecular field ξ s and the order parameter µ s are related through ξ s = c 2D sub (K 1 )µ s , so the freezing condition can be rewritten as For incompressible systems the same equation remains valid when one sets η s = 0 From the discussion in the last section we know that this equation is equivalent to the simple conditionc 2D sub =c c , cf. (42). Comparing (50) with (49), it is easy to realize that the freezing equation in the compressible theory is equivalent to the one in the incompressible limit if one replaces c 2D sub (K 1 )[1 + η(K 1 )] byc 2D sub . Therefore, we can write a critical condition similar to (42) that is valid for a compressible system, where we display explicitly both the T -and the K 1dependences of the correlator. In (51), the correlator depends indirectly on the magnetic field B in the solid phase through K 1 via Consistently, since the magnetic field jumps across the transition (η = 0), the magnetic field in the liquid B l = Φ 0ρ = B/(1 + η) is different from the magnetic field B in the solid phase. Along the melting line B m (T ), K 1 can be written as a function of T only by using (52), i.e. K 1 = K 1 (B m (T )). Hence, at melting, the LHS of (51) can be written as a function of the temperature alone; taking the total derivative d/dT of (51) with respect to the temperature, we obtain We need to compute the derivative where we have used B = B l (1 + η) and the linearized relation between K 1 and G in (52). Inserting Eq. (53) (with the help of (54) and (31)) into (48), we obtain the entropy jump across the transition ∆s = 3µ 2ρ T dB l c 2D sub (T, K 1 ) 1 + η Inserting this expression into (47), we obtain the density jump η CC described by the Clausius-Clapeyron relation where in the last line we have used dB 2 l /4πρT = dΦ 0 B l /4πT = 1 − c 2D sub (0) from (27). The first term in the square brackets in (56) can be rewritten with the help of the freezing condition. As a final result we obtain which is exactly Eq. (44). Thus, we conclude that the saddle point equation (44) is fully consistent with the Clausius-Clapeyron relation.
To obtain an estimate for the density jump across the transition we have performed a numerical minimization of the constrained free energy (30). The system freezes when the solid minimum exhibits the same free energy as the liquid phase, i.e., when δω(η s , µ s ) = 0. In Fig.  4 we present the results of our numerical analysis for various values of T . For each T we show the density jump at the transition; as expected, we find a negative value of η. However, the modulus of the density jump η is always small, between |η| ≈ 10 −4 at low temperatures (large B) and |η| ≈ 10 −2 for large temperatures (low B).
The effect of such a small η on the determination of the melting line and on the value of µ s in the solid phase is negligible. Finally, from (47) we can obtain the value of the entropy jump across the transition. In agreement with previous studies 18 , we obtain a density jump per pancake vortex ∆S/N ≈ 0.4 k B at large fields decreasing to ∆S/N ≈ 0.1 k B at low magnetic fields corresponding to temperatures of order T = 0.2 ε 0 d.

IX. CONCLUSIONS
We have presented an analysis of the freezing transition of the magnetically coupled pancake vortex system within a classical density functional theory. Despite the simplicity of our approach, our results represent a considerable improvement when compared to the predictions of earlier work, particularly in the determination of the melting line. Moreover, we have addressed the problems with the thermodynamic inconsistency which affected earlier work. We showed how to obtain a negative density jump in accordance with the Clausius-Clapeyron equation and the retrograded melting line. Note that our derivation of the Clausius-Clapeyron equation is not limited to the present vortex system but can be generalized to non-magnetic systems.
The techniques described in this paper are easily extended to the study of inhomogeneous situations, as was already done in the analysis of surface effects on the melting transition 20 . The study of artificial surface pinning potentials offers itself as another application of our DFT approach. Such an analysis would shed light on the results of recent experiments described in Ref. 35 where the response of vortices in BiSCCO to a weak perturbation induced by pinning structures created on the sample surface has been investigated.