# Phase-transitions in isotropic extreme type-II superconductors

## I Introduction

Since the discovery of copper oxide based high temperature superconductors (HTSC) , which are of the extreme type-II variety, there has been great interest in their phenomenological phase diagrams. Abrikosov’s mean-field description, which is valid for conventional low temperature superconductors, is expected to be modified by the strong thermal fluctuations in HTSC. Extensive research in both theory, numerical simulations, and experiments over the years have resulted in a general consensus on some of the fundamental issues. The current understanding of HTSC in a uniform magnetic field $`𝐁`$ is as follows. In the absence of any pinning disorder, the low temperature Abrikosov vortex lattice phase melts into a vortex liquid via a first order transition at the temperature $`T_m`$. The Abrikosov vortex lattice phase is characterized by a transverse triangular crystalline order and a finite longitudinal phase coherence.

However, recently there has been some debate about the nature of the vortex liquid which the Abrikosov vortex lattice melts into as temperature increases. Numerous simulations using the 3DXY model , London and Lowest Landau Level approximations have indicated that the vortex liquid is incoherent, i.e. the phase coherence or superfluid density in any direction, is zero. The crystalline order and phase coherence are destroyed simultaneously at the melting transition. This scenario has been supported by experiments on high quality YBCO crystals. Other simulations using the 3DXY model  have suggested that the longitudinal phase coherence persists above the melting transition and only vanishes at a higher “entanglement” temperature $`T_E`$. In this scenario, the vortex liquid at $`T_m<T<T_E`$ would be disentangled with relatively straight vortex lines. For $`T>T_E`$, it was claimed that larger thermal fluctuations causes the vortex lines to be entangled with a concomitant loss of global phase coherence.

Previously, a vortex liquid with nonzero longitudinal phase-coherence, or superfluid density, was found in simulations on an isotropic system. More recently , it has been proposed that for large $`B`$, $`T_m`$ and $`T_E`$ merge into a single transition whereas for small $`𝐁`$ they are well separated. The authors found that for an isotropic system, the crossover between the two regimes is at a filling fraction $`1/18<f<1/6`$, where the precise definition of $`fB`$ will be given shortly. However, these simulations were performed on rather small systems. In the present work, we have performed similar simulations on a much larger system, and found an exclusively incoherent vortex liquid down to $`f=1/60`$. This implies that the only thermodynamically stable vortex liquid phase, is one which has zero longitudinal superfluid density, with full translational and rotational symmetries and zero phase coherence in all directions. We are therefore led to the conclusion that the recently discussed  crossover between the phase-coherent vortex liquid and the phase-incoherent vortex liquid is a numerical artifact.

The organization of this paper is as follows. In Section II we first briefly review the model used in the simulations, and the quantities we calculate. In Section III we discuss our results and their implications. In Section IV we provide a summary and conclusion.

## II Model

We use a uniformly frustrated 3DXY model on a cubic lattice to describe an isotropic, extreme type-II superconductors in a magnetic field. The London model for superconductors can be readily derived from the phenomenological Ginzburg-Landau (GL) model with the approximation that the amplitude of the local complex order parameter $`\psi =|\psi |\mathrm{exp}[i\theta ]`$, is fixed. The Hamiltonian $`H`$ of the London model consists of degrees of freedom in the phase $`\theta (𝐫)`$ and the gauge vector potential $`𝐀_{vp}(𝐫)`$ associated with the magnetic induction $`𝐁`$ inside the system, i.e. $`H=H[\{\theta (𝐫),𝐀_{vp}(𝐫)\}]`$.

