Interplay of spin waves and vortices in the 2D XY model at small vortex-core energy

The Berezinskii-Kosterlitz-Thouless (BKT) mechanism describes universal vortex unbinding in many two-dimensional systems, including the paradigmatic XY model. However, most of these systems present a complex interplay between excitations at different length scales that complicates theoretical calculations of nonuniversal thermodynamic quantities. These difficulties may be overcome by suitably modifying the initial conditions of the BKT flow equations to account for noncritical fluctuations at small length scales. In this work, we perform a systematic study of the validity and limits of this two-step approach by constructing optimised initial conditions for the BKT flow. We find that the two-step approach can accurately reproduce the results of Monte-Carlo simulations of the traditional XY model. In order to systematically study the interplay between vortices and spin-wave excitations, we introduce a modified XY model with increased vortex fugacity. We present large-scale Monte-Carlo simulations of the spin stiffness and vortex density for this modified XY model and show that even at large vortex fugacity, vortex unbinding is accurately described by the nonperturbative functional renormalisation group.


I. INTRODUCTION
More than forty years after its observation in thin 4 He films [1], the Berezinskii-Kosterlitz-Thouless (BKT) transition [2][3][4][5] has registered in recent years an increasing interest in the field of low-dimensional strongly-correlated electron systems. In the last decade, the class of systems in which the BKT transition has been detected has become remarkably wider, including new low-dimensional experimental systems such as quasi-2D layered superconductors [6][7][8][9][10], exciton-polariton systems [11], cold atoms in 2D harmonic traps [12,13] and 2D electron gases [14][15][16] at the interface between insulating oxides in artificial heterostructures. At the same time, interesting new issues appeared also in conventional superconducting (SC) films. For instance, the spatial inhomogeneity of the superconducting order parameter has been observed to become more pronounced near the superconducting-insulator transition (SIT) [17,18], raising the question of its effect on the BKT phase transition [19][20][21].
Theoretically, the BKT scaling is the key to understand the universal features of a wide range of natural phenomena ranging from DNA tangling in biology to pattern formation in complex systems [22][23][24]. The hallmark of such a scaling is the expected discontinuous jump of the phase stiffness J s from a finite value J s (T BKT ) at the BKT temperature T BKT to zero right above it. According to the famous Nelson relation [25], the value J s (T BKT ) is a universal function of the critical temperature itself, i.e., Despite the apparent simplicity of this relation, the predictive power of Eq. (1) to infer the temperature dependence of the phase stiffness and to establish the critical temperature in real microscopic systems is far from obvious. Indeed, in the original derivation based on the celebrated Renormalization Group (RG) BKT equations [4,25], Eq. (1) only accounts for the role of vortex-like excitations to drive the transition, including the possible renormalization effects due to bound vortex-antivortex pairs which proliferate already below T BKT . First, in order to quantify this effect one needs a precise knowledge of the vortex-core energy µ v , which controls the vortex fugacity g = exp(−µ v /k B T ). Indeed, for larger vortex fugacity it is easier to generate vortex pairs at short length scales; these suppress the phase rigidity without destroying it globally and reduce T BKT . A second issue in real systems is the presence of additional non-topological phase excitations that also deplete the phase rigidity. How these additional excitations contribute, cooperate or interfere with vortex-like excitations to determine the global T BKT is far from understood. The paradigmatic example for the difficulty to estimate T BKT from Eq. (1) is provided by the XY model itself, where the BKT transition was initially discussed [2][3][4]. The XY model also admits longitudinal phase fluctuations (or spin waves in the spin language), which are effective at low temperature to deplete the stiffness linearly in temper-ature even while vortex-like fluctuations are thermally suppressed. To deal quantitatively with this effect the approach pursued so far in the literature relied in a two-step procedure: first one derives effective low-energy couplings for the stiffness and fugacity, which then serve as renormalised initial conditions for the BKT flow equations. Since the work of J. Villain within the context of 2D classical planar magnets [26], several efforts have been undertaken to derive the universal BKT scaling directly from the microscopic variables of each model, but no consistent picture has emerged despite the existence of solvable models that display BKT scaling [27]. To account for the non-topological phase modes of the system, in addition to the RG approach, different theoretical techniques have been proposed, such as the self-consistent harmonic approximation [26,28], classical Monte Carlo simulations [29], Gaussian fluctuations [30,31] and the functional renormalization group [32,33].
A second example of the relevance of nonuniversal effects is the application to thin films of superconductors. In this case, it has been suggested [34][35][36][37][38] that the two-step procedure should include the effect of quasiparticle excitations across the superconducting gap on the stiffness, which represent the most relevant source of depletion of the superfluid density in a superconductor. However, even including this effect the comparison with experiments has shown that the vortex-core energy µ v is significantly smaller than expected within the XY model. Indeed, within a microscopic BCS picture for superconductors µ v is expected not to scale with the superfluid stiffness, but rather with the Cooperpair condensation energy E c . On the other hand, within the XY model, the ratio µ v /J s has been estimated [3] long ago to be a constant µ v /J s = π 2 /2 that depends only on the lattice structure, as expected for a model with a single coupling constant. Once more, a quantitative check of the accuracy of such an estimate for the XY model is still lacking. More broadly, the case of superconducting films calls for a deeper understanding of the two-step RG approach for smaller values of the vortex fugacity. In particular, the interplay between transverse and longitudinal phase fluctuations is expected to become more relevant as the lowering of µ v favours thermal vortex-like excitations at short length scales.
What emerges clearly from these examples is that despite the universality of the Nelson criterion (1), the critical temperature T BKT itself is not universal, and it depends on microscopic details which must be properly included in the RG flow to correctly reproduce the temperature dependence of the phase stiffness. In this paper, we will systematically address this issue within the paradigm of a generalized 2D XY model with tuneable vortex-core energy. Our aim is to accurately test the two-step procedure by comparing the numerical results from Monte Carlo (MC) simulations with the RG two-step procedure for increasing values of the vortex fugacity. To this end, we first extract the stiffness and vortex-core energy from low-temperature MC data and then use them as initial conditions for the subsequent RG flow. This procedure turns out to be rather successful for the standard XY model, with an excellent estimate of T BKT , despite some small discrepancy in the temperature evolution of the superfluid stiffness near the transition. For increasing vortex fugacity the accuracy of the original RG equations decreases, but a good estimate of T BKT can still be obtained by considering either higher-order terms in the vortex fugacity [39], via the self-consistent approach suggested by Timm [40] or via the nonperturbative functional renormalization group (FRG). We find that FRG provides an excellent estimate of T BKT , even though it cannot fully capture the temperature dependence of J s (T ) obtained from Monte Carlo simulations. This discrepancy can be an indication of the interplay between longitudinal and transverse phase fluctuations, which cannot be captured by the two-step RG approach.
The paper is organized as follows. In Sec. II we introduce the two-step procedure for the classical XY model, providing a new estimate of the vortexcore energy µ v based on Monte Carlo numerical results. In Sec. III we introduce the modified 2D XY model studied in this work and we present the MC results. In Sec. IV we study the modified XY model via a two-step procedure, where an effective value for both the vortex-core energy µ eff v and the superfluid stiffness J eff are inserted into the BKT flow equations. Finally, we conclude in Sec. V.

