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

More than 40 years after its observation in thin ⁴He films Bishop and Reppy (1978), the Berezinskii-Kosterlitz-Thouless (BKT) transition Berezinsky (1972); Kosterlitz and Thouless (1973); Kosterlitz (1974); José (2013) 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 Corson et al. (1999); Hetel et al. (2007); Bilbro et al. (2011); Yong et al. (2012); Baity et al. (2016), exciton-polariton systems Nitsche et al. (2014), cold atoms in 2D harmonic traps Hadzibabic et al. (2006); Murthy et al. (2015) and 2D electron gases Reyren et al. (2007); Daptary et al. (2016); Monteiro et al. (2017) 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) Lemarié et al. (2013); Sacépé et al. (2011), raising the question of its effect on the BKT phase transition Maccari et al. (2017, 2018); Venditti et al. (2019).

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 Nisoli and Bishop (2014); Seul et al. (1991); Mendoza-Coto and Stariolo (2012). 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 relationNelson and Kosterlitz (1977), the value $J_{s}(T_{BKT})$ is a universal function of the critical temperature itself, i.e.,

$$J_{s}(T_{BKT})=\frac{2T_{BKT}}{\pi}.$$

(1)

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 Kosterlitz (1974); Nelson and Kosterlitz (1977), 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, to quantify this effect one needs a precise knowledge of the vortex-core energy $\mu_{v}$, which controls the vortex fugacity $g=\exp(-\mu_{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 nontopological 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 Berezinsky (1972); Kosterlitz and Thouless (1973); Kosterlitz (1974). 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 temperature 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 renormalized initial conditions for the BKT flow equations. Since the work of J. Villain within the context of 2D classical planar magnets Villain (1975), 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 Lieb (1967). To account for the nontopological phase modes of the system, in addition to the RG approach, different theoretical techniques have been proposed, such as the self-consistent harmonic approximation Villain (1975); Pires (1997), classical Monte Carlo simulations Prokof’ev et al. (2001), Gaussian fluctuations Bighin and Salasnich (2016); Karle et al. (2019) and the functional renormalization group Defenu et al. (2017); Krieg and Kopietz (2017).

A second example of the relevance of nonuniversal effects is the application to thin films of superconductors. In this case, it has been suggested Kamlapure et al. (2010); Mondal et al. (2011a); Benfatto et al. (2007); Mondal et al. (2011b); Yong et al. (2013) 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 $\mu_{v}$ is significantly smaller than expected within the XY model. Indeed, within a microscopic BCS picture for superconductors $\mu_{v}$ is expected not to scale with the superfluid stiffness, but rather with the Cooper-pair condensation energy $E_{c}$. However, within the XY model, the ratio $\mu_{v}/J_{s}$ has been estimated Kosterlitz and Thouless (1973) long ago to be a constant $\mu_{v}/J_{s}=\pi^{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 larger 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 $\mu_{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 Eq. (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 Amit et al. (1980), via the self-consistent approach suggested by Timm Timm (1996) 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 vortex-core energy $\mu_{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 $\mu_{v}^{\text{eff}}$ and the superfluid stiffness $J_{\mathrm{eff}}$ are inserted into the BKT flow equations. Finally, we conclude in Sec. V.

Ever since the seminal papers by Berezinskii, Kosterlitz and Thouless Berezinsky (1972); Kosterlitz and Thouless (1973); Kosterlitz (1974), 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 ($|\vec{S_{i}}|=1$), interacting via a ferromagnetic spin-spin coupling constant $J>0$:

	$\displaystyle H_{XY}$	$\displaystyle=-J\sum_{i,\nu=\hat{x},\hat{y}}\vec{S}_{i}\cdot\vec{S}_{i+\nu}$
		$\displaystyle=-J\sum_{i,\nu=\hat{x},\hat{y}}\cos(\theta_{i}-\theta_{i+\nu}),$		(2)

where $\theta_{i}$ is the angle the $i$th spin forms with a given direction. Within the context of SC films, one can map $\vec{S}_{i}\to|\Delta|e^{i\theta_{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. José (2013) for a review, can be sketched in two main steps:

1. 1.

   The mapping to the Villain model Villain (1975), and

2. 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 Mermin and Wagner (1966). 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, however, drive the BKT phase transition, whose hallmark is the universal jump of the superfluid stiffness Nelson and Kosterlitz (1977). 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}\neq 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}=\pm 1$) play the role of vortices and anti-vortices, respectively. Once restricted to the case $\sum_{i}q_{i}=0$, it reads:

$$H_{V}=-\pi J\sum_{i\neq j}\ln\left(\frac{|\vec{r_{ij}}|}{a}\right)q_{i}q_{j}+2\mu_{v}\sum_{i}\,\nu_{v}({\bf r}_{i})\ ,$$

(3)

where $\nu_{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 $|\vec{r}_{ij}|\gg 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 $\mu_{v}$.

The effect of the transverse phase fluctuation is well captured by the BKT renormalization-group equations Berezinsky (1972); Kosterlitz and Thouless (1973); Kosterlitz (1974), rigorously derived for the Coulomb gas model Eq. (3) at leading order in the vortex fugacity $g\to 0$ José et al. (1977); Amit et al. (1980); Minnhagen (1987); Gulacsi and Gulacsi (1998). The relevant variables are the dimensionless couplings

	$\displaystyle K(0)$	$\displaystyle=$	$\displaystyle\frac{\pi J}{T},$		(4)
	$\displaystyle g(0)$	$\displaystyle=$	$\displaystyle 2\pi e^{-\beta\mu_{v}},$		(5)

expressed in terms of the spin stiffness $J$ and the vortex fugacity $g$.

At long distances, the values of the two couplings $K$ and $g$ are obtained from the numerical solution of the celebrated BKT flow equations Kosterlitz and Thouless (1973); Kosterlitz (1974):

	$\displaystyle\frac{dK}{dl}$	$\displaystyle=$	$\displaystyle-K^{2}g^{2},$		(6)
	$\displaystyle\frac{dg}{dl}$	$\displaystyle=$	$\displaystyle(2-K)g,$		(7)

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}\equiv TK(l\to\infty)/\pi$. The above equations identify two main regimes: a low-temperature regime $K>2$, where the vortex fugacity flows to zero $g\to 0$ while the superfluid stiffness flows to a constant value $K\to K^{*}$, and a high-temperature regime $K<2$, where the vortex fugacity grows $g\to\infty$ and the superfluid stiffness vanishes $K\to 0$. The BKT critical temperature $T_{BKT}$ is defined as the temperature at which the flow reaches the fixed point $K=2$, $g=0$,

$$K(l\to\infty,T_{BKT})=2\implies\frac{\pi J_{s}(T_{BKT})}{T_{BKT}}=2,$$

(8)

and is equivalent to the Nelson criterion Nelson and Kosterlitz (1977). The renormalization of the stiffness from the initial value to a smaller $J_{s}=TK^{*}/\pi$ at $T<T_{BKT}$ depends quantitatively on the value of the vortex-core energy $\mu_{v}$: the larger the initial value of the vortex fugacity $g$ in Eqs. (6) and (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) and (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 noncritical fluctuations into the effective couplings $J_{\mathrm{eff}}$ and $\mu^{\mathrm{eff}}_{v}$, which are then used in place of the bare values $J$ and $\mu_{v}$ as initial conditions for the BKT flow José (2013); Mondal et al. (2011b); Yong et al. (2013); Defenu et al. (2017). 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_{\mathrm{eff}}=J_{s}^{\text{BCS}}(T)$ Mondal et al. (2011b); Yong et al. (2013), within the XY model the effective superfluid stiffness $J_{\mathrm{eff}}(T)$ is reduced with respect to $J$ by longitudinal (spin-wave like) fluctuations already at one loop order. However, since the XY model has only a single energy scale $J$, it is natural to assume that also $\mu_{v}$ scales with $J_{\mathrm{eff}}(T)$. A natural ansatz for the initial values of the couplings is then

	$\displaystyle J_{\mathrm{eff}}(T)$	$\displaystyle=J-\frac{T}{4}\,,$			(9)
	$\displaystyle\mu_{v}^{\mathrm{eff}}(T)$	$\displaystyle=\gamma J_{\mathrm{eff}}(T).$			(10)

The estimate in Eq. (9) has been obtained by including the $(\nabla\theta)^{4}$ contribution in the low energy spin-wave expansion of the XY Hamiltonian (2), computed with perturbation theory around its quadratic Gaussian form Maccari et al. (2019), see Ref. José (2013) for a recent review. For the vortex-core energy, Kosterlitz and Thouless Kosterlitz and Thouless (1973) derived an estimate of $\gamma$ for a square lattice geometry,

$$\gamma^{KT}=\frac{\pi^{2}}{2}.$$

(11)

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 Eq. (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 $\Delta_{x}$ in a given direction,

$$J_{s}^{x}\equiv-\frac{\partial^{2}F(\Delta_{x})}{\partial\Delta_{x}^{2}}\Big{|}_{\Delta_{x}=0},$$

(12)

and it measures the response of the SC film to a transverse gauge field $\bf{A}$. A constant field $\bf{A}$ in a given direction, say $\mathbf{A}=A_{x}$, is indeed equivalent to applying a total phase twist $\Delta_{x}=A_{x}L$ across the sample of length $L$. The superfluid stiffness

$\displaystyle J_{s}^{x}=J_{d}^{x}-J_{p}^{x}$

(13)

can be expressed in terms of the diamagnetic ($J_{d}$) and the paramagnetic ($J_{p}$) response functions

	$\displaystyle J_{d}^{x}$	$\displaystyle=\frac{1}{L^{2}}\Big{[}\Bigl{\langle}\frac{\partial^{2}H}{\partial A_{x}^{2}}\Big{|}_{0}\Bigr{\rangle}\Big{]},$		(14)
	$\displaystyle J_{p}^{x}$	$\displaystyle=\frac{\beta}{L^{2}}\Big{[}\Bigl{\langle}\Bigl{(}\frac{\partial H}{\partial A_{x}}\Big{|}_{0}\Bigr{)}^{2}\Bigr{\rangle}-\Bigl{\langle}\frac{\partial H}{\partial A_{x}}\Big{|}_{0}\Bigr{\rangle}^{2}\Big{]},$		(15)

where $\langle\dots\rangle$ 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,

$$\nu_{v}=\frac{1}{2}\rho_{v}=\frac{1}{2}\biggl{\langle}\frac{1}{L^{2}}\sum_{i}^{N}|v(P_{i})|\biggr{\rangle},$$

(16)

where $\rho_{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 $\nu_{v}(T)$ and $J_{s}(T)$ are determined from Monte Carlo simulations, one can extrapolate an effective vortex-core energy $\mu_{v}^{MC}(T)=\gamma^{MC}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

$${\nu_{v}(T)=Ae^{-2\beta\mu^{\mathrm{eff}}_{v}(T)}},$$

(17)

where $\mu_{v}^{\mathrm{eff}}(T)$ is the temperature dependent vortex-core energy in Eq. (10), while $A$ is a geometrical factor Timm (1996). We fit the MC numerical results for $\nu_{v}$ to the form Eq. (17), see inset of Fig. 1, and find the parameters

	$\displaystyle\gamma^{\mathrm{MC}}=3.56\pm 0.01,$		(18)
	$\displaystyle A^{\mathrm{MC}}=1.87\pm 0.06.$		(19)

The numerical value of $\gamma$ displayed in Eq. (18) is significantly smaller than the long-standing KT estimate Eq. (11). Apparently, the effect of the spin-wave 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^{\mathrm{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_{\mathrm{eff}}(T)$ and $\mu_{v}^{XY}(T)$ given by Eqs . (9),(10),(11) and (18). The results are reported in Fig. 1: with the KT ansatz for $\gamma$ (dashed red line in Fig. 1) the renormalized $J_{s}^{\mathrm{RG}}$ is larger than $J_{s}^{\mathrm{MC}}$ at higher temperatures and leads to a larger value of the critical temperature $T_{BKT}=\pi J_{s}(T_{BKT})/2$. On the other hand, the renormalized $J_{s}^{\mathrm{RG}}$ 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 finite-size effects in the MC numerical results. However, 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 Eq. (9) and neglects the interplay between longitudinal and transverse excitations.

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 $\mu_{v}^{\mathrm{eff}}$. As already stressed, however, within the classical XY model there is no way of tuning the value of $\mu_{v}$ independently from the coupling $J$, as one would need when studying real systems. In particular, the study of thin SC films Mondal et al. (2011b) 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 $\mu_{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),

$\displaystyle H_{XY}^{{\mu}}=-J\sum_{i,\nu=\hat{x},\hat{y}}\cos(\theta_{i}-\theta_{i+\nu})-{\mu}\sum_{i}\big{(}I_{P_{i}}\big{)}^{2}.$

(20)

Here, $I_{P_{i}}$ corresponds to the spin current circulating around a single plaquette $P_{i}$ of area $a^{2}$,

$$\begin{split}I_{P_{i}}=&\sin(\theta_{i}-\theta_{i+\hat{x}})+\sin(\theta_{i+\hat{x}}-\theta_{i+\hat{x}+\hat{y}})\\ +&\sin(\theta_{i+\hat{x}+\hat{y}}-\theta_{i+\hat{y}})+\sin(\theta_{i+\hat{y}}-\theta_{i}).\end{split}$$

(21)

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 $\mu=0$ one recovers the classical XY Hamiltonian (2) while for ${\mu}>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 $\mu_{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.

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 $\mu$. For larger values of $\mu$ (smaller $\mu_{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}-TS_{v}$ dominates already for lower temperature and brings forward the vortex-antivortex 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 Villain (1975). 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\to 0$ limit. At higher temperatures the curves separate, and for larger $\mu$ there is a stronger depletion of the diamagnetic currents due to the interplay between spin-wave 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 Maccari et al. (2019).

The temperature dependence of the superfluid stiffness in Fig. 2(c) reveals another interesting feature for increasing values of $\mu$. Indeed, for smaller values of the vortex-core energy $\mu_{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/\pi$ critical line for $\mu=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 $\mu=\mu^{c}$, above which unbound free vortices are energetically favoured even at zero temperature. As a consequence, for $\mu>\mu^{c}$ the ground state of the system will change from the vortex-vacuum state with $\rho_{v}(T\to 0)=0$ to a vortex-antivortex square-lattice crystal with $\rho_{v}(T\to 0)=1$. A similar effect has been discussed for the 2D Coulomb gas in Refs. Lee and Teitel (1990, 1992); Lidmar and Wallin (1997). For the model (20) we find the critical value $\mu^{c}\simeq 0.15$ and focus henceforth on the regime $\mu<\mu^{c}$. A preliminary study of the vortex crystal phase $\mu>\mu_{c}$ is presented in App. B, while a more complete analysis is deferred to future investigations.

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 $\mu$. As a consequence, while the low-temperature estimate for the superfluid stiffness $J_{\mathrm{eff}}$ (9) remains valid for $\mu\neq 0$, one needs to identify a new renormalized value for $\mu^{\mathrm{eff}}_{v}$.

We will proceed following the same route as before, i.e., by fitting the vortex-pair density $\nu_{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 $\mu$. Moreover, as for the superfluid stiffness $J^{\text{eff}}_{s}$ (9) the presence of low-temperature spin-wave fluctuations will modify the effective value of $\mu$ 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 $\mu$:

$$\mu_{v}^{\text{eff}}(T,\mu)=\gamma^{\mathrm{MC}}J_{\mathrm{eff}}(T)-b_{1}\mu+b_{2}(\mu)T.$$

(22)

For the numerical fits of $\nu_{v}(T)$, we have fixed the values of $A$ and $\gamma^{\mathrm{MC}}$ in Eqs. (18)-(19) and determined the best fits $b_{1}$ and $b_{2}(\mu)$. As shown in Fig.3, they yield very good agreement with the Monte Carlo data for all values of $\mu$.

To assess the validity of the two-step technique for the generalised XY Hamiltonian Eq. (20), one must carefully consider the two underlying approximations:

1. (i)

   the assumption of decoupling between topological configurations and noncritical excitations;

2. (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 renormalize 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. However, 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 Minnhagen (1987).

The main goal of the paper is to prove the possibility to extract accurate nonuniversal quantities ($J_{s}(T)$, $T_{\text{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 Villain (1975); José et al. (1977); Defenu et al. (2017); Benfatto et al. (2013); Mondal et al. (2011b); Bighin and Salasnich (2016); Janke and Nather (1993), 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 $\mu$ 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_{\mathrm{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_{\mathrm{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 $\mu_{v}$. This discrepancy in the prediction from the two-step approach may be partially repaired by the introduction of the non-perturbative flow Eqs. (28) and (29), which yield accurate estimates for $T_{\mathrm{BKT}}$ in the whole $\mu$ 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_{\mathrm{BKT}}$, which is most probably the result of the large interplay between transverse and longitudinal fluctuations occurring at large $\mu$. 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_{\mathrm{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. Lorenzana, N. Dupuis, G. Gori, and A. Trombettoni for stimulating discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via Project-ID 273811115 - SFB1225 (ISOQUANT) and under Germany’s Excellence Strategy “EXC-2181/1-390900948” (the Heidelberg STRUCTURES Excellence Cluster), by the Italian MAECI under the Italian-India collaborative project SUPERTOP-PGR04879, by the Italian MIUR project PRIN 2017 No. 2017Z8TS5B, by Regione Lazio (L. R. 13/08) under project SIMAP, by Sapienza University under project Ateneo 2019 (Grant No. RM11916B56802AFE) and by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine.