Scaling of the magnetic permeability at the Berezinskii-Kosterlitz-Thouless transition from Coulomb gas simulations

A new approach to the Berezinskii-Kosterlitz-Thouless transition in the two-dimensional Coulomb gas model is explored by MC simulation and finite size scaling. The usual mapping of a neutral two-dimensional superconductor in zero magnetic field to a Coulomb gas leads to an unscreened logarithmic interaction between the vortices, and with periodic boundary conditions vortex configurations are always vorticity neutral with an equal number of plus and minus vortices. We demonstrate that relaxing the neutrality condition has certain advantages. It leads to non-neutral vortex configurations that can appear in real systems with open boundary conditions and permits calculation of the compressibility, which for thin film superconductors corresponds to the magnetic permeability. The vortex-number fluctuation has remarkable scaling properties at and below the Berezinskii-Kosterlitz-Thouless transition. The fugacity variable becomes dangerously irrelevant in the low-temperature phase and leads to a multiplicative scaling correction to the mean-square vortex-number fluctuation and to the magnetic permeability. This multiplicative correction strongly affects the scaling properties of the vorticity fluctuation at and below the transition. Consequences of these findings are demonstrated using Monte Carlo simulations. Inclusion of the next-higher order correction to scaling is found to play an important role in the analysis of numerical data for the vortex number fluctuation and permits accurate determination of the critical properties.


INTRODUCTION
The Berezinskii-Kosterlitz-Thouless (BKT) transition [1][2][3] is a paradigm shifting phenomenon that demonstrates how topological fluctuations can create phase transitions without the Landau symmetry breaking mechanism. The topological excitations are quantized vortices and the BKT transition is a pair-unbinding transition, where the vortex correlations change long-distance behavior from algebraic to exponentially decaying. The BKT transition is found in several experimental ultrathin-film systems, e.g., superfluids, superconductors, Josephson-junction arrays, planar magnets, surface roughening, and melting [4].
Theoretical understanding of the BKT mechanism is provided by RG theory [5]. Notable signatures of the BKT transition are the universal jump of the superfluid density [6], the exponentially diverging correlation length, and the nonlinear current-voltage characteristics of ultrathin superconducting films [7,8].
In this paper we focus on a less frequently studied quantity, the magnetic permeability.
We study the scaling properties of the permeability at the BKT transition by Kosterlitz renormalization group theory and Monte Carlo simulations. Andersson and Lidmar [9] have pointed out that the vortex fugacity is dangerously irrelevant for the free vortex density in the low-temperature phase. Here we find that this also applies to the magnetic permeability. We further show that the fugacity enters as a multiplicative logarithmic correction to scaling that makes the permeability vanish at the transition in the thermodynamic limit. These results provide useful additional information that complements the characterization by superfluid transport, and should in principle be experimentally measurable. A multiplicative scaling correction with a similar origin is known for the two-point order parameter correlation function [5] and the XY magnetization susceptibility [10].
The BKT transition has been frequently studied in numerical simulations of XY models and Coulomb gas (CG) models. In this paper we test our predictions by large-scale Monte Carlo (MC) simulations of two-dimensional Coulomb gas models and find good agreement with theory. To locate the transition in the absence of symmetry breaking is somewhat complicated since the average of the order parameter vanishes across the transition and alternative methods must be invoked. The most common quantity considered in simulations of the BKT transition is the superfluid stiffness, or, in the CG language, the dielectric constant. To estimate the location of the transition from numerical data for finite systems a finite-size scaling method has to be used. Corrections to scaling are significant at the BKT transition and need to be included in a finite-size scaling analysis. The standard approach due to Weber and Minnhagen [11] is to include the lowest-order additive logarithmic correction to the universal jump of the superfluid stiffness and make a χ 2 fit of numerical data close to the transition to estimate the critical temperature from the best fit to the universal jump criterion. This approach has been used in a large number of simulation studies of the BKT transition; see Ref. [10] and references therein.
We show that the magnetic permeability µ is also a convenient quantity for locating the transition. We use a simplified version of the Weber-Minnhagen method, by means of an intersection analysis to extract the critical temperature from the finite size correction, without any need to use optimization to find the best fit. We focus on finite-size scaling of the magnetic permeability or, equivalently, the net vorticity fluctuation, which in the CG corresponds to the compressibility, and which is rendered finite by slight modification of the model. Our analysis clearly demonstrates the presence of a multiplicative log correction in the finite size scaling of our MC data and verifies the predictions of the theory. Furthermore, we find that the next order correction to scaling is substantial and demonstrate how it can be effectively included in the finite size scaling analysis. Applying the same analysis to the superfluid stiffness gives consistent results.
The paper is organized as follows. First we discuss scaling of the magnetic permeability at the transition. Then we turn to the two dimensional Coulomb gas model and the modifications needed to calculate the magnetic permeability. Then we describe our Monte Carlo method and the finite size scaling approach including lowest and next-order scaling corrections. Finally the simulation results are presented and discussed.
Here the coupling constant J is the superfluid stiffness, θ is the phase of the superconducting order parameter, A is the vector potential, B = ∇ × A the magnetic flux density, Φ 0 the flux quantum, h is an applied perpendicular magnetic field, and µ 0 is the permeability of free space. Vortices are topological defects in the phase field given by By following standard steps the model can be reformulated and expressed directly in the vortex coordinates [12]. The partition function then becomes where the Hamiltonian for the vortices takes the form of a CG model with interacting charges q n representing the vorticity. Here For finite screening length the BKT transition is in principle destroyed and replaced by a crossover. In real superconducting ultrathin films the screening length is finite but typically very large, comparable or exceeding the system size, meaning that the transition will still appear very sharp. The finite screening length implies that fluctuations in the magnetic flux density are present, leading to a finite magnetic permeability given by where d = 2 is the dimension. In terms of the flux of the B-field through a surface of area L 2 the total net vorticity is m = L 2 B/Φ 0 , and the permeability becomes In the CG language this corresponds to the compressibility κ = (1/L d T )( m 2 − m 2 ). From now on we consider the zero field case with B = m = 0. Below we reformulate the CG model on a lattice and study the permeability with Monte Carlo simulation.