II. THE VORTEX-CORE ENERGY OF THE XY MODEL
Ever since the seminal papers by Berezinskii, Kosterlitz and Thouless [2][3][4], the classical 2D XY model has been the reference model for the study of the BKT phase transition. The model describes the ferromagnetic nearest-neighbor interaction between planar spins with fixed modulus (| S i | = 1), interacting via a ferromagnetic spin-spin coupling constant J > 0: where θ i is the angle the i-th spin forms with a given direction. Within the context of SC films, one can map S i → |∆|e iθi so that Eq. (2) describes only phase fluctuations. By retaining leading terms in the phase gradient one can also identify J with the phase stiffness of the SC. The standard way to treat the XY model analytically, see Ref. [5] for a review, can be sketched in two main steps: 1. The mapping to the Villain model [26], and 2. the study of the BKT renormalization group equations.
The first step consists essentially in writing Eq. (2) as the sum of two decoupled contributions: an effective harmonic Hamiltonian H SW , which accounts for the longitudinal spin-waves excitations, and a vortex Hamiltonian H V , which accounts for the transverse (topological) spin excitations. The spin-wave fluctuations are responsible for the lack of long-range order at any finite temperature [41]. Their anharmonicity reduces the value of the superfluid stiffness J s (T ), which is equal to the bare coupling J at T = 0, continuously in temperature without inducing yet any phase transition.
The topological vortex excitations, on the other hand, drive the BKT phase transition, whose hallmark is the universal jump of the superfluid stiffness [25]. At the critical point (T = T BKT ), the proliferation of free vortices transforms the system from a quasi-ordered phase with zero magnetization and finite phase rigidity (J s = 0) into a disordered phase in which the system is no longer rigid (J s = 0).
As widely discussed in the literature, the Villain model H V is equivalent to the Hamiltonian of a 2D Coulomb gas, where positive and negative charges (q i = ±1) play the role of vortices and anti-vortices, respectively. Once restricted to the case i q i = 0, it reads: where ν v is the vortex-pairs density and a the lattice spacing of the discrete model. The first term in Eq. (3) describes the interaction between two charges at distance | r ij | a from each other, while the second term corresponds to the energetic cost of nucleating vortices at the shortest length scale with vortex-core energy µ v .
The effect of the transverse phase fluctuation is well captured by the BKT renormalization-group equations [2][3][4], rigorously derived for the Coulomb gas model Eq. (3) at leading order in the vortex fugacity g → 0 [39,[42][43][44]. The relevant variables are the dimensionless couplings expressed in terms of the spin stiffness J and the vortex fugacity z = e −βµv . At long distances, the values of the two couplings K and g are obtained from the numerical solution of the celebrated BKT flow equations [3,4]: where l = ln(r/a) is the rescaled length scale. The long-distance (thermodynamic) value of the superfluid stiffness J s is thus determined by J s ≡ T K(l → ∞)/π. The above equations identify two main regimes: a low-temperature regime K > 2, where the vortex fugacity flows to zero g → 0 while the superfluid stiffness flows to a constant value K → K * , and a high-temperature regime K < 2, where the vortex fugacity grows g → ∞ and the superfluid stiffness vanishes K → 0. The BKT critical temperature T BKT is defined as the temperature at which the flow reaches the fixed point K = 2, g = 0, and is equivalent to the Nelson criterion [25]. The renormalization of the stiffness from the initial value to a smaller J s = T K * /π at T < T BKT depends quantitatively on the value of the vortex-core energy µ v : the larger the initial value of the vortex fugacity g in Eqs. (6)-(7), the stronger the relative suppression of the stiffness due to bound vortex-antivortex pairs which nucleate at small length scales. The flow Eqs. (6)-(7) account only for the effect of vortex excitations, so that any other excitation that contributes to renormalize the superfluid stiffness or the vortex-core energy in temperature must be introduced by hand. This is the key idea behind the two-step approach, which consists in incorporating the effect of non-critical fluctuations into the effective couplings J eff and µ eff v , which are then used in place of the bare values J and µ v as initial conditions for the BKT flow [5,32,37,38]. While for the case of 2D superconducting films the depletion of the superfluid stiffness is mainly due to quasiparticle excitations, well accounted for by the BCS expression J eff = J BCS s (T ) [37,38], within the XY model the effective superfluid stiffness J eff (T ) is reduced with respect to J by longitudinal (spin-wave like) fluctuations already at one loop order. On the other hand, since the XY model has only a single energy scale J, it is natural to assume that also µ v scales with J eff (T ). A natural ansatz for the initial values of the couplings is then The estimate in Eq. (9) has been obtained by including the (∇θ) 4 contribution in the low energy spinwave expansion of the XY Hamiltonian (2), computed with perturbation theory around its quadratic Gaussian form [45], see [5] for a recent review. For the vortex-core energy, Kosterlitz and Thouless [3] derived an estimate of γ for a square lattice geometry, This estimate, however, was obtained within the 2D Coulomb gas model Eq. (3), and it has not yet been verified to what extent it corresponds to the actual value of the vortex-nucleation energy in the full XY model Eq. (2). Monte Carlo simulations offer an opportunity for this verification: the superfluid stiffness J s obtained numerically from the cosine model (2) can indeed be compared with the RG estimate implemented via the two-step procedure. Besides, from the numerical calculation of the vortex-pair density one can determine the vortex-core energy of the model and compare it with the Eq. (11).
The superfluid stiffness J s is defined as the second derivative of the free energy with respect to a phase twist ∆ x in a given direction, and it measures the response of the SC film to a transverse gauge field A. A constant field A in a given direction, say A = A x , is indeed equivalent to applying a total phase twist ∆ x = A x L across the sample of length L. The superfluid stiffness can be expressed in terms of the diamagnetic (J d ) and the paramagnetic (J p ) response functions where . . . stands for the thermal average. The explicit expressions of J d and J p are reported in Appendix A. Apart from the superfluid stiffness, we have also measured the vortex-pair density of the system, where ρ v is the vortex density and v(P i ) = +1(−1) if the plaquette P i is occupied by a (anti-)vortex and 0 otherwise. Once ν v (T ) and J s (T ) are determined from Monte Carlo simulations, one can extrapolate an effective vortex-core energy µ M C v (T ) = γ M C J s (T ) and compare it with Eq. (11). Indeed, assuming that the system at low T is formed by noninteracting vortex pairs, one can introduce the following BKT ansatz for the low-temperature pair density where µ eff v (T ) is the temperature dependent vortexcore energy in Eq. (10), while A is a geometrical factor [40]. We fit the MC numerical results for ν v to the form (17), see inset of Fig. 1, and find the parameters The numerical value of γ displayed in Eq. (18) is significantly smaller than the long-standing KT estimate Eq. (11). Apparently, the effect of the spinwave excitations and their relation to the presence of vortex pairs, which is not properly accounted for in the original derivation, is to lower the energetic cost to nucleate a vortex pair.
Then, we compare the superfluid stiffness J MC s (T ), obtained by Monte Carlo simulations of the full XY model, with the one obtained by solving the flow Eqs. (6)- (7) with the initial condition J eff (T ) and µ XY v (T ) given by (9),(10),(11) and (9),(10), (18). The results are reported in Fig. 1: with the KT ansatz for γ (dashed red line in Fig. 1) the renormalized J RG s is larger than J MC s at higher temperatures and leads to a larger value of the critical temperature T BKT = πJ s (T BKT )/2. On the other hand, the renormalized J RG s obtained using the vortex-core energy (18) (blue line in Fig. 1) gives a very good agreement both with the MC critical temperature and with the whole temperature dependence of the superfluid stiffness. The discrepancy between the two at high temperature T > T BKT can be easily explained in terms of finitesize effects in the MC numerical results. On the other hand, the deviation observed at low temperatures (T < T BKT ) could be given to the approximation which truncates high order spin-wave contributions to the depletion of the superfluid stiffness (9) and neglects the interplay between longitudinal and transverse excitations.

III. THE MODIFIED XY MODEL
So far, we have reviewed the standard analytical treatment of the XY model, solving the BKT renormalization-group equations via a two-step procedure that we have improved with a new estimate for µ eff v . As already stressed, however, within the classical XY model there is no way of tuning the value of µ v independently from the coupling J, as one would need when studying real systems. In particular, the study of thin SC films [37] has revealed a much smaller value of the vortex-core energy than expected from the XY model, a challenge for the study of the BKT transition in the limit of high vortex fugacity.
For this purpose, we introduce a modified version of the original XY model in which µ v can be tuned independently from the spin stiffness J. The bare vortex-core energy in Eq. (10) can be modified by adding an extra potential term to the Hamiltonian in Eq. (2), Here, I Pi corresponds to the spin current circulating around a single plaquette P i of area a 2 , The advantage of choosing an additional potential term of the form in Eqs. (20), (21) is that it is a local term, it is derivable with respect to the phase difference, and it has a direct physical interpretation in terms of currents. For µ = 0 one recovers the classical XY Hamiltonian (2) while for µ > 0 the additional potential term favours the presence of spin currents in the system. This enhances the nucleation of vortices and reduces the resulting vortex-core energy µ v . In the following, we will first discuss the Monte Carlo numerical study of the modified XY model. Then, we will analyse it theoretically using the RG equations for the BKT transition, implemented via the two-steps procedure.

IV. MONTE CARLO SIMULATIONS
In Fig. 2, the MC numerical results for the temperature dependence of J d , J p and J s are reported for different values of µ. For larger values of µ (smaller µ v ) the BKT critical temperature decreases, since at lower values of the vortex-core energy the entropic term in the free energy of a single vortex F v = E v − T S v dominates already for lower temperature and brings forward the vortexantivortex unbinding.
As argued above, the decoupling between transverse and longitudinal phase fluctuations cannot be rigorously applied to the XY model as it only holds for its quadratic approximation, the Villain model [26]. At low temperatures the curves for the diamagnetic response function J d collapse onto each other and thus demonstrate the decoupling between spin-wave and vortex excitations in the T → 0 limit. At higher temperatures the curves separate, and for larger µ there is a stronger depletion of the diamagnetic currents due to the interplay between spinwave excitations and vortices, caused by the anharmonic terms in the XY coupling, see Fig. 2(a). The paramagnetic response functions in Fig. 2(b) display a very slow increase in the low-temperature regime as they only receive contributions from high-order powers in the phase gradient [45].
The temperature dependence of the superfluid stiffness in Fig. 2(c) reveals another interesting feature for increasing values of µ. Indeed, for smaller values of the vortex-core energy µ v it becomes more and more difficult to distinguish between the BKT universal jump of J s at T = T + BKT and its rapid downturn, which starts already at a lower temperature. The 2T /π critical line for µ = 0.125 appears, in fact, to be halfway through the jump, rather than at its onset. This aspect can be very important for a correct interpretation of experimental data, as the observation of a BKT critical line halfway across the jump can be considered an indication for a low vortex-core energy within the system. At small vortex-core energies, the energy-entropy balance suggests the presence of a critical value µ = µ c , above which unbound free vortices are energetically favoured even at zero temperature. As a consequence, for µ > µ c the ground state of the system will change from the vortex-vacuum state with ρ v (T → 0) = 0 to a vortex-antivortex squarelattice crystal with ρ v (T → 0) = 1. A similar effect has been discussed for the 2D Coulomb gas in Refs. [46][47][48]. For the model (20) we find the critical value µ c 0.15 and focus henceforth on the regime µ < µ c . A preliminary study of the vortex crystal phase µ > µ c is presented in App. B, while a more complete analysis is deferred to future investigations.

V. RG STUDY OF THE MODIFIED XY MODEL
Let us proceed with the study of the generalised XY model (20) by means of the two-step RG approach. As in Sec. II, one can treat the modified XY Hamiltonian (20) in the spirit of the Villain approximation by separating the spin-wave from the topological excitations. It is easy to verify that H SW remains unchanged, while H v acquires a new term proportional to µ. As a consequence, while the lowtemperature estimate for the superfluid stiffness J eff (9) remains valid for µ = 0, one needs to identify a new renormalized value for µ eff v . We will proceed following the same route as before, i.e., by fitting the vortex-pair density ν v from MC numerical simulations of the modified XY model with Eq. (17). In this case, however, the vortex-core energy will depend on both J and µ. Moreover, as for the superfluid stiffness J eff s (9) the presence of low-temperature spin-wave fluctuations will modify the effective value of µ as well, whose temperature dependence can be considered in first approximation to be simply linear. In light of these considerations, we have used the following ansatz for the vortex-core energy at finite µ: For the numerical fits of ν v (T ), we have fixed the values of A and γ MC in Eqs. (18)- (19) and determined the best fits b 1 and b 2 (µ). As shown in Fig.3, they yield very good agreement with the Monte Carlo data for all values of µ.
In order to assess the validity of the two-step technique for the generalised XY Hamiltonian (20), one must carefully consider the two underlying approximations: (i) the assumption of decoupling between topological configurations and non-critical excitations; (ii) the BKT flow equation that is perturbative in the vortex fugacity.
The hypothesis (i) implies that critical and noncritical fluctuations occur on well-separated scales, the firsts occurring only on short-distances. The noncritical fluctuations are indeed present at all length scales and they contribute to renormalise the couplings in play even at long distances. However, to properly include these contributions, they should be inserted within the RG flow equations themselves. Something that, as we discussed, is not feasible so far. On the other hand, the limitation to small vortex fugacity can be easily circumvented by including additional terms in the perturbative expansion or by employing a self-consistent RG scheme [43].

A. The critical temperature
Both assumptions (i) and (ii) are satisfied for the Villain model, and the computation of the critical temperature from the BKT flow equations with J eff = J and µ eff v = µ v = Jπ 2 /2 perfectly reproduces high-precision MC results [49]. In the XY model case with µ = 0 one can efficiently compute the renormalization of the effective superfluid stiffness due to spin waves via functional RG, meanfield or low-temperature expansion, yielding estimates for the BKT critical temperature [32,50,51] within 5% of the numerically exact MC value [52][53][54]. Finally, the computation of the effective superfluid stiffness for 2D Fermi gases via Landau's quasiparticle-excitation formula produces a consistent picture for the dependence of the BKT temperature on the scattering length [30].
This continues to work in the modified XY model (20): as the chemical potential µ grows, topological fluctuations at increasingly smaller scales enhance the interplay between longitudinal and transverse excitations and seem to undermine the validity of both assumptions (i) and (ii). Yet, the results obtained for the BKT temperature as a function of µ are in fairly good agreement with the MC results, see Fig. 4.
Applying the finite-size scaling procedure described in Ref. [55,56] to our MC data, we obtain a reliable estimate for the critical temperature T BKT (black points in Fig. 4), which nicely reproduces the high-precision results [54,56] at µ = 0 (black star). The theoretical estimates for T BKT obtained by inserting the effective couplings into the initial conditions for the BKT flow Eqs. (6) and (7), produce a consistent picture for T BKT up to µ 0.15, where the nature of the ground state changes and vortex fluctuations become relevant at T = 0. The results from the two-step approach with the traditional RG equations are shown as red circles in Fig. 4. It is worth noting that the RG analysis performed using the effective couplings Eqs. (23)-(24) estimated using the coefficients in Tab. I greatly improves the accuracy of the two-step approach in the µ = 0 case, with respect to the traditional KT initial condition µ v = π 2 2 J eff . As the coupling µ increases, the RG results increasingly deviate from the MC estimates, in agreement with the expectations of larger corrections to the Coulomb gas approximation. At this stage, it is not clear whether these larger deviations are due to an increase in the interplay between longitudinal and transverse phase fluctuations or to higher-order vortex fugacity corrections to the traditional BKT flow equations.
Higher-order terms in the vortex fugacity can be introduced into the BKT flow equations using different approaches. As a first trial, we employed the modified RG equations described in Ref. [40] using the effective couplings discussed in Eqs. (23) and (24). This analysis slightly improves the accuracy of the predicted critical temperatures especially at high µ values, as shown with the blue diamonds in Fig. 4. In order to further address the issue of higher-order fugacity corrections to the BKT flow in Eqs. (6) and (7) in the next section we perform the two-step approach treating the Coulomb gas problem via the functional RG approach.

B. Functional renormalization-group flow for vortex unbinding
The functional RG approach (FRG) is based on the possibility to write an exact RG equation for the effective action [57][58][59][60], which may then be solved by projecting it onto a restricted theory (coupling) space with a chosen ansatz [61,62]. The main difference between the FRG scheme and the traditional approach by Wilson lies in the introduction of a momentum-dependent infrared regulator R k (p) in order to remove the divergence of the propagator close to the critical point and to allow for the deriva-tion of nonperturbative flow equations as the cutoff scale k is lowered.
In principle, the use of a nonperturbative framework should allow one to study the RG flow for the XY model beyond the quadratic (Villain) limit without the need for the two-step approach described earlier in this section. Nevertheless, efforts to construct an RG flow capable of describing the emergence of topological excitations from microscopic physics were hindered by the difficulty to reproduce the line of fixed points characteristic of the BKT transition [42,63,64] without the use of ad-hoc techniques [65] [66].
Therefore, we will perform the FRG study of the unbinding transition within the low-energy sine-Gordon model [67,68]. The flow equations for the vortex fugacity and the stiffness can be conveniently rewritten as with cutoff k = 1/r = exp(− )/a, inverse mass w k = 1/(2K k ), fugacity u k = g k /π and the inverse propagator P k = w k p 2 + R k . The flow Eqs. (25) and (26) have been derived within the derivative expansion approximation, and their nonperturbative character is apparent from the presence of arbitrary powers of the coupling u k . The flow eqs. (25) and (26) have been shown to reproduce the salient features of the BKT transition irrespectively of the choice of the regulator function [67,68], and to consistently reproduce the results of Zamolodchikov's c-theorem [69] in the sine-Gordon model and its generalisations [70,71]. Moreover, when more advanced functional truncations are considered, the FRG approach yields accurate predictions for the excitation spectrum of the sine-Gordon model [72]. For explicit computations it is necessary to specify the form of the regulator function. One straightforward choice is the power-law regulator where a and b are two free parameters, which are chosen based on an optimisation criterion. However, the common criteria to optimise the regulator function within the FRG approach do not apply to the computation of non-universal quantities [73,74]. Furthermore, the universal features of the BKT transition, such as the jump of the superfluid stiffness, are reproduced by the flow eqs. (25) and (26) 6) and (7) (dashed lines) and the FRG flow Eqs. (28) and (29), both with the initial couplings given in Eq. (23) and (24). From top left to bottom right we show the results for µ = 0, 0.025, 0.05, 0.075, 0.1, 0.125.
independently from the regulator function R k [67] and they cannot be used as a criterion for the choice of the regulator. Therefore, we are going to use a different criterion for the choice of the optimal values of the parameters a and b. In order to obtain the perturbative BKT flow in Eqs. (6) and (7), one has to assume a short-distance regularisation, which traditionally relies on considering the Coulomb gas charges as hard disks of finite radius [3]. As discussed above, the RG flow produced by such phenomenological regularisation has been shown to produce also quantitatively reliable result for small vortex fugacities, see Fig. (4). Then, we choose the parameters a and b in the regulator function (27) in such a way that the flow Eqs. (25) and (26) reproduce Eqs. (6) and (7) at leading order in the vortex fugacity u k .
Such optimisation procedure yields the optimal parameter b = 1, which leads to the analytical flow equations whereũ k = u k /k 2 and a = 1/ √ 6π 2 in order to reproduce the BKT flow at leading order. The estimate of the critical temperatures obtained by the nonperturbative flow eqs. (28) and (29) with initial conditions Eq. (23) and (24) are shown as green squares in Fig. 4. We find that the FRG approach produces better results than the traditional BKT flow or the one presented in Ref. [40] and always agrees with MC within uncertainties, apart from µ = 0.125 very close to the vortex lattice transition.

C. The superfluid stiffness
Apart from the prediction for the critical temperature T BKT , it is interesting to verify the reliability of the two-step approach in modelling different thermodynamic quantities. The most relevant quantity for our analysis is the superfluid stiffness J s (T ), which is obtained in our model by renormalising the low-temperature effective stiffness J eff in Eq. (9). The results for the superfluid stiffness obtained both by the traditional RG flow and the FRG generalisations are shown in Fig. 5. As expected from our optimisation procedure the results of the RG flow eqs. (6) and (7) are very close to the FRG ones for small vortex fugacities µ = 0, 0.025, 0.05, see first three panels in Fig. 5. As µ grows, the FRG results for the superfluid show a much better agreement with the MC data than the perturbative RG.
One crucial feature which emerges from the comparison between the FRG and the RG estimates for the superfluid stiffness with respect to the numerical data is that the functional shape of the MC J s (T ) below T BKT could not be correctly reproduced by any of these approaches. Although the FRG results show the correct unbinding point, they still slightly overestimate the value of J s below T BKT . The stronger depletion of the MC J s (T ) approaching T BKT from below appears to persist in the thermodynamic limit, while above the critical temperature the MC J s (T ) are reduced to reproduce the discontinuous jump in the thermodynamic limit. In conclusion, the pre-critical descent seems to be a direct indication of the enhanced interplay between spin-wave fluctuations and vortices at high vortex fugacities.

VI. CONCLUSIONS
The main goal of the paper is to prove the possibility to extract accurate nonuniversal quantities (J s (T ), T BKT ) for superfluid systems close to the vortex unbinding transition using the RG description for the Coulomb gas. The procedure to connect the low-energy Coulomb gas description with the microscopic model by effective bare initial conditions in the BKT flow has been widely tested on the square-lattice XY model [26,30,32,37,42,49,51], but its accuracy for large vortex fugacities had still to be investigated. In order to prove the feasibility of such a two-step procedure we employed MC simulations to study a modified version of the XY model, where the vortex-core energy can be tuned by a parameter µ in the Hamiltonian (20). The values of the superfluid stiffness (J s ), the diamagnetic (J d ) and paramagntic (J p ) currents, see Fig. 2, and the unbinding temperature T BKT have been numerically derived and compared with the theoretical predictions obtained by inserting the initial conditions given by Eqs. (4) and (5) with the effective couplings in Eqs. (23) and (24) into the BKT and FRG flows.
A summary of our analysis can be found in Fig. 4, where the MC value of T BKT is compared with the one obtained by the two-step procedure using different RG approaches. As expected, the reliability of the RG flow in Eqs. (6) and (7) decreases for lower vortex core energies µ v . This discrepancy in the prediction from the two-step approach may be partially repaired by the introduction of the nonperturbative flow Eqs. (28) and (29), which yield accurate estimates for T BKT in the whole µ range.
Conversely, the study of the superfluid stiffness reported in Fig. 5 shows that even the FRG result cannot capture the pronounced depletion in the superfluid stiffness occurring below T BKT , which is most probably the result of the large interplay between transverse e longitudinal fluctuations occurring at large µ. This interplay cannot be captured by the two-step approach, and hence this feature remains unmatched in the RG analysis. Surprisingly, these missing features do not preclude accurate estimates of the critical temperatures via the FRG approach.
In conclusion, our analysis proves that the BKT flow equations Eqs. (6) and (7) represent more than just a low-energy description of the universal behaviour of 2D superfluid systems and may more generally be used to construct a quantitative theory for the unbinding transition once a suitable value for the bare superfluid stiffness J eff has been identified. For this purpose, we employed the numerical simulations for the vortex density in the system to determine the effective vortex core energy via the low-temperature ansatz in Eq. (17), which, once inserted into the RG flow equations, provides a consistent picture for vortex unbinding.
Acknoledgements. We thank J. We have performed MC simulations of the Hamiltonian in Eq. (20). The linear system sizes considered are L = 8, 16, 32, 64, 128, 256 with periodic boundary conditions. For each value of µ ∈ {0, 0.025, 0.05, 0.075, 0.1, 0.125} and size, we have run 10 4 MC steps and discarded the transient regime occurring in the first 2 × 10 3 steps. Each MC step consists of five canonical Metropolis spin flips of the whole lattice, followed by ten micro-canonical overrelaxation sweeps of all spins. The thermalization at low temperatures has been sped up by a temperature annealing procedure. Finally, the observables measured have been averaged both over the canonical ensemble (thermal average) and over five independent samples.
The explicit expressions for the diamagnetic (J d ) and the paramagnetic (J p ) current (alongx) for the model Eq. (20) studied are For each size L of the system the Nelson criterion [25] J s /T BKT = 2 π has been applied to calculate the finite-size unbinding transition temperature T BKT (L). Then, for each value of µ, the thermodynamic critical temperature T BKT (∞) is extrapolated by means of the finite-size scaling analysis based on the behaviour of the BKT correlation function close to criticality: ξ A exp(c/ √ t), with t the reduced temperature t = (T − T BKT )/T BKT . The relation used to fit our data and estrapolate T BKT is [56] with b and c fitting paramenters. The trend of β BKT (L) as a function of (ln(bL)) −2 is shown for each µ and for the best fitted parameters in Fig. 6. In Sec. IV of the main text, it has already been explained that we expect a transition at µ ≈ 0.15 between two different zero-temperature ground states: the vortex vacuum with ρ v (T = 0) = 0 for µ < µ c , and a vortex crystal with ρ v (T = 0) = 1 for µ > µ c . Even if the study of the vortex crystal phase was not the main scope of the present paper, we would like to present here some results at µ = 0.15, 0.2, which have been excluded from the analysis in the main text as they belong to a different critical scenario.
The numerical analysis of the melting transition at µ 0.15 is reported on Fig. 7. The superfluid stiffness in Fig. 7(a) displays qualitatively the same behaviour as in the vortex-unbinding regime. The presence of a different ground state is reflected in the zero-temperature value of the superfluid stiffness which increases with µ, as expected since vortex fluctuation are in this case relevant also at low temperatures. On the other hand, the superfluid stiffness still presents a sharp jump at the melting point of the 2D crystal. Even if the melting transition in temperature is expected in two dimensions to belong to the BKT universality class, the present data do not have enough precision to investigate whether this jump is consistent with a BKT transition or rather has the characteristic of a continuous transition.
Nonetheless, more insight can be gained by inspecting the behaviour of the vortex density for the two cases µ = 0.15, 0.2 shown in Fig. 7(b). The two curves display a rather different behaviour: in the µ = 0.15 case one has ρ v = 1 at T 0, indicating the crystal phase; as T increases the vortex density drops (almost) discontinuously to ρ v < 0.5 and, then, slowly increases toward a high-temperature value in agreement with the one observed for the disordered phase at µ = 0.125. Conversely, the vortex density at µ = 0.2 decreases monotonically (and more smoothly) as a function of the temperature. Partial justification for this behaviour can be found by comparing this behaviour with the one of the lattice Coulomb gas studied in Ref. [47]. There, the vortex crystal at low temperatures appears for fugacities g 0.4 and melts either with a first order transitions for g ≈ 0.4 or with a second order phase transition. Our results suggest that the case µ = 0.15 belongs to the first scenario, while µ = 0.2 to the second one.
The validity of the Coulomb gas description for this model even in the crystal phase is supported by the behaviour of the critical temperature as a function of µ. In the superfluid region T BKT decreases with increasing µ until one encounters the critical threshold µ c ≈ 0.15, above which the critical tem-perature of the system starts increasing with µ. The resulting curve is nonmonotonic but perfectly continuous, as in the neutral Coulomb gas case. This behaviour is confirmed by Fig. 7(c), where the critical temperatures for the superfluid and the crystal melting transition are reported vs. µ. Further numerical analysis would be necessary to clarify possible differences arising in the modified XY model, due to the coupling of vortices and spin waves, with respect to the Coulomb gas case [47].