Vortex supersolid in the XY model with tunable vortex fugacity

In this paper, we investigate the XY model in the presence of an additional potential term that independently tunes the vortex fugacity favouring their nucleation. By increasing the strength of this term and thereby the vortex chemical potential µ, we observe significant changes in the phase diagram with the emergence of a normal vortex-antivortex lattice as well as a superconducting vortex-antivortex crystal (lattice supersolid) phase. We examine the transition lines between these two phases and the conventional non-crystalline one as a function of both the temperature and the chemical potential. Our findings suggest the possibility of a peculiar tricritical point where second-order, first-order, and infinite-order transition lines meet. We discuss the differences between the present phase diagram and previous results for two-dimensional Coulomb gas models. Our study provides important insights into the behaviour of the modified XY model and opens up new possibilities for investigating the underlying physics of unconventional phase transitions.


INTRODUCTION
Systems displaying multiple forms of long-range order in their ground state have always fascinated physicists for their potential to exhibit a complex phase diagram.Different from simpler systems, they can host multiple phase transitions and reveal new intermediate phases between the ground state and the high-temperature phase.Apart from multi-component systems, such as multiband superconductors or bosonic mixtures, also single-component systems can present a similar scenario.The two-dimensional (2D) Coulomb gas (CG) model is a paradigmatic example.
The 2D CG is an effective model for superconducting (SC) and superfluid vortices which, in two dimensions, are equivalent to logarithmically interacting charges.In the limit of small vortex fugacity, the model undergoes a Berezinskii-Kosterlitz-Thouless (BKT) [1][2][3] transition separating a low-temperature phase, where vortices and antivortices are tightly bound in pairs, from a high-temperature phase where free vortices proliferate and lead to a discontinuous vanishing of the condensate phase rigidity.
As the vortex fugacity g increases above a critical value g c , however, the low-temperature phase of the system undergoes a first-order phase transition from a vortex-vacuum superfluid to a vortex-antivortex superfluid crystal, which additionally breaks the discrete Z 2 symmetry associated with the two energetically equivalent checkerboard configurations of the lattice.As a result, in this regime the ground state exhibits two coexisting orders: a quasi-longrange order of the superfluid phase, characterized by a finite superfluid stiffness J s , and a long-range positional order, characterized by a finite Ising order parameter for the staggered vorticity M stag .Establishing how such a vortex supersolid melts into the disordered high-temperature phase has been a topic of great interest.The phase diagram of the 2D Coulomb gas at large vortex fugacity has been extensively investigated both for discrete lattice models [4][5][6] and in the continuum limit [7].In the presence of a discrete underlying grid, it was shown [4][5][6] that at large vortex fugacity, the system undergoes two distinct phase transitions with an intermediate non-superfluid phase where the discrete Z 2 symmetry is spontaneously broken.
Addressing this problem within a 2D XY model has proven to be much more challenging.A ground state formed by a Z 2 vortex supersolid can be realized, in this model, by applying a uniform transverse magnetic field to the system with half a magnetic flux quantum crossing each plaquette of the spin lattice.The resulting model is the well-known fully frustrated XY (FFXY) model.Over the years this has been the subject of extended theoretical discussions, with a series of conflicting analytical and numerical results about the number of phase transitions and their nature [8].Finally, in 1996 Olsson [9] numerically demonstrated the presence of two phase transitions that are very close together, with the BKT critical temperature, T BKT , slightly smaller than the Ising critical temperature, T I , associated with the vanishing of M stag .The theoretical argument for the observed splitting was afterwards provided by Korshunov [10].The continuous nature of the Ising transition ensures that, when approaching T I from below, the proliferation of Ising domain walls with a net polarization continuously decreases both M stag and J s .Hence, there are in general two possible scenarios that describe the melting of a ground state with coexisting superfluidity and staggered vortex structures: 1) the system exhibits a preemptive first-order phase transition with J s and M stag vanishing discontinuously at the same critical temperature; 2) the system undergoes two phase transitions with T BKT < T I .Indeed, as soon as domain-wall excitations reduce J s below the BKT critical value J s (T BKT ) = 2T BKT /π, vortex-antivortex pairs unbind and J s drops discontinuously to zero.The FFXY model exhibits the second scenario, as confirmed also by more recent numerical studies [11,12].Yet, although the ground state of the FFXY model shares the same orders and symmetries as that of the 2D CG model at large vortex fugacity g, neither the FFXY nor the classical XY model allows for a systematic study of the phase diagram as a function of g.The XY model is, indeed, a single-coupling model where the value of the vortex fugacity cannot be tuned independently but is rather fixed by the value of the spin-exchange coupling J.
In the present work, we face this challenge by studying the phase diagram of the modified XY model that we introduced in a previous work [13], where the vortex fugacity can be tuned independently and in a direct way without changing the relevant interactions at play [14].By employing large-scale Monte Carlo simulations we assess the phase diagram of the model and show that the system undergoes a single first-order phase transition with T BKT = T I for a finite range of values of the vortex fugacity g c < g < g * , while for g > g * the two phase transitions split apart with T BKT < T I .The quantitative numerical characterisation of a BKT transition at large but finite vortex fugacity, which goes beyond the traditional BKT picture with a line of fixed points at zero fugacity, is relevant in numerous physical systems, including two-dimensional Kondo lattices [15,16], and recently in the description of the metal-insulator transition in disordered 2D materials [17].In thin superconducting films, a finite density of vortex-antivortex pairs can be induced at low temperatures by spatially correlated-disorder [18,19], while stable configurations of vortex supersolids can be realized via magnetic pinning arrays [20,21] or superconductor/ferromagnet hybrid structures [22].The formation and melting of a vortex-antivortex lattice in superfluid 4 He films can be observed by the presence of a transverse mode that can exist only in the crystalline phase, and the vortex fugacity can be tuned by additional 3 He atoms [23].More recent realisations include ultracold fermionic gases [24] and polariton fluids [25].High vortex fugacities may also emerge in long-range interacting systems.Indeed, generic power-law couplings 1/r α may disrupt the BKT in d = 2 by increasing the vortex fugacity [26,27].It is worth noting that 1/r 2 interactions induce BKT scaling also in several d = 1 models [28,29].

THE MODEL
The model studied in this work is a modified version of the original XY model with an extra potential term added to tune the vortex fugacity independently from the ferromagnetic coupling J.The Hamiltonian of the modified XY model, introduced in our previous work [13], reads: with I Pi the spin current circulating around a unit plaquette P i of area a 2 = 1, For µ = 0, Eq.( 1) is the classical XY model, where the value of the vortex fugacity is fixed by the bare spin stiffness J. On the other hand, by considering nonzero values of µ one can independently tune g to either favour for µ > 0, or disfavour for µ < 0, the vortex nucleation in the system.Thus, by increasing µ > 0, the value of the vortex-core energy µ v ∝ −µ decreases and, in turn, the value of the vortex fugacity g = 2πe −βµv increases.The energy-entropy balance for the proliferation of free vortices suggests that the BKT critical temperature decreases as the value of µ increases.At the same time, it is also apparent that there exists a critical value µ = µ c at which the ground state of the system undergoes a first-order phase transition from a superfluid with vanishing vortex density ρ v (T → 0) → 0 ("vortex vacuum") to a vortex-antivortex superfluid crystal with ρ v (T → 0) → 1 [5].
While in our previous work [13] we focused on the regime µ < µ c here we will investigate the phase diagram of the model (1) for µ > µ c .In this regime, the ground state is a checkerboard configuration of vortices and antivortices, forming a squared lattice.Being it superfluid, the ground state is a "supersolid", or more precisely a "lattice supersolid", as the presence of an underlying square grid reduces the translational symmetry broken by the crystal from a continuous to a discrete Z 2 symmetry.[30,31].As a function of µ, we will determine the value of the two critical temperatures: T BKT , at which a superfluid quasi-condensate forms, and T I , at which a charge-ordered state forms, that is described by a real Z 2 order parameter associated with the two possible staggered magnetizations of the vortex-antivortex lattice.This systematic investigation will enable us to assess the phase diagram of the system and to establish, for each value of µ, whether the system displays two separate phase transitions, or a single preemptive first-order phase transition where both the superfluid stiffness J s and the staggered magnetization M stag jump discontinuously to zero at the same critical temperature T BKT = T I .

MONTE CARLO SIMULATIONS
We assess the phase diagram of the model (1) in the regime µ > µ c via large-scale Monte Carlo (MC) simulations.This allows us to properly account for the non-trivial interactions between the different topological phase excitations at play, which include vortices, Ising-like domain walls between the two possible values of M stag , and kink-antikink excitations along the domain walls [32].
We studied the model (1) on a discrete square grid of spacing a = 1 and size N = L × L, for different values of the linear size L. Details of our MC simulations can be found in the Supplementary Materials.
To assess the values of the BKT critical temperature, we computed the superfluid stiffness J ν s , which measures the response of the system to a phase twist ∆ ν along a given direction ν.This can be thought of in terms of twisted boundary conditions, θ i+Lν = θ i + ∆ ν , reabsorbed via a gauge transformation in a new set of variables θ i = θ i − r i,ν ∆ ν /L, with periodic boundary conditions.For a superconducting film, it corresponds to the response to a transverse gauge field A and it signals the onset of perfect diamagnetism, i.e., the well-known Meissner effect.J s is defined as: and has two contributions the diamagnetic (J ν d ) and the paramagnetic (J ν p ) response functions where . . .stands for the thermal average over the MC steps.The explicit expressions of J ν d and J ν p are reported in the Appendix of [13].In this work, we have computed the superfluid response along ν ≡ x and in what follows we will simply refer to J s ≡ J x s .When increasing the temperature below T BKT , the superfluid stiffness continuously decreases mainly due to the presence of non-topological phase excitations, such as spin waves and domain walls with a net polarization [10].As soon as T BKT is reached, the proliferation of free vortices becomes entropically favoured and J s discontinuously jumps to zero.According to the Nelson-Kosterlitz criterion [33], at the critical point J s and T BKT are linked via the universal relation: J s (T BKT ) = 2T BKT /π, which ultimately allows for the determination of the critical temperature.
In this work, we assess the value of T BKT by the BKT finite-size scaling of the superfluid stiffness [34]: where L 0 is chosen to give the best crossing point at finite temperature (see also Supplementary Materials S2).The BKT finite-size scaling of J s for µ = 0.3 > µ c is reported in Fig. 1(a), where we found L 0 = 10.5.
On the other hand, in order to assess the Z 2 Ising critical temperature T I associated with the melting of the vortex-antivortex crystal, we define a vortex ordering parameter as the staggered magnetization: where i labels the unitary plaquette of the spin lattice located at (x i , y i ).The vortex charge q i is obtained by computing the phase circulation around each unitary plaquette, being:  10).The crossing point locates T I , indicated here with a dashed-dotted line.(c) Olsson's plot [9] for different values of the system size L. At the BKT critical point, while the superfluid stiffness jumps from J s (T − BKT ) = 2T − BKT /π to J s (T + BKT ) = 0, the staggered magnetization is observed to increase with L and reaches a finite value in the thermodynamic limit.This is an additional confirmation that T BKT < T I in this case.The error bars are computed via a standard bootstrapping resampling method.Where not visible, the error bars are smaller than the point symbols.
where the phase difference along each bond is defined so as to lie between the interval (−π, π].The vortex charge takes the values q i = 0, +1, −1, respectively, if a vortex, an antivortex, or zero vortices are located at the centre of the i-th plaquette.A vortex-antivortex crystal is characterised by M stag = ±1, according to the two possible equivalent configurations of the vortex-antivortex checkerboard.To determine the value of T I , we analyse the finite-size scaling of the Binder cumulant U stag associated with the staggered magnetization: In the high-temperature limit the Binder cumulant approaches U stag (T T I ) → 1 and in the low-temperature limit U stag (T T I ) → 0.3, while at the critical point it is expected to assume a universal value independent on the system size [35].The finite-size scaling of the Binder cumulant is reported in Fig. 1(b) for µ = 0.3.
At this value of the vortex chemical potential µ = 0.3, we found two distinct and yet very close critical temperatures with T BKT = 2.0040 ± 0.0003 slightly smaller than T I = 2.01595 ± 0.00004.
As a further numerical confirmation of the splitting between the two phase transitions, we follow the scheme proposed by Olsson [9].Olsson's scheme consists in extracting a set of temperatures T L for different system sizes L, which are defined as the temperatures where the superfluid stiffness crosses the 2T /π BKT critical line, i.e., J s (T L , L) = 2T L /π.By increasing the size L, T L decreases and approaches the thermodynamic limit T L→∞ → T BKT from above.If the two phase transitions are separated with T BKT < T I , the value of the staggered magnetization M stag (T L , L) at T L should increase with increasing system size L and eventually reach a nonzero value in the thermodynamic limit.This is precisely what we observe in this case, as reported in Fig. 1(c).At the temperatures T L , indicated by a dashed vertical line, the value of M stag (T L , L) increases, confirming that T BKT < T I .To establish the phase diagram of the model (1), we repeated the same analysis for different values of µ.
When approaching the critical value µ c below which the ground state is a vortex-vacuum superfluid, we find that the separation between the two phase transitions reduces until they eventually merge into a single first-order phase transition at µ = µ * > µ c .In particular, while down to µ = 0.2 (see Figs. S2-S5 of the Supplementary Materials) we still find evidence of a splitting between the two transitions, at µ = 0.175 our numerical simulations suggest that the system undergoes a single first-order phase transition.
The numerical evidence for a single first-order transition is threefold.The first indications in this sense are the failure of the BKT scaling Eq. ( 7) for the superfluid stiffness (see Fig. S6(a)) and the pronounced peaks in the Binder cumulant in the proximity of the critical point (see Fig. S6(b)) [36].
Second, an unambiguous demonstration of first-order phase transition at µ = 0.175 is provided by the presence of two peaks in the energy-density distribution P (E/N ) at the critical point.As reported in Fig. 2, at µ = 0.175 the minimum value P (E min /N ) of the distribution between the two peaks vanishes by increasing the system size L (see Fig. 2  Third, for a more quantitative analysis of the order of the transition, we looked at the finite-size scaling of the maximum value C max v of the specific heat at the critical temperature.The specific heat C v being defined as: where E is the total energy of the system.For a second-order phase transition, , where d = 2 is the spatial dimension of the system and ν = 1 is the critical exponent.Conversely, when the transition is of first order, for the Ising model in two dimensions the specific-heat peak scales as the volume of the system [36], i.e., C max v ∝ L d .For µ = 0.175, 0.2, 0.3, we have extracted the value of C max v at different system sizes L and derived the exponent C max v ∝ L y via a linear fit of the data in a log-log plot (see Fig. 3).For µ = 0.3, this analysis yields y = 0.2 ± 0.01 (see Fig. 3(a)), in good agreement with the value y = 0 expected in two spatial dimensions for a Isinglike second-order phase transition.For smaller µ, instead, we observe a more divergent behaviour with y = 1.21 ± 0.02 at µ = 0.2 (see Fig. 3(b)) and, ultimately, y = 1.93 ± 0.02 for µ = 0.175, which is consistent with a first-order phase transition (see Fig. 3(c)).
Taken together, these findings consistently indicate the presence of a critical value 0.175 ≤ µ * < 0.2 at which the two phase transitions merge into a single first-order transition.At the same time, they also suggest the presence of a tricritical point 0.175 ≤ µ tric < 0.2 at which the Z 2 second-order Ising transition becomes first order.Our data seem to indicate that for the modified XY model µ tric ≡ µ * .At present, however, we cannot rule out the possibility that, although they are very close, µ tric > µ * .The complete phase diagram of the model ( 1) is shown in Fig. 4(a).For µ < µ c the BKT critical temperatures are those derived in our previous work [13].In the regime µ c < µ ≤ 0.175, the critical temperatures of the first-order phase transition have been computed by a finite-size scaling analysis of the temperatures corresponding to the specific-heat peak C max v (T c , L) (see Fig. S7).According to Fig. 4(a), for µ < µ c the system exhibits a single BKT phase transition from a quasi-long-range ordered superconducting state to a disordered one.By increasing the value of µ at low temperatures, the vortex fugacity increases until, at 0.14 < µ c < 0.145, the system undergoes a first-order phase transition [5] from a vortexvacuum superconductor to a vortex supersolid which additionally breaks the Z 2 discrete symmetry associated with the two possible vortex-antivortex crystal configurations.
By increasing the chemical potential above the critical value µ c , we find that up to a value of µ * > µ c there exists a single first-order transition line separating the vortex-antivortex SC crystal from the high-temperature disordered state.For µ > µ * , instead, the two phase transitions split apart with T BKT < T I .In this regime, a new intermediate phase appears where the system is a non-superconducting vortex-antivortex crystal spontaneously breaking the Z 2 symmetry associated with the charge ordering.
Differently from the 2D Coulomb gas counterpart [4], however, the region of the phase diagram hosting this new phase is quite small and the two transitions remain close for all values of µ studied.Nonetheless, the splitting between the two transitions ∆T c = T I − T BKT increases almost linearly with µ (see Fig. 4(b)).Via a linear fit of ∆T c vs µ, we also extracted an estimate of µ * at which the two transitions merge.The obtained value µ * = 0.192 is consistent with the analysis reported above.

CONCLUSIONS
In this study, we conducted a comprehensive numerical investigation of the modified XY model by introducing a plaquette term to control the fugacity of vortices.Our findings reveal that as the vortex fugacity increases, the low-temperature superfluid BKT state turns into a vortex supersolid with finite superconducting density and charge ordering.At low temperatures, this state emerges from the superconducting vacuum via a first-order phase transition.However, as the temperature increases, a complex phase diagram emerges.At temperatures T 1 and chemical potential µ ≤ 0.14, a BKT transition line branches out of the first-order line, and vortex unbinding destroys the superconducting order.The transition line separating this new disordered state from the superconducting crystal remains first order up to µ * ≈ µ tric , while for larger µ an increasing temperature leads to the vanishing of superfluid order via the BKT mechanism, followed by the melting of the normal vortex-antivortex crystal into the disordered state via an Ising-like second-order line, as shown in Fig. 4.
Our results are consistent with the analysis conducted in Ref. [5] for the two-dimensional Coulomb gas, but two  1).The light-blue and green areas indicate the two possible low-temperature states of the system.For µ ≤ 0.14 this is a vortex-vacuum superconducting state (green area), while for µ ≥ 0.145 it turns into a vortex-antivortex superconducting crystal (light-blue area).The BKT critical points (green triangles) in the region µ < µ c of the phase diagram are those derived in our previous work [13] and separate the vortex-vacuum SC state from the disordered high-temperature state.In the region µ c < µ < µ * the system undergoes a single first-order transition (red dots) from a vortex-antivortex SC crystal to a disordered non-SC state.important differences stand out: 1. First, the area between the two transition lines separating the superconducting crystal from the normal crystal and the disordered state is extremely small and only grows linearly by increasing the chemical potential.
2. Second, the branching point of the second BKT line coincides within our numerical precision with the tricritical point µ tric , where the first-order line meets the second-order Ising transition.
These differences may be attributed to the intrinsic differences between the two Hamiltonians, particularly to the fact that the topological excitations, i.e., the vortices, are coupled to the low-energy spin waves in the XY model, while this interaction is neglected in the Coulomb gas representation of the problem.Additionally, while our study focuses primarily on the superfluid stiffness J s , Ref. [5] characterizes the superconductor by the inverse dielectric constant.These two quantities are closely related in the traditional XY model with µ = 0, but the same relation does not hold in this study, where the plaquette term in the Hamiltonian (1) gives an explicit contribution to the superfluid stiffness.
In conclusion, resolving the nature of the unconventional tricritical point, where the first-and second-order lines meet with the infinite-order BKT line, requires the derivation of an improved BKT flow equation that can capture the mechanism of defect unbinding at finite fugacity.Such a theoretical framework should be able to capture both BKT scaling and the second-order transition line within the same formalism, and its development represents the most significant future direction of this work.At the same time, research on chiral magnets with strong Dzyaloshinskii-Moriya interactions (see [14] and references therein), as well as experimental realizations of cold-atom systems with related phase diagrams [37] or spin-torque interactions [38], can provide a complementary experimental route to investigate the nature of such unconventional tricritical point.