SCALING RESULTS
In 3d bulk type-II superconductors the magnetic permeability changes from zero in the Thus, in general the permeability is expected to vanish as a power law at criticality and the net vorticity fluctuations to approach a constant at the transition. However, as will be discussed next this turns out not to be the case at the lower critical dimension d = 2.
To construct corrected scaling relations that apply for d = 2 a more accurate treatment is needed of vortex interactions and vortex number fluctuations. The scaling properties are described by renormalization group (RG) theory. The RG flow for the Coulomb gas is most easily expressed in terms of reduced temperature and fugacity variables defined as where l = ln b and b is a rescaling factor. The resulting RG flow obeys where the constant C is determined by the initial conditions. Below T c we have C 2 > 0 and above T c we have C 2 < 0. The BKT transition occurs for T = T c and then the flow obeys x + y = 0. Explicit solutions are given by where b 0 is determined by the initial condition for the RG flow. At and below the transition, the RG flow ends on the critical line x ≤ 0, y = 0.
According to the naive scaling Eq. (7), the mean square vortex density should approach a constant for b → ∞. This is not correct at the BKT transition where instead the mean square vorticity is expected to be proportional to the renormalized fugacity that scales to zero as b → ∞. Consider the renormalization of the magnetic permeability, The naive scaling would hold only if the right hand side tends to a constant as b → ∞, but this is not the case when y → 0, i.e., at and below T c . This is seen by explicit calculation of the partition function for m = 0, ±1 which yields for small z and thus goes to zero at the fixed line. For the permeability this gives µ ∼ y/L 2 . This modifies the naive scaling result m 2 ∼ constant and demonstrates that the fugacity gives a multiplicative correction to scaling. The multiplicative scaling correction makes the net vorticity asymptotically approach zero at the BKT transition. Stopping the RG flow at b ≈ L gives the finite size scaling formulas Right at T c this gives a multiplicative logarithmic correction, while below T c the naive power-law scaling is changed into one with a temperature dependent exponent.
The rest of the paper will study these relations by finite size scaling of data from Monte Carlo simulation. In the analysis of MC data the scale factor b will be taken to be the finite system size L, and the initial condition b 0 is an UV cutoff that corresponds to the vortex core radius L 0 .

MC SIMULATION OF THE LATTICE COULOMB GAS
To test the modified scaling relations described above we performed large-scale Monte Carlo simulations of a Coulomb gas model for vortex fluctuations. The lattice twodimensional Coulomb gas (CG) model is defined by [14,15] where q i = 0, ±1, ... is the Coulomb gas charge, or equivalently vorticity, on lattice site i of a square lattice with L × L sites with periodic boundary conditions. Here we only consider the case of no net applied magnetic field. The lattice Coulomb interaction is given by where k µ = 2πn µ /L, and we set the lattice spacing to a = 1. µ v is the vortex chemical potential, E c = −µ v is the vortex core energy, and z = e βµv is the fugacity. The particle number or total vorticity is N = i |q i | = N + + N − , and the net vorticity is m = i q i = where β = 1/T . We studied a few different values of the vortex chemical potential µ v . All where λ is the screening length. Here B(r) correspond to the magnetic flux through a plaquette in units of Φ 0 . Since V (k = 0) = ∞, uniform fluctuations in B(k = 0) are accompanied by fluctuations in the net vorticity such that q(k = 0) = B(k = 0). Integrating out the fluctuations in B(k) leads to a screened Coulomb interaction given by For finite λ the self energy is finite and therefore charged configurations appear so that m 2 can be studied. We considered two different models of screening. The first model includes The value of the screening length λ needs to be selected in a special way to guarantee that the CG has a BKT transition. The problem is that the screened models described The main focus of this paper is to study how the BKT transition is seen in the finite size scaling properties of MC data for various quantities. For the neutral case with m = 0 (and λ = ∞), the BKT transition can conveniently be located from simulation data using the Nelson-Kosterlitz universal jump [6] of the superfluid density ρ s at the transition. In the CG language this corresponds to a universal jump in the dielectric response function 1/ = 2πρ s /ρ 0 that will be considered here. At the BKT transition temperature T c the dielectric response function −1 jumps from the finite value to zero at the transition. To calculate the dielectric function in a MC simulation of the CG it is useful to add a polarization term to the energy [16,17] where P = (P x , P y ) = i q i r i is the polarization. The dielectric function is then obtained from the polarization fluctuation by [16,17] The We refer to one sweep through the system as L × L update attempts to insert dipoles or single charges. Other types of MC moves like moving or removing particles or pairs can also be used and potentially improves convergence properties but we settled with the moves described above since they gave satisfactory convergence of the simulation. We found that 10 3 initial sweeps to establish equilibrium was sufficient. In equilibrium measurements of the observables were done after each sweep. A total of 10 4 terms were collected to form averages. The runs were repeated about 10 2 times until sufficiently small statistical errors had been obtained. Error bars were estimated from the standard deviation of the results from different runs. Single histogram reweighting was used to obtain data at a range of nearby temperatures from simulations done at T = 0.2115 [18].

