Disordered hyperuniformity in superconducting vortex lattices

PHYSICAL REVIEW RESEARCH 2, 033133 (2020) Disordered hyperuniformity in superconducting vortex lattices José Benito Llorens,1Isabel Guillamón ,1Ismael G. Serrano ,2,3Rosa Córdoba ,2,3,4Javier Sesé ,2,3 José María De Teresa ,2,3M. Ricardo Ibarra ,2,3Sebastián Vieira,1Miguel Ortuño,5and Hermann Suderow 1 1Laboratorio de Bajas Temperaturas y Altos Campos Magnéticos, Unidad Asociada UAM, CSIC, Departamento de Física de la Materia Condensada, Instituto Nicolás Cabrera and Condensed Matter Physics Center, Universidad Autónoma de Madrid, 28049 Madrid, Spain 2Laboratorio de Microscopías Avanzadas (LMA), Instituto de Nanociencia de Aragón (INA), Universidad de Zaragoza, 50018 Zaragoza, Spain 3Departamento de Física de la Materia Condensada and Instituto de Ciencia de Materiales de Aragón (ICMA), Universidad de Zaragoza-CSIC, 50009 Zaragoza, Spain 4Instituto de Ciencia Molecular, Universitat de València, Catedrático José Beltrán 2, 46980 Paterna, Spain 5Departamento de Física, CIOyN, Universidad de Murcia, Murcia 30071, Spain (Received 24 February 2020; revised 11 May 2020; accepted 22 June 2020; published 24 July 2020) The current carrying capability of type II superconductors under magnetic fields is determined to a large extent by the interaction of superconducting vortices with pinning centers. Vortices are arranged in lattices with varying degrees of disorder depending on the balance between the intervortex interactions and the pinning strength. We analyze here vortex arrangements in disordered vortex lattices of different superconducting systems, single crystals (Co-doped NbSe2, LiFeAs, and CaKFe4As4), and amorphous W-based thin films (with critical temperatures Tcfrom 4 K to 35 K and critical fields from 3.4 T to more than 90 T). We calculate for each case the structure factor and number variance and compare to calculations on an interacting set of partially pinned particles. We find that random density fluctuations appear when pinning overcomes interactions and show that the suppression of density fluctuations is correlated to the presence of interactions. We discuss the results within the framework of hyperuniform distributions and find that all studied lattices follow a similar increase of the number variance with the defect density. DOI: 10.1103/PhysRevResearch.2.033133 I. INTRODUCTION Particles occupying sites of a random lattice present density fluctuations at all length scales. It has been proposed that increasing interparticle interactions reduces long-range density fluctuations, deviating from random behavior [1,2]. The absence of density fluctuations occurs in all ordered lattices, including those with very-long-range order such as moiré patterns and quasicrystals. Density fluctuations can be described by discussing power laws in the structure factor S(k) and the number variance σ2(R). The term “hyperuniformity” was coined to gather ordered as well as disordered lattices with reduced density fluctuations [1,2]. It has been shown that disordered hyperuniform lattices can be created by artificially designing disordered patterns without density fluctuations at large length scales [1–12]. It is not yet fully clear however how to link density fluctuations with interactions in a disordered hyperuniform system. Recently, disordered hyperuniform behavior has been found in the contact number between subsystems of particles [7,8], although the origin of the suppressed density fluctuations remains under debate [13]. Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Here we address the behavior of vortices in type II superconductors. Vortices are whirlpools of currents, each carrying a flux quantum [14,15]. They repel each other and often form ordered lattices with intervortex distances a0∝1 √H. Intervortex interactions are screened above the penetration depth λ, which is most often much larger than the intervortex distance a0[16–19]. In thin films, λstrongly increases with decreasing thickness and vortices interact through their stray field, which leads to a long range Coulomb-like 1 rinteraction [20–22]. Vortices are pinned at defects in superconducting crystals and these defects are often randomly distributed. How the vortex lattice reacts to pinning depends on pinning strength and distribution as well as stiffness of the vortex lattice and determines the capability of type II superconductors to carry large currents. The degree of order in superconducting vortex lattices is determined by the balance between pinning and interaction strengths and is related to the critical current [14,15]. The question we address here is to what extent intervortex interactions can eliminate density fluctuations while allowing the formation of a disordered lattice due to pinning. The influence of vortex arrangements with suppressed density fluctuations in the current carrying capability of a superconductor was analyzed theoretically in Ref. [11]. Arrangements of vortices in cuprate superconductors at very small magnetic fields where weak intervortex interactions are weak were discussed in Ref. [23]. We analyze here vortex lattices in bulk superconductors and in thin films obtained in Refs. [24–28] and provide 2643-1564/2020/2(3)/033133(9) 033133-1 Published by the American Physical Society JOSÉ BENITO LLORENS et al. PHYSICAL REVIEW RESEARCH 2, 033133 (2020) additional data in a thin film. We mostly focus on lattices with large amounts of disorder. We calculate S(k) and σ2(R) and discuss the power laws as a function of kand R.We discuss data on a nanostructured W-based thin film (W-film-1) with strong pinning and polycrystalline vortex lattice arrangements that are disordered at high magnetic fields. We show results in Co-doped 2H-NbSe2, a system with strong pointlike pinning centers that also leads to polycrystalline hexagonal lattices which increase the level of disorder when increasing magnetic fields [25]. We discuss a W-based thin film with weak one-dimensional (1D) disorder potential created by the discommensuration between the vortex lattice and a linear nanostructure. This film shows ordered vortex lattices consistent with a Bragg glass at low magnetic fields that are disordered at high magnetic fields (W-film-2) [24,29]. We consider that the ordered lattices are representative of the behavior expected in 2H-NbSe2or other materials with wellordered lattices [30–34]. Finally, we discuss two iron-based superconductors with strongly disordered lattices and strong pointlike intrinsic pinning at all magnetic fields, LiFeAs and CaKFe4As4[26,27]. We find a gradual appearance of density fluctuations with increasing disorder in the lattice that successfully describes all studied vortex lattices. II. METHODS To analyze interactions in vortex lattice images we first have to recall the meaning of the spatial dependence of the structure factor and number variance. We show in Figs. 1(a) and 1(b) images including matrices of points distributed on a two-dimensional square. Points are positioned randomly [Fig. 1(a)] and following a hyperuniform distribution [Fig. 1(b)]. Both images are square and have the same size, with area A, and the same number of points N. In a random distribution [Fig. 1(a)], there are density fluctuations at all length scales. In a hyperuniform distribution [Fig. 1(b)], density fluctuations disappear at large length scales. We have generated the hyperuniform distribution of Fig. 1(b) by starting with points arranged in a square lattice of constant a0=√A/N. We then add to each lattice point a vector  r with random coordinates, whose absolute value r<2a0[inset of Fig. 1(b)]. The result is a random distribution which is spatially uniform for length scales larger than the intercell distance. It is important to realize that there are no further signatures of the ordered lattice in Fig. 1(b) other than the large-scale uniformity. The calculation of the number variance is very useful to discuss density fluctuations [1,2,11]. It is given by σ2(R)= N2(R)−N(R)2, where N(R) is the number of points inside a circle of radius R. The variance σ2(R) increases with R as a power law σ2(R)∝Rβ. In a random distribution of points [Fig. 1(a)], σ2(R) grows as the dimension (i.e., with the area, according to the large number law) so that β=2 [Fig. 1(c)]. On the other hand, in a hyperuniform distribution [Fig. 1(b)], σ2(R) grows as the dimension minus one (i.e., with the perimeter) with β=1 [Fig. 1(d)]. Ordered lattices provide a strong decay for small wave vectors and an oscillating σ2(R) with a dip each time Ris somewhat smaller than integer multiples of the average interparticle distance [35]. FIG. 1. Set of points distributed in a 2D space, obtained by assigning some pixels the value 1 following (a) a random distribution and (b) a hyperuniform distribution. Note that both distributions of points are fully disordered. In the inset of (b) we show schematically how we created the hyperuniform distribution. We started from a square lattice (gray crosses) with lattice constant a0and added to each lattice point a vector with random coordinates  r=(x,y) and magnitude r<2a0. This gives the distribution of black dots. (c) and (d) Variance of the distribution σ2as a function of the radius normalized to a0, calculated as explained in the text. Here σ2(R) follows a power law with the exponent β(red dotted line). Result are shown for (c) the random distribution shown in (a) and (d) the hyperuniform distribution shown in (b). The structure factor S(k) normalized to one for large values of the reciprocal space vector k is shown for (e) the random distribution shown in (a) and (f) the hyperuniform distribution shown in (b). Here S(k) also follows a power law with exponent α[red dotted lines in (e) and (f)]. To calculate the structure factor S(k) we use the Fourier transform of the image of the vortex positions and make the radial average over the Fourier transform. In the Fourier space, the structure factor S(k) decreases at small kwith the power law S(k)∝kα. In the random distribution α=0 [Fig. 1(e)], while in the hyperuniform distribution α>1 [Fig. 1(f)]. There is a relation between the exponents of S(k) and those of σ2(R). When 0 <α<1, β=2−α, and when α>1, β remains locked at 1. Ordered lattices show a Bragg peak at the reciprocal lattice wavelength and a strongly decaying S(k) for small k. 033133-2 DISORDERED HYPERUNIFORMITY IN SUPERCONDUCTING … PHYSICAL REVIEW RESEARCH 2, 033133 (2020) In the Appendix we illustrate in greater detail how to obtain S(k) and σ2(R). We first Delaunay triangulate all vortex positions, eliminating missing bonds in the sides of the image and carefully verifying that every vortex was correctly identified, as in Ref. [24]. We then use the resulting matrix of points to perform the calculation of S(k) and σ2(R) and obtain the exponents as described above (see Fig. 1). In the Appendix we describe in detail how to calculate σ2(R) from an image of a disordered vortex lattice (Fig. 4in Appendix). As the results in W-film-1 have not been published, we provide a few experimental details. We acquire the image in zero-field cooled conditions at 100 mK using the system described in Refs. [24,28,36]. The sample has been made using a focused ion beam assisted deposition and has a composition similar to the composition of W-film-2, which has perfectly ordered lattices in a large range of magnetic fields, described in Ref. [24]. The critical temperature is 5 K [37,38]. However, contrary to W-film-2, here the substrate has strong random thickness modulations at a length scale of about 200 nm, which considerably enhance pinning. The vortex lattice rearranges accordingly, showing a polycrystalline pattern, which we discuss in the Appendix (Fig. 5). We perform the numerical simulation using the following procedure, which tries to mimic the dynamics of real vortices as much as possible. Vortices are represented by interacting points. We distribute Npoints at random on a square sample with periodic boundary conditions. We choose N=2000, which is a number large enough to avoid undesirable boundary effects and small enough to converge quickly under our dynamics. We fix a certain percentage of points (between 10% and 50%) in their initial random positions, while the rest points are free to move and tend to minimize their potential energy. We use a 1/rinteraction, which converges faster than the usual intervortex interaction and produces a similar spatial distribution. We do not take into account the fact that the intervortex interaction is screened above λ, although this should not greatly influence the comparison to the experiment. The procedure selects a point at random and chooses a possible new position displaced by a small distance from the original position in a random direction. If the new position minimizes the system energy, the point is moved to the new position; otherwise is kept in the old position. These energy relaxing displacements are repeated iteratively. After a number of displacements equal to 1000 times the number of points, changes in positions are minute, since the system has reached a metastable state, and we stop the algorithm. Once a final configuration has been reached, we run a Delaunay triangulation algorithm to determine the nearest neighbors of any point and calculate the proportion of points with a number of neighbors different from 6, in order to be able to compare directly with experiment. Finally, we calculate S(k) and σ2(R), obtaining βfrom the slope of σ2(R), as discussed previously for the experimental results. III. RESULTS In Fig. 2we show our results. We start with lattices in Codoped 2H-NbSe2[Figs. 2(a) and 2(b)][25]. We find S(k)∝kα with α⩾1 and β=1. We then show results in the amorphous W-based film (W-based-1) in Figs. 2(c) and 2(d). We find a FIG. 2. The left column shows the structure factor S(k) obtained for vortex lattices in many different materials as a function of the reciprocal lattice vector kin units of 2π a0in each image. The right column shows the variance σ2(R) obtained in the same materials. A few images of vortex positions are shown between graphs. (a) and (b) Results in Co-doped NbSe2with data taken from Ref. [25]. The lateral size of the images is 375 nm. (c) and (d) Results obtained in a W-based thin film with strong pinning (W-film-1). The lateral sizes of the images are 1700, 1000, and 750 nm. (e) and (f) Results in LiFeAs with data from Ref. [26]. The lateral size of the images is 500 nm. (g) and (h) Results in a highly ordered W-based thin film (W-film-2) with a very weak 1D disorder potential, with data from Ref. [24]. The lateral sizes of the images are 1000 nm (2.5 and 3 T), 550 nm (4 and 5 T), and 500 nm (5.5 T). (i) and (j) Results in pure CaKFe4As4 from Ref. [27]. The lateral sizes of the images are 400 and 470 nm. Points are joined by lines as a guide. We plot all data in logarithmic scale and provide the power-law dependences with exponents α=1 for S(k)andβas shown in the legends of the figures for σ2(R). See also Table I. 033133-3 JOSÉ BENITO LLORENS et al. PHYSICAL REVIEW RESEARCH 2, 033133 (2020) TABLE I. Superconducting parameters of the systems studied here, critical temperature Tcand critical magnetic field Hc2. Pinning is by point defects in Co-NbSe2[25] and by strong surface corrugation in W-film-1 and is in both cases strong. In the W-film-2, there is an extremely weak interaction between thevortex lattice and the film, due to a weak linear nanofabricated structure [24]. In LiFeAs pinning is by point defects [26]. In CaKFe4As4pinning is very strong [27,39– 42]. The exponent βof the variance σ2(R) is obtained from the data in Fig. 2. Random behavior is characterized by β=2 and no crystalline order. The deviation from random towards disordered hyperuniform is seen by βdecreasing from 2 towards 1 in the absence of crystalline order. Crystalline order shows the presence or absence of a visible oscillating pattern in σ2(R) (last row). Parameter Co-NbSe2W-film-1 W-film-2 LiFeAs CaKFe4As4 Tc(K) 5.7 4 4 17.5 35 Hc2(T) 3.4 6.4 6.4 15 ≈90 β1 1 1.3 1.1–1.5 1.7 oscillating σ2(R) yes yes yes-no yes-no no very similar behavior, with α⩾1 and β=1. Both of these systems thus show close to hyperuniform behavior. They also both show that there is short-range hexagonal order in all studied images at length scales well above a0. We can see this in the oscillations appearing in σ2(R) close to integers of a0. The resulting polycrystalline vortex lattices are shown in the Appendix (Figs. 5and 6). In Figs. 2(e) and 2(f) we show results in W-film-2, from Ref. [24]. At small magnetic fields, hexagonal order is nearly perfect with just a few dislocations. This provides oscillations close to integers of the lattice constant and β≈1. When increasing the magnetic field, the number of defects in the lattice increases, particularly above about 4 T. The corresponding angular and positional correlation functions show the proliferation of topological defects in the lattice [24]. At 5 T the lattice has no long-range positional nor orientational order [24]. Here we show that S(k) and σ2(R) remain with the same power-law dependences, with β=1.1 and α⩾1 for magnetic fields below or equal to 5 T. At 5T, although the lattice has neither long-range positional nor orientational order (see Ref. [24]), the length scale for orientational order is sufficiently large to provide β≈1, i.e., near disordered hyperuniform behavior. When reaching 5.5 T, the decay length for orientational order goes from about five times a0down to a couple of a0[24]. At the same time, the oscillations in σ2(R) vanish totally at 5.5 T. However, there is also a strong deviation from hyperuniformity, with an increase of βto 1.6. Thus, the onset of strong disorder leads to a tendency to form a random distribution of vortices. In Figs. 2(g) and 2(h) we show results in LiFeAs, from Ref. [26]. The vortex lattice is highly disordered above about 2 T, with no clear hexagonal patterns observed at any length scale. As discussed in Ref. [26], the structure factor has a square-shaped orientational dependence, which shows that overall there is a tendency of the vortex lattice to lock its shape to the square crystal lattice. Note that oscillations in σ2(R) are much less pronounced than in other cases, although these are clearly visible at 2 T. The coefficient αis slightly smaller than one and βis close to one at 2 T but increases FIG. 3. (a) Exponent of σ2(R), β, as a function of the standard deviation normalized by the intervortex distance SD a0. (b) Exponent βas a function of the defect density. We show results (circles) for the systems of Fig. 2, with colors similar to those in Fig. 2: Green is Co-NbSe2, orange is W-film-1, blue is LiFeAs, red is W-film-2, and violet is CaKFe4As4. Open circles indicate results where σ2(R) has an oscillatory behavior, showing the presence of hexagonal order. Closed circles show situations with a fully disordered vortex lattice and no oscillations in σ2(R). Gray points show the result of the numerical simulation discussed in the text (the gray dashed line is a guide to the eye). with the magnetic field. This situation is close to a disordered hyperuniform arrangement. Note that, despite the presence of disorder, the vortex interaction is not negligible since there the tendency of the disordered vortex lattice is to show square arrangements at small distances [26]. Next we analyze highly disordered vortex lattices in CaKFe4As4,fromRef.[27]. We find [Figs. 2(i) and 2(j)]a considerable deviation from disordered hyperuniform behavior, with βclose to or larger than 1.5 and αless than 1. We find no signatures of oscillations in σ2(R). We summarize all results in Table I. The actual positions of vortices in all discussed lattices are given in the Appendix (Fig. 7). We can now plot (Fig. 3)βas a function of the standard deviation (SD) in the nearest-neighbor vortex positions normalized to the intervortex distance a0,SD a0, and the density of defects in the vortex lattice (we call a defect a vortex with coordination number different from 6). If we start from an ordered lattice, we are close to β=1 and SD a0as well as the defect density close to zero. We see that when there are oscillations in σ2(R) (open points in Fig. 3), β=1 although the SD a0and the defect density can be quite large. Notice that there are no data with β=1 and SD a0larger than about 30–40% of the intervortex distance. This ressembles a Lindemann criterion. Above a certain fluctuation amplitude (SD a0), the ordered lattice is unstable. When we have a randomly disordered vortex lattice, we expect βtending towards 2 and large values for SD a0and of the defect density. This indeed occurs for the fully disordered lattices of CaKFe4As4, with LiFeAs at 11 T and with Wfilm-2 at 5.5 T (closed circles in Fig. 3). For fully disordered hyperuniform or close to hyperuniform behavior, we expect β close to 1. Either SD a0or the density of defects or both should be large. In LiFeAs at 2 T we observe βclose to 1 and a small SD a0but a large number of defects. As shown in Ref. [26], the disordered lattice is locked to the crystal lattice, following its 033133-4 DISORDERED HYPERUNIFORMITY IN SUPERCONDUCTING … PHYSICAL REVIEW RESEARCH 2, 033133 (2020) orientation. In W-film-2 at 5 T we observe similar parameters, βclose to 1 and a small SD a0but a large number of defects. There are oscillations in σ2(R) highlighting that orientational order is maintained up to several a0. Thus, the appearance of disordered hyperuniformity, or the decrease in density fluctuations in disordered lattices, is linked to the presence of interactions, either in the form of locking to the crystal lattice (LiFeAs) or in short-range orientational order (W-film-2). We have performed a numerical simulation to get lowenergy configurations of interacting vortices with a fixed proportion of fixed vortices, as described in Sec. II.InFig.3(gray points) we show the results for the exponent βas a function of the standard deviation normalized by the intervortex distance [Fig. 3(a)] and of the defect density [Fig. 3(b)]. The agreement with experimental results is fairly good, giving support to the idea that the degree of hyperuniformity is associated with the number of defects in the lattice. Numerical simulations also serve to obtain an indication of the number of defects in terms of the number of pinning points included in the calculation. In the Appendix (Fig. 8) we show how these two quantities are related. We note that the number of defects is clearly correlated with the number of pinning sites and it is a bit larger than the latter, especially at low densities of pinning sites. In vortex lattices, the number of defects is a consequence of the balance between the pinning energy and the intervortex interactions. For example, in LiFeAs, disorder appears relatively far from Hc2, with a vortex lattice that interacts strongly through the crystalline lattice, as shown by the fourfold symmetry in the structure factor discussed in Ref. [26]. This leads to βclose to 1, with a large number of defects, although SD a0is maintained to relatively small values. When increasing the magnetic field there is a strong tendency towards random behavior. In W-film-2 with βclose to 1, strong disorder appears at 5 T (which is 0.78Hc2), although with a smaller defect density and with enough orientational order to present maxima in σ2(R). At a slightly larger magnetic field, at 5.5 T (which is 0.85Hc2), the lattice of W-film-2 yields to random disorder and βincreases to 1.5, with practically the same SD a0 but a larger number of defects. We note that the model produces lattices that have very small values of SD a0in Fig. 3(a). These lattices do not appear in the experiment (no points below SD a00.1). When the vortex lattice still has short-range positional order (open circles in Fig. 3), β=1. In this group of data, we have lattices that show defects and polycrystalline arrangements (W-film-1 and Co-NbSe2, shown in Fig. 3), as well as ordered lattices with a very small number of defects (open red circles of W-film2inFig.3). In the latter the positional correlations decay exponentially with distance [24], which explains why the SD a0remains above 0.1. Fully ordered hexagonal lattices (or the vortex Bragg glass with algebraically decaying positional correlation) provide β≈1 for close to zero SD a0in the representation of Fig. 3[24,29,43–46]. IV. DISCUSSION The vortex lattice arrangements are a consequence of the balance between elastic and pinning energies. In Co-doped NbSe2and in the W-film-1 thin film, pinning is strong but structured, leading to hexagonal vortex clusters observed at all magnetic fields. On the other hand, in CaKFe4As4pinning centers are so strong and randomly distributed that the vortex lattice is essentially randomly disordered in the whole range of magnetic fields studied [27,39–42]. In W-film-2, where the disorder potential is very weak, the vortex lattice is only disordered when it is very soft, very close toHc2, but then the vortex distribution shows a strong tendency to disorder randomly, because intervortex interactions are very weak. In LiFeAs we have qualitatively the same behavior as in CaKFe4As4,but with weaker pinning. Furthermore, the intervortex interaction with nonlocal contributions due to the influence of the crystal lattice symmetry [18,47] is still important and responsible for decreasing density fluctuations. The vortex lattice of LiFeAs depends strongly on the temperature range where the magnetic field is applied. There are measurements showing hexagonal vortex lattices in the same magnetic field range [48], whereas the ones we have used here [26] and neutron scattering experiments [49] provide disordered lattices. Notice that the disordered lattices discussed here are locked to the crystal lattice. This is a rather peculiar combination of long-range interaction and disorder. Locking can be explained by nonlocal corrections to the London model that favor a fourfold vortex lattice [18,47]. It is relevant to note that the only disordered lattices with close to hyperuniform behavior (β≈1, points in Fig. 3corresponding to LiFeAs at small magnetic fields and to W-film-2 at 5 T) have relatively small SD a0. In the case of LiFeAs, data follow closely the calculations. Therefore, we expect that fully disordered vortex lattices will not fall to the behavior of lattices showing fluctuations in σ2(R) as a consequence of short-range order (open circles with β≈1 in the representation of Fig. 3), but rather follow the smooth increase of βwith the amount of disorder predicted by the model. Calculations show that the vortex glass can present disordered hyperuniformity in a range of magnetic field and temperatures in the presence of strong repulsive interactions and quenched disorder [11]. Those authors proposed a phase diagram with the close to disordered hyperuniform behavior in between the Bragg glass and the random vortex glass. Our results confirm indeed the presence of this intermediate state and show that it can be obtained as a balance between interaction and pinning. In a recent work, vortex lattices at very small magnetic fields were analyzed in view of their hyperuniform properties. Those authors analyzed images of the high-Tccuprate superconductor Bi2Sr2CaCu2O8+δwith magnetic Bitter decoration in the presence of disorder [23]. At small magnetic fields, vortices are very far apart and their mutual repulsion is small [50,51]. Vortex arrangements are then strongly influenced by their interaction with pinning centers. Furthermore, at high temperatures, close to the transition to the normal state, the vortex lattice melts, leading to the vortex liquid which is a dynamic tangle of vortices [52,53]. Therefore, experiments at small magnetic fields are made by cooling from the liquid phase, which results in quenched vortex arrangements [54]. The authors of Ref. [23] concluded that long-wavelength fluctuations are systematically suppressed in the vortex lattice at small magnetic fields, as a consequence of the hydrodynamic 033133-5 JOSÉ BENITO LLORENS et al. PHYSICAL REVIEW RESEARCH 2, 033133 (2020) (b)(a) FIG. 4. Example for the calculation of σ2(R)andS(k)inavortex lattice image. (a) Image obtained in CaKFe4As4at 4 T. Vortices are the yellow colored patches. The vortex lattice is disordered. (b) Black points show the vortex positions of (a). These are obtained by calculating the center of mass of each white patch in (a), obtained as described in the text. In the inset we show the Fourier transform of the image of the vortex positions. To calculate σ2(R), we produce randomly located circles with varying radius R, such as those shown schematically in different colors. Circles with the same radii are shown in the same color. properties of the liquid phase, which leads to deviations from fully random vortex distributions. V. CONCLUSION In summary, we have analyzed the conditions for the formation of disordered hyperuniform vortex lattices in superconductors at high magnetic fields. The vortex lattice shows a tendency away from random behavior and towards hyperuniformity when the number of pinning centers is between 30% and 40%. We find that the length scale of the interaction plays a minor role in determining the strength of density fluctuations. Instead, the balance between pinning and vortex lattice stiffness controls density fluctuations, with a continuous variation between density fluctuations and disorder in the lattice. We show that the decreased density fluctuations require intervortex interactions. We conclude that we can identify emergent correlations in a vortex lattice using the structure factor S(k) and number variance σ2(R) and show that disordered vortex lattices are characterized by suppressed density fluctuations. ACKNOWLEDGMENTS We are very grateful to Tetsuo Hanaguri for sharing with us data in raw form. We also acknowledge discussions with Charles Reichhardt, who directed our attention to this problem, and with Yanina Fasano and Maria Iavarone. This work was supported by the Spanish State Agency for Research (Grants No. FIS2017-84330-R, No. MAT2017-82970-C2-1R, No. MAT2017-82970-C2-2-R, No. RED2018-102627-T, No. RYC-2014-15093, and No. CEX2018-000805-M), the Aragón Regional Government (Construyendo Europa desde Aragón) through Projects No. E28_20R with European Social Fund funding, the Comunidad de Madrid through program NANOMAGCOST-CM (Grant No. S2018/NMT-4321), and EU Program No. CA16218 (Nanocohybri). I.G. acknowledges support from the European Research Council (e) FIG. 5. Vortex lattice in the W-film-1 thin film for different values of the magnetic field: (a) 0.03 T, (b) 0.07 T, (c) 0.09 T, (d) 0.11 T, and (e) 5 T. Vortices are shown as black regions. Red bars are 120 nm long. The position of each vortex, identified by the method mentioned in the text, is shown by blue dots. Vortices with fewer than five nearest neighbors are shown as green dots and with seven nearest neighbors as white dots. Pairs of such vortices provide one dislocation. Black lines provide the Delaunay triangulation of the vortex lattice. Hexagons show the orientation of the lattice in different parts of the images. Lattices showing different orientations are separated by yellow dashed lines. 033133-6 DISORDERED HYPERUNIFORMITY IN SUPERCONDUCTING … PHYSICAL REVIEW RESEARCH 2, 033133 (2020) FIG. 6. Results in Co-doped NbSe2,fromRef.[25], for (a) 2.5 T, (b) 2.7 T, and (c) 3.3 T. We represent vortices, triangulation, lattice orientation, and red scale bars as in Fig. 5. FIG. 7. Vortex positions in all the images analyzed in this work. (a) Vortex positions in Co-doped NbSe2, from Ref. [25]. (b) Results obtained in a W-based thin film with strong pinning (W-film-1). (c) Results in a highly ordered W-based thin film (W-film-2) with a very weak 1D disorderpotential, with data from Ref. [24].(d) Results in LiFeAs, obtained from data from Ref. [26]. (e) Results in pure CaKFe4As4, from Ref. [27]. Magnetic fields are given at each panel. Scale bars in red are 200 nm long. PNICTEYES (Grant Agreement No. 679080). M.O. acknowledges support from Fundación Séneca through Grant No. 19907/GERM/15. We also acknowledge SEGAINVEX at UAM. R.C. acknowledges the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/PR19/11700008. APPENDIX To calculate σ2(R) and S(k) we start by finding vortex positions following Ref. [24]. We maximize the contrast in the image, inserting a threshold that gives a clear view of vortices as single colored and extended disks. We calculate the center of mass of each disk and use this to identify the position of each vortex. This leads to the matrix of points shown in Fig. 4(b). We then make the Fourier transform to find S(k). To calculate σ2(R) we follow Ref. [11]. We generate circles of size Rcentered at randomly generated positions and increase Rfrom the average intervortex distance a0to nearly the size of the image. We use two conditions. First, circles have to be complete and within the image. Second, circles cannot FIG. 8. Density of defects (vortices with coordination number different from 6) vs the density of fixed points (pinned vortices) in the matrix positions from the calculations described in the text. 033133-7 JOSÉ BENITO LLORENS et al. PHYSICAL REVIEW RESEARCH 2, 033133 (2020) overlap. We schematically show a few circles in Fig. 4(b).In each circle, we count N(R), the number of vortices inside the circle, and obtain σ2(R)=N2(R)−N(R)2by averaging over many circles. When Ris small, we obtain in one trial many circles spanning the whole image. When Ris large, we obtain just a few circles. We make the calculation in such a way as to increase the number of random tries giving the center of the circles with R, taking care that we average over at least 100 circles for all R. In W-film-1, vortices get pinned by differences in the thickness of the thin film, as described in detail in Refs. [37,55]. In Figs. 5(a)–5(e) we show results with increasing magnetic fields in W-film-1. We see that the vortex lattice remains with regions showing hexagonal order at all magnetic fields. We observe that regions with different vortex lattice orientation are separated by regions with a large number of dislocations. The smallest ordered regions appear at relatively large magnetic fields. A further increase of the magnetic field towards Hc2leads to a disordered random configuration, shown in Fig. 5(e), where we can still identify hexagonally ordered regions of finite size. This lattice is accordingly still hyperuniform, but with a large density of defects (about 40%). In Fig. 6we show the vortex lattice in Co-doped NbSe2. Regions with hexagonal order are observed at all magnetic fields [Figs. 6(a)–6(c)] and the number of defects can be quite large, of about 40% [25]. For completeness we show all vortex positions used to calculate the structure factor and number variance in Fig. 2 in Fig. 7. We also show the density of defects vs the pinned points from the model calculations in Fig. 8. [1] S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003). [2] S. Torquato, Phys. Rep. 745, 1 (2018). [3] W. Man, M. Florescu, E. P. Williamson, Y. He, S. R. Hashemizad, B. Y. C. Leung, D. R. Liner, S. Torquato, P. M. Chaikin, and P. J. Steinhardt, Proc. Natl. Acad. Sci. USA 110, 15886 (2013). [4] Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C. Corbo, and S. Torquato, Phys.Rev.E89, 022721 (2014). [5] R. Kurita and E. R. Weeks, Phys.Rev.E84, 030401(R) (2011). [6] R. Dreyfus, Y. Xu, T. Still, L. A. Hough, A. G. Yodh, and S. Torquato, Phys. Rev. E 91, 012302 (2015). [7] D. Hexner, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 121, 115501 (2018). [8] D. Hexner, P. Urbani, and F. Zamponi, Phys.Rev.Lett.123, 068003 (2019). [9] C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. Lett. 106, 178001 (2011). [10] J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo, Phys. Rev. Lett. 115, 108301 (2015). [11] Q. Le Thien, D. McDermott, C. J. O. Reichhardt, and C. Reichhardt, Phys. Rev. B 96, 094516 (2017). [12] I. A. Sadovskyy, Y. L. Wang, Z.-L. Xiao, W.-K. Kwok, and A. Glatz, Phys. Rev. B 95, 075303 (2017). [13] H. Diamant, J. Club Condens. Matter Phys. (2019). [14] G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Rev. Mod. Phys. 66, 1125 (1994). [15] E. H. Brandt, Rep. Prog. Phys. 58, 1465 (1995). [16] J. Friedel, P. G. De Gennes, and J. Matricon, Appl. Phys. Lett. 2, 119 (1963). [17] P. G. de Gennes and J. Matricon, Rev. Mod. Phys. 36, 45 (1964). [18] V. G. Kogan, Phys. Rev. B 24, 1572 (1981). [19] A. I. Buzdin, A. S. Mel’nikov, and A. V. Samokhvalov, J. Supercond. Novel Magn. 26,2853 (2013). [20] J. Pearl, Appl. Phys. Lett. 5, 65 (1964). [21] E. H. Brandt, Phys.Rev.B48, 6699 (1993). [22] E. Herrera, I. Guillamón, J. A. Galvis, A. Correa, A. Fente, S. Vieira, H. Suderow, A. Y. Martynovich, and V. G. Kogan, Phys. Rev. B 96, 184502 (2017). [23] G. Rumi, J. Aragón Sánchez, F. Elías, R. Cortés Maldonado, J. Puig, N. R. Cejas Bolecek, G. Nieva, M. Konczykowski, Y. Fasano, and A. B. Kolton, Phys.Rev.Res.1, 033057 (2019). [24] I. Guillamón, R. Córdoba, J. Sesé, J. M. De Teresa, M. R. Ibarra, S. Vieira, and H. Suderow, Nat. Phys. 10, 851 (2014). [25] M. Iavarone, R. Di Capua, G. Karapetrov, A. E. Koshelev, D. Rosenmann, H. Claus, C. D. Malliakas, M. G. Kanatzidis, T. Nishizaki, and N. Kobayashi, Phys. Rev. B 78, 174518 (2008). [26] T. Hanaguri, K. Kitagawa, K. Matsubayashi, Y. Mazaki, Y. Uwatoko, and H. Takagi, Phys.Rev.B85, 214505 (2012). [27] A. Fente, W. R. Meier, T. Kong, V. G. Kogan, S. L. Bud’ko, P. C. Canfield, I. Guillamón, and H. Suderow, Phys. Rev. B 97, 134501 (2018). [28] I. Guillamón, H. Suderow, P. Kulkarni, S. Vieira, R. Córdoba, J. Sesé, J. D. Teresa, M. Ibarra, G. Shaw, and S. Bannerjee, Physica C 503, 70 (2014). [29] T. Giamarchi and P. Le Doussal, Phys. Rev. Lett. 72, 1530 (1994). [30] H. F. Hess, R. B. Robinson, R. C. Dynes, J. M. Valles, and J. V. Waszczak, Phys. Rev. Lett. 62, 214 (1989). [31] I. Guillamón, H. Suderow, S. Vieira, L. Cario, P. Diener, and P. Rodière, Phys. Rev. Lett. 101, 166407 (2008). [32] O. Fischer, M. Kugler, I. Maggio-Aprile, C. Berthod, and C. Renner, Rev. Mod. Phys. 79, 353 (2007). [33] H. Suderow, I. Guillamón, J. G. Rodrigo, and S. Vieira, Supercond. Sci. Technol. 27, 063001 (2014). [34] A. M. Troyanovski, J. Aarts, and P. H. Kes, Nature (London) 399, 665 (1999). [35] M. A. Klatt, J. Kim, and S. Torquato, Phys.Rev.E101, 032118 (2020). [36] H. Suderow, I. Guillamon, and S. Vieira, Rev. Sci. Instrum. 82, 033711 (2011). [37] I. Guillamón, H. Suderow, S. Vieira, A. Fernández-Pacheco, J. Sesé, R. Córdoba, J. M. D. Teresa, and M. R. Ibarra, New J. Phys. 10, 093005 (2008). [38] R. Córdoba, T. I. Baturina, J. Sesé, A. Y. Mironov, J. M. De Teresa,M.R.Ibarra,D.A.Nasimov,A.K.Gutakovskii,A.V. Latyshev, I. Guillamón, H. Suderow, S. Vieira, M. R. Baklanov, J. J. Palacios, and V. M. Vinokur, Nat. Commun. 4, 1437 (2013). [39] S. J. Singh, M. Bristow, W. R. Meier, P. Taylor, S. J. Blundell, P. C. Canfield, and A. I. Coldea, Phys. Rev. Mater. 2, 074802 (2018). 033133-8 DISORDERED HYPERUNIFORMITY IN SUPERCONDUCTING … PHYSICAL REVIEW RESEARCH 2, 033133 (2020) [40] S. Pyon, A. Takahashi, I. Veshchunov, T. Tamegai, S. Ishida, A. Iyo, H. Eisaki, M. Imai, H. Abe, T. Terashima, and A. Ichinose, Phys. Rev. B 99, 104506 (2019). [41] N. Haberkorn, M. Xu, W. R. Meier, J. Schmidt, S. L. Bud’ko, and P. C. Canfield, Phys.Rev.B100, 064524 (2019). [42] S. Ishida, A. Iyo, H. Ogino, H. Eisaki, N. Takeshita, K. Kawashima, K. Yanagisawa, Y. Kobayashi, K. Kimoto, H. Abe, M. Imai, J.-i. Shimoyama, and M. Eisterer, npj Quantum Mater. 4, 27 (2019). [43] T. Klein, I. Joumard, S. Blanchard, J. Marcus, R. Cubitt, T. Giamarchi, and P. Le Doussal, Nature (London) 413, 404 (2001). [44] T. Giamarchi and P. Le Doussal, Phys. Rev. B 52, 1242 (1995). [45]D.S.Fisher,M.P.A.Fisher,andD.A.Huse,Phys.Rev.B43, 130 (1991). [46] J. Kierfeld, T. Nattermann, and T. Hwa, Phys. Rev. B 55, 626 (1997). [47] V. G. Kogan, M. Bullock, B. Harmon, P. Miranovi´ c, L. Dobrosavljevi´ c-Gruji´ c, P. L. Gammel, and D. J. Bishop, Phys. Rev. B 55, R8693 (1997). [48] S. S. Zhang, J.-X. Yin, G. Dai, H. Zheng, G. Chang, I. Belopolski, X. Wang, H. Lin, Z. Wang, C. Jin, and M. Z. Hasan, Phys. Rev. B 99, 161103(R) (2019). [49] D. S. Inosov, J. S. White, D. V. Evtushinsky, I. V. Morozov, A. Cameron, U. Stockert, V. B. Zabolotnyy, T. K. Kim, A. A. Kordyuk, S. V. Borisenko, E. M. Forgan, R. Klingeler, J. T. Park, S. Wurmehl, A. N. Vasiliev, G. Behr, C. D. Dewhurst, and V. Hinkov, Phys.Rev.Lett.104, 187001 (2010). [50] R. Prozorov, A. F. Fidler, J. R. Hoberg, and P. C. Canfield, Nat. Phys. 4, 327 (2008). [51] J. B. Llorens, L. Embon, A. Correa, J. D. González, E. Herrera, I. Guillamón, R. F. Luccas, J. Azpeitia, F. J. Mompeán, M. García-Hernández, C. Munuera, J. A. Sánchez, Y. Fasano, M. V. Milosevic, H. Suderow, and Y. Anahory, Phys. Rev. Research 2, 013329 (2020). [52] R. Cubitt, E. M. Forgan, G. Yang, S. L. Lee, D. M. Paul, H. A. Mook, M. Yethiraj, P. H. Kes, T. W. Li, A. A. Menovsky, Z. Tarnawski, and K. Mortensen, Nature (London) 365, 407 (1993). [53] E. Zeldov, D. Majer, M. Konczykowski, V. B. Geshkenbein, V. M. Vinokur, and H. Shtrikman, Nature (London) 375, 373 (1995). [54] M. Marchevsky, P. Kes, and J. Aarts, Physica C 282–287, 2083 (1997). [55] I. G. Serrano, J. Sesé, I. Guillamón, H. Suderow, S. Vieira, M. Ibarra, and J. M. De Teresa, Beilstein J. Nanotechnol. 7,1698 (2016). 033133-9