For an isotropic extreme type-II superconductor, the penetration depth $`\lambda `$ is much larger than the coherence length $`\xi `$ such that the GL parameter $`\kappa =\lambda /\xi \mathrm{}`$. This means that the magnetic fields surrounding the vortex lines strongly overlap with one another giving a spatially smooth $`𝐁`$. This condition is ensured in the regime where $`𝐁>𝐁_{c1}`$. In other words, the fluctuations of $`𝐀_{vp}(𝐫)`$ on the length scale of $`\xi `$ are negligible compared to the fluctuations of $`\theta (𝐫)`$. Therefore, we can further simplify the problem by dropping the degrees of freedom of $`𝐀_{vp}(𝐫)`$ from the Hamiltonian $`H`$, and fix $`B`$ equal to the external applied magnetic field. The resulting Hamiltonian $`H=H[\{\theta (𝐫)\}]`$ is the 3DXY model. This $`\kappa \mathrm{}`$ (or frozen gauge) approximation has been widely used as a phenomenological model for superconductors. Note that, within this model, the system has no magnetic flux lines, since there can be no tubes of confined magnetic flux when $`\lambda \mathrm{}`$. The system only exhibits vortex lines.

In order to perform numerical simulations on the resulting model, we discretize the model on a 3D cubic lattice with grid spacing $`\xi `$. The dimensionless Hamiltonian of this lattice model is given by 

$`H[\{\theta (𝐫)\}]=J_0{\displaystyle \underset{𝐫,\alpha =x,y,z}{}}\mathrm{cos}\left(j_\alpha [\theta (𝐫)]\right)`$ (1)   
where $`J_0=\mathrm{\Phi }_0^2\xi /4\pi ^2\mu _0\lambda ^2`$ is the isotropic coupling energy, $`\mathrm{\Phi }_0=\mathrm{}\pi /e`$ is the flux quanta and $`\mu _0`$ is the permeability of a vacuum in SI units. The dimensionless vector $`𝐫`$ labels the position of an arbitrary grid point. We define the gauge-invariant phase difference $`j_\alpha (𝐫)`$ as

$`j_\alpha (𝐫)\theta (𝐫+𝐞_\alpha )\theta (𝐫){\displaystyle \frac{2\pi \xi }{\mathrm{\Phi }_0}}{\displaystyle _𝐫^{𝐫+𝐞_\alpha }}𝑑𝐫^{}𝐀_{vp}(𝐫^{})`$ (2)   
where $`𝐞_\alpha `$ is the unit vector along the $`\alpha `$-axis. The convention is that $`j_\alpha (𝐫)`$ “flows” from the grid at $`𝐫`$ to $`𝐫+𝐞_\alpha `$. In this model, it is natural to define a dimensionless temperature $`\stackrel{~}{T}=k_BT/J_0`$. The magnetic induction $`𝐁`$ is conveniently represented by the filling fraction

$`f={\displaystyle \frac{B\xi ^2}{\mathrm{\Phi }_0}}.`$ (3)   
In this paper, we perform a Monte Carlo simulation on an isotropic system with $`V=N_x\times N_y\times N_z`$ number of grid points. The aim of our numerical simulations is to identify possible thermodynamic phases (i.e. vortex lattice or liquid, etc.) and phase transitions (first order or continuous) or crossovers associated with them. Specifically, we are interested in calculating the internal energy, specific heat, structure factor, and helicity moduli. These thermodynamic quantities and their physical significance will be discussed next. In addition, we introduce and define a quantity $`O_L`$ which we denote the vortex-path probability, and discuss some of its physical implications.

A first order transition is indicated by a $`\delta `$-function anomaly in the specific heat, equivalently a discontinuous jump in the internal energy. On the other hand, the hallmark of a continuous transition is a jump in the specific heat, modified by fluctuation contributions to the anomaly, and a continuous internal energy. The internal energy per site is obtained by averaging $`H`$ in Eq. (1) over the thermal equilibrium states, normalized by the total number of grids, i.e. $`E=<H>/V`$. We define a dimensionless specific heat per site $`C`$ using the standard fluctuation theorem of a classical system with Gibbs distribution:

$`C={\displaystyle \frac{<H^2><H>^2}{V(k_BT)^2}}`$ (4)   
A convenient and widely used quantity to probe the global phase coherence of the system is the helicity modulus, which is proportional to the second derivative of the free energy associated with Eq. (1) with respect to an infinitesimal phase twist. On a square lattice, the dimensionless helicity modulus $`\mathrm{{\rm Y}}_\alpha `$ can be written as

$`\mathrm{{\rm Y}}_\alpha ={\displaystyle \frac{1}{V}}<{\displaystyle \underset{𝐫}{}}\mathrm{cos}j_\alpha {\displaystyle \frac{1}{\stackrel{~}{T}}}\left|{\displaystyle \underset{𝐫}{}}\mathrm{sin}j_\alpha \right|^2>.`$ (5)   
$`\mathrm{{\rm Y}}_\alpha `$ measures the “stiffness” of $`\theta `$ with respect to twisting along the $`\alpha `$-direction. If $`\mathrm{{\rm Y}}_\alpha >0`$ then $`\theta `$ is stiff in the $`\alpha `$-direction, or more precisely, there is global phase coherence or superconducting response in that direction. By tracking the temperature at which $`\mathrm{{\rm Y}}`$ goes to zero, one can determine the superconducting-normal metal transition temperature of the system. In the mixed phase $`\mathrm{{\rm Y}}_{}>0`$ and $`\mathrm{{\rm Y}}_{}=0`$, where the subscripts $``$ and $``$ denote the directions longitudinal and transverse to $`𝐁`$, respectively.

The structure factor probes the transverse crystalline order of the vortex system. We adopt a conventional definition :

$`S(𝐤_{})={\displaystyle \frac{1}{(fV)^2}}<\left|{\displaystyle \underset{𝐫}{}}v_{}(𝐫)\mathrm{exp}(i𝐤_{}𝐫)\right|^2>`$ (6)   
where $`𝐤_{}`$ is a two-dimensional reciprocal vector and $`v_\alpha (𝐫)`$ is the local vorticity measured on the dual square lattice grid , composed of the centers of every direct unit cell. A crystalline ordered phase is characterized by $`S(𝐤_0)>0`$ (or Bragg peaks) where $`𝐤_0`$ are the discrete set of reciprocal lattice vectors associated with the crystal structure. On the other hand, $`S(𝐤)`$ for a phase with full rotational invariance exhibits ring patterns.

At a fixed temperature, the equilibrium configurations are generated by making random changes to $`\theta (𝐫)`$ at each grid-point via the the Metropolis algorithm. This is equivalent to randomly changing all the six $`j_\alpha `$ attached to each grid-point. To ensure conservation of vorticity in each unit cell, $`j_\alpha `$ has to be in the range $`[\pi ,\pi )`$. Addition of $`\pm 2\pi `$ shall be administrated to bring $`j_\alpha `$ back into range at every Monte Carlo step. This procedure introduces vortex-loops into the system. Such loops are the elementary topological excitation of the model. Periodic boundary conditions (PBC) are imposed on $`j_\alpha `$ such that $`j_\alpha (𝐫+N_\mu e_\mu )=j_\alpha (𝐫)`$ for $`\mu =(x,y,z)`$. More details about the Monte Carlo procedure may be found in Refs..

We define $`O_L`$ as the probability of finding a directed vortex path threading the entire system transverse to the induction $`B`$, without using the PBC along the field direction . It is obtained by computing the number $`N_V`$ of times we find at least one such path threading the system in any direction $`B`$ in $`N_P`$ different phase-configurations, normalized by $`N_P`$, i.e. $`O_L=N_V/N_P`$. The fact that $`O_L=0`$ implies that there is no connected path of vortex segments that threads the entire system in the transverse direction, without using PBC along the field direction several times. Now, let $`N_L^\alpha `$ ($`\alpha [x,y,z]`$) denote the areal density of connected vortex paths threading the system in any direction, including the direction parallel to the induction. It is clear that in the Abrikosov vortex lattice phase $`O_L=0`$, and $`N_L^z=B/\mathrm{\Phi }_0`$, while $`N_L^x=N_L^y=0`$. Thus, $`N_L^\alpha `$ is a conserved quantity at fixed induction $`B`$. On the other hand, $`O_L=1`$ implies that $`N_L^{x,y}>0`$, and the total number of vortex paths threading the system in any direction scales with system size, but undergoes thermal fluctuations. Therefore, $`N_L^\alpha `$ is no longer a conserved quantity.