RESULTS
As a first step we discuss how the BKT transition temperature can be estimated from MC data using the universal jump criterion. In this calculation we use the CG model without any fluctuations in the net vorticity so that m = 0 throughout the simulation. Figure 1 shows MC data for 1/ obtained by evaluating Eq. to the presence of big finite size effects, no clear indication of a jump is seen in the data for finite system sizes. Instead, as the system size increases we expect the data curves to slowly approach the characteristic square root cusp at T c and undergo the universal jump [7]. As seen in the figure the approach to the asymptotic behavior is very slow which complicates the estimation of T c .
The problem with the slow approach to the universal jump is overcome by the Weber-Minnhagen finite size scaling form of the approach to the universal jump given by [11] where C = −2 ln b 0 is an unknown constant. This form follows from Kosterlitz RG equations and gives the leading additive logarithmic finite-size correction to the universal jump value.
A common approach to estimating the transition temperature T c is to minimize the RMS deviation between numerical data and the finite size form over variations in both T c and C.
Here we propose a simpler method. If the universal jump value is assumed to be correct, solving Eq. (22) for the constant C gives A plot of MC data curves for 1 1/4 T −1 −2 ln L vs T for different system sizes L is thus expected to produce a system-size independent intersection point at (T c , C). This procedure involves no parameter fitting and straightforwardly produces an accurate estimate.  where subscript L denotes data for system size L. Hence this quantity is expected to be independent of system size at the transition. Figure 5 shows the corresponding MC data curves. An intersection point for large system sizes is found at T c ≈ 0.2115 which agrees with the value obtained for the dielectric function in Fig. 2. However, the intersection quality is not nearly as good and finite size effects are substantial despite including the correction to scaling. In a sense it is not unexpected that the approach to the transition is slower in this case than in the neutral case since we introduced another length scale λ ∝ L in the problem.
In addition the intersection point gives the estimate b 0 ≈ e −B/A ≈ 0.151, which is similar to the value 0.135 found from dielectric constant above. Also note that a similar subtraction method enables a direct test of the size of the universal jump in the dielectric function.
Since the intersection points in Fig. 5 drift with system size, it is motivated to attempt to include the next order correction in the finite size scaling analysis. The form of the next order correction for large length scales is known and for the dielectric constant becomes [10,19,20] 1/4 T c = 1 + 1/2 ln L/L 0 + const · ln ln L/ ln 2 L. Adopting this form of the correction to the charge fluctuation at T c gives where A, B, C are unknown constants. Compared to the case of fitting the dielectric function where the known universal jump could be used, the charge fluctuation involves one more unknown constant, which complicates fitting to numerical data. A possible approach is to look for the parameters that give the best χ 2 -fit to the data. Instead we construct a simple intersection plot that contains the same information. Equation (26) gives where f L = C ln ln L/ ln 2 L. Subtraction of this quantity for pairs of system sizes (2L, L)

eliminates B and gives
Hence a plot of this quantity is expected to produce a system size intersection at T c when the optimal value of C is used. Finally we investigate the scaling properties of m 2 below T c . According to Eqs. (6), (13) the size dependence here becomes a power law with a temperature dependent exponent given by m 2 ∼ L 2x R (T ) . Figure 7 shows MC data points for m 2 plotted vs L for temperatures T < T c . The solid lines are power law fits to the MC data points. Figure 8 shows the fitted power law exponent x plotted vs 1 − 1/4T . The theoretical prediction discussed above is   Furthermore, the next-order scaling correction is needed to get an accurate estimate of the critical point and we implemented this calculation in an intersection analysis. We found excellent agreement between the different methods for locating the transition. Finally we showed that Kosterlitz RG theory leads to a power-law dependence of the charge fluctuation with system size which is to a good approximation consistent with our MC data.
In summary we consider non-neutral vortex fluctuations as a useful quantity for studying the BKT transition and show that it has a novel multiplicative log correction to scaling at the transition. It would be interesting to look for such properties in experiments by for example magnetic permeability measurements in effectively two-dimensional superconductors and in cold atom systems.