Uniformly frustrated XY model: strengthening of the vortex lattice by intrinsic disorder

In superconducting films, the role of intrinsic disorder is typically to compete with superconductivity by fragmenting the global phase coherence and lowering the superfluid density. Nonetheless, when a transverse magnetic field is applied to the system and an Abrikosov vortex lattice forms, the presence of disorder can actually strengthen the superconducting state against thermal fluctuations. By means of Monte Carlo simulations on the uniformly frustrated XY model in two dimensions, we show that while for weak pinning the superconducting critical temperature $T_c$ increases with the applied field $H$, for strong enough pinning the experimental decreasing dependence between $T_c$ and $H$ is recovered with a resulting more robust vortex lattice.

In superconducting films, the role of intrinsic disorder is typically to compete with superconductivity by fragmenting the global phase coherence and lowering the superfluid density. Nonetheless, when a transverse magnetic field is applied to the system and an Abrikosov vortex lattice forms, the presence of disorder can actually strengthen the superconducting state against thermal fluctuations. By means of Monte Carlo simulations on the uniformly frustrated XY model in two dimensions, we show that while for weak pinning the superconducting critical temperature Tc increases with the applied field H, for strong enough pinning the experimental decreasing dependence between Tc and H is recovered with a resulting more robust vortex lattice.

I. INTRODUCTION
The melting transition of two-dimensional (2D) Abrikosov vortex lattice (VL) has always attracted a significant experimental and theoretical interest. The interplay between magnetic field, random pinning and phase fluctuations makes the phase diagram of the system rich of different phases of matter ranging from the 2D Bragg glass to the 2D vortex glass [1][2][3], from the isotropic vortex liquid to the more recently observed hexatic liquid phase [4,5]. In the absence of disorder, at low temperatures and low vortex densities, the vortex lattice exhibits a quasi-long range order. As the temperature or the vortex density increases, topological defects such as dislocations and disclinations form, eventually leading to the melting of the vortex lattice according to the Berezinskii-Kosterlitz-Thouless theory [6][7][8] afterwards refined by Nelson, Halperin and Young [9][10][11]. In this frame, the role of disorder is on the one hand to favour the formation of such a lattice defects turning the ground state from a quasi ordered Bragg glass [12,13] to a disordered vortex glass. On the other hand, by acting as a pinning potential for vortices it also prevents the VL from sliding throughout the system and destroying the global superconducting (SC) phase coherence. Thus, differently from the Meissner state, where disorder competes with superconductivity, in the mixed state its role is less straightforward. Here, indeed, a true superconducting state exhibiting zero resistivity in the limit of zero current can only settle if certain degree of disorder is present within the system. Moreover, as recently shown [14], vortex thermal fluctuations responsible for the fast depletion of the superfluid stiffness can be even reduced by the presence of a strong vortex pinning.
From a theoretical perspective, one possible way to investigate such interplay is by studying the classical XY model in the presence of both disorder and transverse magnetic field. As an effective model for superconducting phase fluctuations, the 2D XY model has been over the years applied to the study of disordered and inhomogeneous SC films [15][16][17][18][19][20][21][22] as well as to the 2D vortex lattice melting for different values of the applied magnetic field [23][24][25][26][27][28][29][30]. The combined effect of quenched disorder and transverse magnetic field has been so far addressed by including in the XY model spatially random quenched magnetic fluxes, but with a resulting zero net magnetic field [2,15,31,32].
In this manuscript, we will instead consider the effect of quenched disorder on the melting of the 2D vortex lattice which forms when a nonzero net transverse magnetic field is included within the XY model. Starting from the consideration that the presence of a lattice introduces a pinning effect at low temperature even in the absence of disorder, our Monte Carlo simulations show that apart from disordering the VL ground state, the presence of inhomogeneities makes the vortex lattice stronger against thermal fluctuations with respect the homogenous case according with the experimental observations [14]. The presence of a very week pinning can indeed lead to a very counter-intuitive observation, namely the increase of the SC critical temperature by the increase of the applied magnetic field. However, when the disorder is sufficiently strong the experimentally observed dependence between T c and H is recovered.

II. MATERIALS AND METHODS
We consider a two-dimensional XY model in the presence of random pinning and transverse magnetic field. Its Hamiltonian read: where θ i is the SC phase at site i of a L x × L y lattice, J µ i are the random couplings between the neighbouring sites i and i + µ and the phase shift arXiv:2110.02275v1 [cond-mat.supr-con] 5 Oct 2021 F µ i accounts via the minimal substitution for the presence of a finite transverse magnetic field Hẑ = ∇ × A. Being: with Φ 0 = hc/2e the quantum flux, and A µ i is the vector potential along the bond connecting two neighbours spins i and i + µ. The sum of F µ i going counterclockwise around any closed path C of bonds on the lattice is 2π times the number of magnetic flux quanta f C penetrating the path: Since in the following we will always consider H to be uniform in space, we will refer to the intensity of the applied field in terms of the flux quanta f penetrating through an unitary plaquette P : In literature , one usually refers to this case as the uniformly frustrated XY model with frustration f . Indeed, the phase shift F µ i in the cosine argument of (1) adds frustration to the system by rendering the ground state no longer ferromagnetic: at T = 0, the phases θ i 's instead to be all equal, will vary from site to site trying to minimise the new gauge-invariant phase (θ i − θ i+µ + F µ i ). Consequently, the value of f will correspond to the level of such frustration, determining the inhomogeneous space structure of the ground state itself. Specifically, the ground state of the uniformly frustrated XY model will consist of a periodic configuration of vortices in the phase angle θ i , whose number is directly proportional to f . The number of vortices of the ground state for a given value of f , can be easily derived by rewriting the charge neutrality condition in terms of the new gauge invariant phase: Since in the present study we will consider periodic boundary conditions, a good gauge choice for the vector potential A is the Coulomb gauge: ∇ · A = 0. For simplicity, we will consider: so that: Finally, it is important to highlight that with this choice not all the values of f will be allowed. Indeed, periodic boundary conditions, together with the Coulomb gauge A y = 2πf x, give rise to the constraint: Hence, for a given value of L x the smallest frustration we can introduce within the system is: In the present work, we have studied the model Eq.
(1) on a square lattice with periodic boundary conditions for a linear size of L x = L y = L = 64 and different values of f . In our Monte Carlo (MC) simulations, we have used a local Metropolis algorithm, needed to probe the correct canonical distribution of the system, combined with a micro-canonical Overrelaxation algorithm. Specifically, each Monte Carlo step consists of 5 Metropolis spin flips of the whole lattice, followed by 10 Over-relaxation sweeps of all the spins. To help the correct thermalization of the system at lower temperatures, we have used a Simulated-Annealing procedure. For each run we have made 50000-75000 MC steps, measuring the main observables with a frequency of 5 steps, after having discarded the firsts 25000. Finally, the averages have been computed over 5 independent runs for the clean case and over 10 different realization of quenched disorder for the disordered case.
In the present work, together with the ground state of the vortex lattice, we have studied the SC transition as function of the frustration f and the level of disorder. To address this issue and measure the SC phase coherence of the system, we have computed the superfluid stiffness J µ s which accounts for the linear response of the system to an infinitesimal twist ∆ µ of the gauge invariant phase along a given direction µ. As such, it is defined as the second derivative of the free energy with respect to ∆ µ at ∆ µ = 0: being it finite in the SC phase and zero in the normal phase. Its expression for the model Eq.
(1) reads: where β is the inverse temperature and . . . stays both for the Monte Carlo thermal average and for the average over the independent runs.

III. RESULTS
In what follows, we will present the numerical results obtained via Monte Carlo simulations both for the clean and for the disordered case, considering different values of f ∈ [0, 1 2 ].