We propose the following scenario to interpret the change in $`O_L`$ and $`N_L^\alpha `$. Number conservation uniquely identifies a $`U(1)`$-symmetry, and hence the low-temperature phase of the vortex-system (the dual of the phase-representation of the superconductor) exhibits explicit $`U(1)`$-symmetry, since $`O_L=0`$. At high temperatures $`O_L=1`$, $`N_L^\alpha `$ is not conserved, and the $`U(1)`$-symmetry is broken. A $`U(1)`$-symmetric phase cannot be analytically continued to a $`U(1)`$-nonsymmetric one. The change in $`O_L`$ from $`0`$ to $`1`$ could therefore signal a phase-transition, in this case involving breaking a global $`U(1)`$-symmetry, in finite as well as zero magnetic field. However, to substantiate such a claim, one needs to argue that $`O_L`$ is related to a local order parameter of the system. To this end, we note that it is possible to transcribe the vortex-part of the Ginzburg-Landay theory in the phase-only approximation in such a way that the vortex-part of the theory is specified in terms of a local, complex matter field $`\varphi (𝐱)`$, and that the theory then explicitly exhibits a $`U(1)`$-symmetry . In the Lattice London Model, corresponding to the Villain-approximation to the Ginzburg-Landau theory, this symmetry is therefore only implicit, or “hidden”. The probability of finding a connected vortex-path starting at a point $`𝐱`$ and ending at point $`𝐲`$, $`G(𝐱,𝐲)`$, is given in terms of the two-point correlation function of the matter field $`\varphi (𝐱)`$, $`G(𝐱,𝐲)=<\varphi ^{}(𝐱)\varphi (𝐲)>`$ . $`O_L`$ may be viewed as a special case of $`G`$, and in the thermodynamic limit corresponds to $`lim_{|𝐱𝐲|\mathrm{}}G(𝐱,𝐲)`$. If $`G(𝐱,𝐲)`$ is non-zero in this limit, this suggests the possibility of having $`<\varphi (𝐱)>0`$, and hence a broken $`U(1)`$-symmetry. Although this does not constitute a proof that $`O_L`$ is connected to a local order parameter whose expectation value is associated with a broken $`U(1)`$-symmetry, it seems to be suggestive of such a phase-transition existing even at finite magnetic field. Note that the above local matter field $`\varphi `$ appears to be the dual field of a complex order parameter appearing in a somewhat different independent approach to the same problem by Tešanović .

## III Results and Discussion

In this section, we discuss the results of our simulations on an isotropic system with size of $`120^3`$ grid points, which is the largest to date on an isotropic system. The system is subdivided into multiple sections, and the Monte Carlo procedure is implemented in parallel across the sections using 3D “black and white” decomposition. The filling fraction considered are $`1/f=20,40,60`$, and $`\mathrm{}`$. The system is cooled from high temperatures. For each temperature, a typical run consists of 120000 Monte Carlo sweeps across the whole lattice, $`1/3`$ of that were used for equilibration. Near phase-transitions, up to 600000 sweeps were used.

In Fig. (1), we present the results for the zero-field case, $`f=0`$, where we show the specific heat $`C`$, the superfluid stiffness $`\mathrm{{\rm Y}}`$, and the vortex-path probability $`O_L`$ as functions of temperature. Note how the specific heat anomaly, the vanishing of the superfluid stiffness and $`O_L`$, all coincide with $`\stackrel{~}{T}_c2.20`$. The physical interpretation of $`O_L=1`$ is that the effective long-wavelength vortex-line tension vanishes . This claim is substantiated by calculating the probability distribution $`D(p)`$ of vortex loops as a function of perimeter $`p`$, at various temperatures. We may fit this distribution to the form 

$`D(p)`$ $`=`$ $`Ap^\alpha e^{\beta \epsilon p};T<T_c`$ (7)   
$`=`$ $`Ap^\alpha ;TT_c`$ (8)   
where $`\epsilon `$ is an effective, temperature dependent, long-wavelength line tension, $`\beta =1/k_BT`$, and $`\alpha =5/2`$. The results are shown in Fig. (2), demonstrating that at $`T=T_c`$, a purely algebraic decay is realized, implying $`\epsilon =0`$. Below $`T_c`$, the exponential decay is well fitted, showing that $`\epsilon 0`$. The inset of Fig. (2) shows $`\epsilon `$ as a function of temperature, as obtained in our simulations. A similar method of extracting $`\epsilon `$ for the zero-field case using a similar form for $`D(p)`$ (without he power-law prefactor) has previously been used by Li and Teitel, in Ref. .

The results for the quantities $`S(𝐤)`$, $`\mathrm{{\rm Y}}`$ and $`C`$ for $`1/f=20`$ are shown in Fig. (3). Similar qualitative features are also found for the cases $`1/f=40,60`$. First, $`S(𝐤_{})`$ exhibits six-fold Bragg peaks at low temperatures (not shown). This indicates that the low temperature phase is a triangular vortex lattice phase. The destruction of the vortex lattice structure is marked by a melting transition at $`\stackrel{~}{T}_m1.34`$ where $`S(𝐤_1)`$ drops sharply to zero, where $`𝐤_1`$ is the wavevector of one of the six first order Bragg peaks. The sharpness of the drop in $`S(𝐤_1)`$ strongly suggests that the transition is first order. This is confirmed by the appearance of a $`\delta `$-function like peak in $`C`$ at the same temperature $`\stackrel{~}{T}_m`$. Coincidentally, $`\mathrm{{\rm Y}}_{}`$ which is finite at low temperatures, also drops sharply to zero at $`\stackrel{~}{T}_m`$. The isotropic system exhibits longitudinal superconductivity below $`\stackrel{~}{T}_m`$, but not above. Moreover, based on our above discussion in Section II, we interpret the rise in the quantity $`O_L`$ from $`0`$ to $`1`$ as signalling that the effective long-wavelength vortex-line tension vanishes. Thus, the vortex liquid phase is divided into two regions in phase-space. In one region, the vortex liquid is phase-incoherent i.e. has no longitudinal superfluid density but has finite vortex-line tension. We propose that this phase exhibits a vortex-associated $`U(1)`$-symmetry. In the other region the vortex liquid is phase-incoherent but has zero vortex-line tension. We propose that this phase exhibits a broken $`U(1)`$-symmetry.

Based on this, one can conclude that a first order melting transition of the Abrikosov vortex lattice exists in an isotropic system in the absence of pinning. The entire vortex liquid phase is incoherent, i.e. a vortex liquid phase with no longitudinal superconductivity. Note that the same conclusion has been reached in earlier simulations on anisotropic systems. There have been earlier reports of a disentangled vortex liquid, i.e. vortex liquid with non-zero longitudinal phase coherence. These results have been obtained in similar simulation. In Ref. , it was argued that a phase-coherent vortex liquid should be most easily observed in isotropic systems for $`f1/18`$. However, these simulations were performed on comparatively small systems, typically of sizes no larger than $`24^3`$. Our present results are based on much larger systems, $`120^3`$. We have observed a slight difference in the temperatures at which $`S(𝐤)`$ and $`\mathrm{{\rm Y}}_{}`$ vanish in simulations on smaller system or lower number of sweeps. Therefore, it may be concluded that the existence of a vortex liquid with a nonzero longitudinal superfluid density is a numerical artifact of small system sizes and/or insufficient simulation time . It is well understood that the $`\lambda `$-transition is driven by proliferation of thermally excited loops of all sizes. In recent work by some of us, it was proposed that the same mechanism is driving the first order melting transition at low $`B`$ and a newly discovered continuous transition involving the breaking of U(1) symmetry at large $`B`$.

Similar results for $`1/f=40`$ and $`60`$ enable us to propose a simple phase diagram for an isotropic extreme type-II superconductor in the absence of disorder (see Fig. (4)). Immediately below and above the $`\stackrel{~}{T}_m`$ line, the phases of the vortex system are identified as the Abrikosov vortex lattice and the incoherent vortex liquid, respectively.

In a pin-free system, one would expect $`\mathrm{{\rm Y}}_{}`$ to be zero at all temperature. In this case, the numerical lattice, on which simulations are performed, effectively pins the vortex lines from moving in the transverse plane and counteracts the Lorentz forces on them in the presence of a transverse applied current. However, the pinning is overcome by thermal fluctuations at higher temperature and the vortex lines are depinned at temperature $`\stackrel{~}{T}_d0.62`$. Fortunately, we see that $`\stackrel{~}{T}_d\stackrel{~}{T}_m`$, which means that near $`\stackrel{~}{T}_m`$, the vortex lines, and the melting process, are completely free from the effects of the numerical grid. Therefore, the features of $`C`$, $`\mathrm{{\rm Y}}_{}`$ and $`S(𝐤)`$ at $`\stackrel{~}{T}_m`$ are genuine thermodynamic effects.

## IV Conclusion

We have performed simulations of the uniformly frustrated 3DXY model on a large isotropic system ($`120^3`$ grid points) for a variety of filling fractions $`1/f=20,40,60`$ and $`\mathrm{}`$. We found a first order melting transition in this isotropic system for all the three non-zero values of $`f`$ considered. The longitudinal phase coherence and triangular crystalline order of the Abrikosov vortex crystal are simultaneously destroyed at the melting transition. Above the melting temperature, the incoherent vortex liquid is the only thermodynamic phase. We have demonstrated that previous claims of the existence of disentangled vortex liquid is due to performing simulations using insufficient system sizes and simulation times.

We have shown that the effective vortex-line tension vanishes precisely at the zero-field superconducting transition. The loss of superfluid stiffness, the loss of line tension, and the abrupt change in the connectivity of the vortex tangle, as signalled by the change in the quantitity $`O_L`$ across the system, all coincide in this case. A similar change in connectivity across the vortex system takes place at finite magnetic field. The results of the Ref. , the present paper, and in particular those of Ref.  strongly indicate that this change in connectivity is sharp in the limit of large systems, thus indicating the loss of number-conservation of connected vortex paths threading the system. Since the finite probability of finding a connected vortex path threading the system in a direction other than the magnetic field may be tied to the finite expectation value of a local complex matter field , this lends further support to the argument that the change in $`O_L`$ signals the breaking of a $`U(1)`$-symmetry . At the very least this proposition appears to be intriguing enough to warrant further investigation.

We finally caution the reader that we so far have not been able to detect any anomaly in specific heat at the suggested new finite-field transition inside the vortex liquid, for the isotropic case. Even the anomaly at the first order melting transition is weak in the isotropic case, and is difficult to bring out in simulations. It may be that considerably larger systems are needed for the isotropic case in order to see signals in the specific heat of the suggested new transition due to the small amount of entropy in the transition. This is the reason why the anomaly at the first order vortex lattice melting transition is difficult to observe in simulations. It is conceivable that increasing the anisotropy of the system should bring out the anomaly clearer, if it exists. This indeed is the case for the anomaly at the vortex lattice melting transition. A weak anomaly associated with the putative $`U(1)`$-transition, may in fact have been observed for the anisotropic case in Ref. .

## Acknowledgment

We thank J. S. Høye, P. Olsson, S. Teitel, and Z. Teśanovič for discussions and consrtructive remarks. This work is supported by the Norges Forskningsråd under grants 115577/410, 110566/410 and 110569/410. We acknowledge the use of the Cray T3E-600 at the Norwegian Supercomputing Project of NTNU, Trondheim, Norway.