FIG. 1 :
FIG. 1: Monte Carlo results for the case µ = 0.3.(a) Determination of T BKT from finite-size scaling of the superfluid stiffness J s renormalized according to the BKT scaling Eq. (7) for L = 64, 96, 128, 192, 256 (from top to bottom).The best crossing point is obtained with L 0 = 10.5.As expected, it lies on the 2T /π critical line (continuous black line).The dotted line indicates the extracted BKT critical temperature.(b) Determination of the Ising critical temperature T I from finite-size scaling of the Binder cumulant U stag defined in Eq.(10).The crossing point locates T I , indicated here with a dashed-dotted line.(c) Olsson's plot[9] for different values of the system size L. At the BKT critical point, while the superfluid stiffness jumps from J s (T − BKT ) = 2T − BKT /π to J s (T + BKT ) = 0, the staggered magnetization is observed to increase with L and reaches a finite value in the thermodynamic limit.This is an additional confirmation that T BKT < T I in this case.The error bars are computed via a standard bootstrapping resampling method.Where not visible, the error bars are smaller than the point symbols.

FIG. 2 :
FIG.2: Evidence for a first-order transition: the energy-density distribution P (E/N ) is shown for different system sizes L (N = L 2 ) at the temperature corresponding to the specific-heat peak.While for (a) µ = 0.175 there are two peaks indicating a first-order transition, for (b) µ = 0.2 in the thermodynamic limit a single peak emerges consistent with a continuous transition second-order transition.

FIG. 4 :
FIG. 4: (a) µ − T phase diagram of the model (1).The light-blue and green areas indicate the two possible low-temperature states of the system.For µ ≤ 0.14 this is a vortex-vacuum superconducting state (green area), while for µ ≥ 0.145 it turns into a vortex-antivortex superconducting crystal (light-blue area).The BKT critical points (green triangles) in the region µ < µ c of the phase diagram are those derived in our previous work[13] and separate the vortex-vacuum SC state from the disordered high-temperature state.In the region µ c < µ < µ * the system undergoes a single first-order transition (red dots) from a vortex-antivortex SC crystal to a disordered non-SC state.Finally, for µ > µ * the two phase transitions separate and an intermediate non-SC state with a finite Z 2 crystalline order appears.To highlight the splitting of the two critical temperatures in this region of the phase diagram, in panel (b) we report the value of T I − T BKT as a function of the chemical potential µ.A linear fit of the data (dashed grey line) gives an estimate µ * = 0.192 of the value at which the two phase transitions merge.The error bars are computed via a standard bootstrapping resampling method.Where not visible, the error bars are smaller than the point symbols.
FIG. 4: (a) µ − T phase diagram of the model (1).The light-blue and green areas indicate the two possible low-temperature states of the system.For µ ≤ 0.14 this is a vortex-vacuum superconducting state (green area), while for µ ≥ 0.145 it turns into a vortex-antivortex superconducting crystal (light-blue area).The BKT critical points (green triangles) in the region µ < µ c of the phase diagram are those derived in our previous work[13] and separate the vortex-vacuum SC state from the disordered high-temperature state.In the region µ c < µ < µ * the system undergoes a single first-order transition (red dots) from a vortex-antivortex SC crystal to a disordered non-SC state.Finally, for µ > µ * the two phase transitions separate and an intermediate non-SC state with a finite Z 2 crystalline order appears.To highlight the splitting of the two critical temperatures in this region of the phase diagram, in panel (b) we report the value of T I − T BKT as a function of the chemical potential µ.A linear fit of the data (dashed grey line) gives an estimate µ * = 0.192 of the value at which the two phase transitions merge.The error bars are computed via a standard bootstrapping resampling method.Where not visible, the error bars are smaller than the point symbols.

S3.
FIG.S2: Monte Carlo results for the case µ = 0.275.(a)Finite-size scaling of the superfluid stiffness J s renormalized according to the BKT scaling Eq.(7).(b) Finite-size scaling of the Binder cumulant U stag .(c) Finite-size Olsson's plot[9].It shows that at the BKT critical point, the staggered magnetization is finite in the thermodynamic limit.Thus confirming that T BKT < T I .