A. Clean case
Let us start with the clean case where J µ i = 1 ∀i, µ so that: For f = 0, the model (11) is the classical XY model which undergoes to the well known Berezinskii-Kosterlitz-Thouless [6][7][8] (BKT) transition at T = T BKT . As already mentioned, for finite values of f the phase transition is no longer driven by the unbinding of vortex-antivortex pairs, but rather by the melting of the vortex lattice which naturally forms when a transverse magnetic field is applied to the system. In the clean case, despite the absence of disorder, the vortex lattice is pinned at low temperatures by the presence of the underlying square grid that defines the array of Josephson junctions. By acting as a periodic pinning potential for vortices, such a square grid can become particularly relevant for large values of f eventually determining the symmetry of the vortex lattice itself. The most emblematic example is the limiting case f = 1/2, where in the ground state the vortex lattice assumes a checkerboard ordered pattern. The model Eq.(11) with f = 1/2 is also known as the fully frustrated XY model (FFXY) and as such it has been extensively studied (see [30] and reference therein) with particular focus on its critical behaviour and the nature of the phase transitions it undergoes.
In the present work, we will instead focus on smaller value of f , where the periodic pinning potential does not induce such peculiar vortex lattice configuration. The resulting vortex lattice found numerically in the ground state are shown in Fig.1. For all the values of f considered, at low temperatures the vortices form an ordered and pinned lattice whose symmetry partially depends also on the commensurability with the underlying numerical grid.
At finite but not maximal f different kind of phase transitions can occur [25,29], and their nature is still unclear in most of the cases. Without pretending to address this issue, in the present work we focus on the dependence of the critical temperature T c , at which the superfluid stiffness vanishes, on the applied magnetic field f . The numerical results for different values of f are reported in Fig.2.
We can immediately notice that the critical temperature T c at which the system loses its phase coherence is strongly suppressed for smaller values of f , completely at odd with the usual experimental observations [2,4,5,[33][34][35][36][37]. The general trend seems, indeed, to be a proportionality between T c and f : lower critical temperatures for smaller frustration. The observed trend has been already discussed in [29], where Alba et al. have shown that in the limit of small frustration f = 1/n (and n 1), the critical temperature T c decreases with the increase of n: T c ∼ 1/n → 0 as n → ∞.
Apart from specific cases where the commensurability with the underlying square grid particularly strengthens or weakens the vortex-lattice structure, by lowering the applied magnetic field the pinning of the vortex lattice becomes weaker than in the case where the vortex density is higher, with a resulting decrease of the critical temperature with f . By increasing the pinning potential via the introduction of disorder, however, the scenario is reversed and, in agreement with experimental observations, with increasing applied magnetic field the critical temperature decreases. We will discuss the numerical results of the disordered case in the following section.

B. Disordered case
We now consider the case of a transverse magnetic field applied to a disordered SC film. In the present manuscript, we will use as disordered coupling constants J µ i in Eq.(1) the inhomogeneous local stiffness derived from the (quantum) XY pseudo-spin 1/2 model in random transverse field (RTF) [38,39]. This disorder has been shown to be appropriate to model disordered superconductors with a non-trivial spatial structure [39][40][41], as well as to account for the experimental observation of a rather broad BKT transition around the critical temperature T BKT [21,42], at which a sharp jump of the superfluid-density would be expected for zero disorder [43]. Previous studies at zero magnetic field [21,42] have shown that such a spatiallycorrelated disorder with large enough low-stiffness puddles is crucial to induce an anomalous vortex nucleation and, consequently, a smearing of the BKT transition. On the other hand in the present case where vortices are induced by a finite transverse magnetic field we do not expect that the results will depend crucially on the choice of disorder, so our finding should be general also for different disorder realizations. The level of disorder is here labeled by W/J (see [21] for more details) and it is taken to be quenched in temperature.
Let us start by considering a relatively weak disorder level W/J = 4. Consistently with the experimental observations of [14], our MC numerical results reveal that the presence of disorder leads to a modification of the ground-state vortex lattice, enlightening further its underlying mechanism. As shown in Fig.3, the core of the vortices is indeed pinned by the inhomogeneity of the local stiffness, which makes them moving towards those regions with lower J µ i so to gain in energy by minimising the Hamiltonian(1). To highlight the correlation in space between low-couplings regions and the vortex lattice deformation, in Fig. 3 we have superimposed the vortex lattice to the couplings map, obtained by computing over each plaquette the average value of the local stiffness J µ i . The presence of disorder not only modifies the vortex lattice ground state, but also impact the superfluid response of the system as function of the applied magnetic field. In Fig.4, we report the temperature dependence of the superfluid stiffness J s (T ) obtained for different values of f . As one can see, the presence of disorder restores the experimentally measured dependence between T c and f by rendering more robust the vortex lattice against thermal fluctuations. In this regard, it is quite impressive to notice that, compared to the homogeneous case (see Fig.3), the critical temperature T c corresponding to the lowest value of f = 1/64 has increased by a factor of ten by the effects of the inhomogeneity.
The strengthening of the vortex lattice due to disorder is even more pronounced when looking at stronger disorder regimes. In Fig.5, we report the superfluid stiffness trend in temperature for the same values of f but with a disorder level of W/J = 10. Looking for instance at the lowest value of the field (f = 1/64), the critical temperature is reduced, because of the field, just by half with respect to the zero-field value, while at weak disorder (W/J = 4) it was five times smaller.
In order to highlight such increase of robustness, as effect of the increase of the intrinsic disorder, we have reported in Fig.6 (a) the extrapolated values of the critical temperature as function of the applied field, for the two disorder regimes considered. In Fig.6(b), this effect is made even more evident by rescaling the curves of T c by their value in absence of magnetic field.

IV. DISCUSSION
By means of Monte Carlo simulations on the uniformly frustrated XY model, we have shown that the presence of disorder, mimicked via random couplings, modifies the ordered vortex lattice of the ground state. For strong enough disorder, the energy of the system is no longer minimized by an ordered vortex-lattice structure, but rather by a disordered structure in which the vortex cores are located where the local stiffness of the superfluid is lower so as to reduce the energy cost associated with a phase twist.
At the same time, we have shown that in the presence of a very weak pinning as in the case with no disorder, where the only pinning potential is due to the presence of the underlying square grid, the increase of the applied magnetic field does not reduce the superconducting critical temperature, but in most cases it contributes to its increase, in stark contrast with the experimental observations [2,4,5,[33][34][35][36][37]. The presence of disorder not only restores the experimentally observed dependence between T c For each values of f = 0 and disorder level, the critical temperature Tc has been estimated from the superfluid stiffness trend in temperature Js(T ) as the temperature at which Js vanishes starting from low temperatures. For the case f = 0, we used instead the Nelson-Kosterlitz universal relation [43].
and f , but also acts by making the ground-state vortex lattice more robust against thermal fluctuations.
Comparing two different levels of disorder, we have indeed shown that the superconducting critical temperature is suppressed much less with respect to the zero-field case in the strong disorder regime than in the weak disorder regime.
In conclusion, despite separately both disorder and magnetic field act on the SC thin film by suppressing the superfluid stiffness, when a transverse magnetic field is applied to the system, the presence of disorder help to prevent the destruction of the SC state by thermal fluctuations.