Hyperuniform vortex patterns at the surface of type-II superconductors

A many-particle system must posses long-range interactions in order to be hyperuniform at thermal equilibrium. Hydrodynamic arguments and numerical simulations show, nevertheless, that a three-dimensional elastic-line array with short-ranged repulsive interactions, such as vortex matter in a type-II superconductor, forms at equilibrium a class-II hyperuniform two-dimensional point pattern for any constant-$z$ cross section. In this case, density fluctuations vanish isotropically as $\sim q^{\alpha}$ at small wave-vectors $q$, with $\alpha=1$. This prediction includes the solid and liquid vortex phases in the ideal clean case, and the liquid in presence of weak uncorrelated disorder. We also show that the three-dimensional Bragg glass phase is marginally hyperuniform, while the Bose glass and the liquid phase with correlated disorder are expected to be non-hyperuniform at equilibrium. Furthermore, we compare these predictions with experimental results on the large-wavelength vortex density fluctuations of magnetically decorated vortex structures nucleated in pristine, electron-irradiated and heavy-ion irradiated superconducting BiSCCO samples in the mixed state. For most cases we find hyperuniform two-dimensional point patterns at the superconductor surface with an effective exponent $\alpha_{\text{eff}} \approx 1$. We interpret these results in terms of a large-scale memory of the high-temperature line-liquid phase retained in the glassy dynamics when field-cooling the vortex structures into the solid phase. We also discuss the crossovers expected from the dispersivity of the elastic constants at intermediate length-scales, and the lack of hyperuniformity in the $x\,-y$ plane for lengths $q^{-1}$ larger than the sample thickness due to finite-size effects in the $z$-direction.


I. INTRODUCTION
Hyperuniform point patterns, defined by a complete suppression of density fluctuations in the largewavelength limit [1], have attracted a great interest in the last years. Such behavior can spontaneously emerge, following either equilibrium or non-equilibrium protocols, in disordered ground states, glass formation, jamming, Coulomb systems, spin systems, photonic and electronic band structure, localization of waves and excitations, self-organization, fluid dynamics, number theory, stochastic point processes, integral and stochastic geometry, photoreceptor cells, and even the immune system [1][2][3]. Hyperuniform systems are proposed to be distinguishable states of matter characterized by special properties [2]. Besides, these properties can be technologically exploited directly, or indirectly, by coupling a given system with an hyperuniform pattern. The fabrication of such patterns, either in a controlled or a self-assembled way, is hence also of interest from an applied point of view [3].
Point patterns formed by vortex matter nucleated in type-II superconductors have been a paradigmatic soft condensed matter system for studying basic questions on statistical physics, such as the statics and dynamics of elastic manifolds in random media and glassy phases in general. This system has also been a playground to understand the rich interplay between elasticity, quenched disorder, thermal fluctuations, driving forces, anisotropic and finite-size effects, either at equilibrium or out of it [4]. Since the mean vortex density can be easily controlled by changing the applied field H, vortex matter systems are particularly suitable for studying ordering and density fluctuations at microscopic scales, in different equilibrium or non-equilibrium liquid, solid, or glassy phases. However, the occurrence of hyperuniformity in vortex matter has not been experimentally addressed yet.
Recent theoretical studies on the effect of hyperuniform pinning arrays [5,6] on vortex matter report an isotropic enhancement of the critical currents in comparison with a non-hyperuniform distribution of point pins. The work of Ref. citeLeThien2017 also predicts that vortex matter may exhibit disordered hyperuniformity in the presence of hyperuniform or random pinning arrays. In particular, hyperuniform vortex configurations are proposed for rigid vortices repelling with short-range interactions in presence of a Poisson distribution of point pins, in a narrow region between the Bragg glass and the vortex liquid phases. This last result is in contrast with the behavior expected for the same rigid vortex system at thermal equilibrium in the absence of disorder. Indeed, the fluctuation-compressibility theorem forbids hyperuniformity [2] in a system with a non-divergent compression modulus for nearly-uniform deformations. On the other hand, hyperuniformity is expected for clean thin-film superconductors with the penetration length much larger than the sample thickness and logarithmic vortex-vortex repulsion. Ref. [5] shows that this expectation is satisfied also in the presence of Poisson-distributed point pins.
Interestingly, as we show in this work, threedimensional systems made of many non-rigid interacting lines directed along the z-direction are generically expected to follow an hyperuniform point pattern for any constant-z cross section of the embedding space. This behavior arises mainly from the constraint that elastic lines can not start nor terminate inside the medium, and can be understood considering general hydrodynamic arguments valid in the solid and liquid phases, even including weak pinning in the liquid [7,8]. These arguments show that, even though interactions between lines are short-ranged in all directions, the effective two-dimensional bulk modulus of the (compressible) three-dimensional system smoothly diverges in the large-wavelength limit. Hence, at thermal equilibrium, density fluctuations smoothly vanish and the configurations are hyperuniform for any particular constant-z slice. We show that by adding heuristic arguments these assertions can be extended to predict the density fluctuations in the Bragg glass, Bose glass, and liquid vortex phases with correlated disorder generated by columnar defects (CD). We find that disorder modifies in a nontrivial way the large-wavelength density fluctuations of the ideal clean system. In the former case hyperuniformity becomes marginal, while in the correlated disorder case hyperuniformity is destroyed in the glassy and liquid equilibrium phases. These hydrodynamic predictions are different from the strong type of hyperuniformity numerically predicted in Ref. 5 for logarithmic and short-range interactions in a three-dimensional system of rigid vortices in presence of quenched disorder.
Motivated by the above predictions and open questions we experimentally study large wavelength vortex density fluctuations in magnetically decorated vortex structures over extended fields-of-view (thousands of vortices) in pristine, electron-irradiated and heavy-ion irradiated (namely with CD) Bi 2 Sr 2 CaCu 2 O 8+δ superconducting samples. In the observable spatial range we systematically find, for all samples and vortex densities probed, an effective hyperuniform behavior close to the one predicted for the line liquid in equilibrium under weak disorder. We argue that this result can be explained considering that larger wavelength density fluctuations have also a slower dynamics, and are thus effectively arrested by a realistic field-cooling process with finite temperature sweep-rate. We also show that dispersivity of the elastic constants is experimentally relevant, as well as finitesize effects in the direction of the applied field, which should ultimately kill the predicted asymptotic hyperuniformity at in-plane scales of the order of the superconductor thickness.
The paper is organized as follows. In Section II A we define the main observables that we will use to study hyperuniformity in magnetically decorated vortex structures. In Section II B we review the hydrodynamic arguments supporting the emergence of hyperuniformity at any z-constant slice of a three-dimensional system of repelling elastic lines for different phases. In Section III we study the large-wavelength density fluctuations in magnetically decorated vortex structures nucleated in Bi 2 Sr 2 CaCu 2 O 8+δ samples with different types of disorder. In Section IV we discuss the theoretical interpretation of these results, and report new predictions expected for vortex matter in ad hoc experiments that might shed more light into this issue.

II. DENSITY FLUCTUATIONS: PHENOMENOLOGICAL THEORY
In this section we describe the physical magnitudes that we study and the predictions obtained for the different vortex phases using a hydrodynamic approach, scaling arguments and numerical simulations of simple models.
A. Cross-section hyperuniformity L x y r n (z) FIG. 1. Simulation-snapshot of an array of fluctuating directed elastic lines (green) modeling vortices in a type-II superconductor with the magnetic field applied in the zdirection. Two constant-z cross sections are highlighted, the top and an inner layer of the sample. The hydrodynamics of the elastic lines predicts a class-II hyperuniform twodimensional point pattern at each constant-z cross section (full circles), in spite of repulsion between lines being shortranged. The pattern frozen at the top layer can be accessed experimentally via Bitter magnetic decoration experiments (see Section.III) .
We will be interested in the two-dimensional largewavelength density fluctuations at a constant-z crosssection of a three-dimensional array of vortex lines nucleated in a real sample of thickness L. The average magnetic induction at equilibrium, B ≡ Bẑ, determines the number of vortices per unit area in any slice, n 0 ≡ B/Φ 0 , with Φ 0 = 2.07·10 −7 G.cm 2 the flux quantum. We model this system of individual vortex lines by considering directed elastic lines. We will also assume, as an approx-imation, that vortex lines do not present overhangs nor pinch-off loops, and that at a given instant the nth line can be appropriately described by a parameterized position r n (z) ≡ (x n (z), y n (z)) (see Fig.1 for a schematic representation).
The structure factor of a given frozen two-dimensional point pattern with N 1 points at a constant-z cross section is defined as with q = (q x , q y ). For thick systems in the z-direction L → ∞, and the bulk value of S(q, z) is expected to be independent of z, namely S(q, z) ≡ S(q) [9]. A constant-z cross section pattern is considered hyperuniform if its two-dimensional large-wavelength (smallwave-vector) density fluctuations are suppressed lim q→0 S(q) = 0. (2) We will be particularly interested in the common cases where the structure factor vanishes isotropically as a power-law near the reciprocal space origin, with α ≥ 0 a characteristic exponent. In the asymptotic limit N 1, this property can be translated to the variance of the number N (R) of particles observed inside an hyper-spherical window with radius R, σ 2 N (R) = N 2 (R) − N (R) 2 . In the latter . . . denotes average over randomly-distributed windows and the expression is valid in the large-R limit. The variance is predicted to scale as [1] For instance, a Poisson distribution of points is not hyperuniform since the number variance is extensive with the window volume, i.e. σ 2 N ∼ R d . A hyperuniform system is considered to be class-I if α > 1, class-II if α = 1, and class-III if 0 < α < 1. A large and diverse list of physical and mathematical systems for each class is given in

B. Hydrodynamic arguments
The equilibrium hydrodynamics of three-dimensional flux-line liquids have been thoroughly discussed in the past [7,8]. In this section we first briefly review the basic physics and also add some new predictions regarding large-wavelength isotropic density fluctuations in the Bragg glass, the Bose glass (samples with correlated disorder) and liquid phases. This theoretical framework is relevant for our discussion of hyperuniformity at constant-z cross sections since the two-dimensional structure factor describing large-scale density fluctuations can be obtained from the three-dimensional one by simple integration over the z-direction wave-vectors. This yields, at equilibrium, an effective (single-layer) two-dimensional compression modulus [10].
If the penetration depth of a superconductor is much smaller than the thickness of the sample, λ ab L, (see Fig.1) vortex lines repel each other with a roughly exponential dependence on distance for scales of order λ ab or larger. We neglect the contribution of long-range surface forces between vortex tips in a finite sample [10]). In the absence of anisotropy in the x − y plane, and for the intermediate vortex densities we are interested in, these interactions tend to form a triangular Abrikosov lattice well characterized by dispersive and axially-symmetric elastic modulii of tilt c 44 (q, q z ), compression c 11 (q, q z ), and shear c 66 (q, q z ). With increasing temperature, this solid phase melts into a liquid phase of fluctuating lines with c 66 (q, q z ) = 0. However, in the liquid c 11 (q, q z ) and c 44 (q, q z ) remain finite and can be reinterpreted as the elastic modulii of a viscoelastic liquid of an axially symmetric phase of elastic lines. The three-dimensional elastic modulii depend on material parameters and may be complicated functions of field and temperature due to the several microscopic length-scales that come into play [11], as discussed in Sec.IV. For the following discussion, it is only relevant to consider that c 11 (q, q z ) and c 44 (q, q z ) remain finite and axially-symmetric in the q → 0, q z → 0 limit of quasi-uniform compressions and tilting deformations. Strikingly, we will show that, at thermal equilibrium, every constant-z cross section becomes incompressible, even though the three-dimensional system is compressible.
The starting point for the elastic line hydrodynamics of line liquids is the Landau free energy functional, where r = (x, y), c 44 (r − r , z − z ) and c 11 (r − r , z − z ) are the non-local tilt and compression modulii [12] with Fourier transforms c 11 (q, q z ) and c 44 (q, q z ), the two-dimensional vortex density fluctuations at layer z around its mean value n 0 and the two-dimensional tangent field density for a collection of N 1 vortex-lines positioned at r j (z) at a given constant-z cross section. We assume a slab geometry with a thickness L in the z-direction and an area A in the x − y plane such that n 0 = N/A ≡ B/Φ 0 (see Fig.1). The last term in F describes the coupling of the vortex density with the pinning potential V D (r, z), which can have different correlations. In this paper we will be interested in the cases of uncorrelated (strictly speaking short-range correlated) isotropic point disorder and the long-range correlated disorder associated to CD.
Density fluctuations around the average can be measured in the reciprocal space by δn (q, z) = n (q, z) − n 0 Aδ q,0 where q and q z are the wave-vectors in the inplane directions and in the average-z direction of the lines. Fluctuations are conveniently quantified by the full three-dimensional structure factor where . . . and . . . denote averages over disorder and thermal fluctuations, respectively. These averages can be calculated, at equilibrium, from Eq.(5) with the constraint imposing the continuity of vortex lines [7]. It is usually agreed [7,13,14] that the elastic modulii of the liquid phase c 11 and c 44 should be similar to those of the solid or glassy phase, an assumption we will follow in general. We will also neglect the small renormalization of the modulii in the presence of uncorrelated disorder [15]. However, for correlated disorder generated by CD, c 44 is drastically renormalized and diverges at the Bose glass transition due to the broken statistical symmetries. Moreover, also c 11 is expected to get strongly renormalized and to diverge in the so-called Mott glass phase for the field such that the number of vortices equals that of CD (matching field). To be aware of these important differences between these two types of disorder, we will denote the elastic constants in the case of correlated disorder asc 44 andc 11 .
Density fluctuations are essentially controlled by compression modes [14]. Therefore, although the free energy of Eq.(5) is presented for the liquid phase (c 66 = 0), in the ideal clean case (V D = 0) this expression can be applied to compute the density fluctuations of the solid phases as well. In real cases with disorder, for the solid phases the situation is more subtle since the coupling to the pinning potential strongly depends on the periodicity of the vortex structure a 0 . This implies that the c 66 , not included in the description of Eq.(5), does play a role. Nevertheless, within the elastic approximation, hydrodynamic density fluctuations in the solid phase with disorder can be directly related to the prediction made for the displacement field (see discussion in Sec. II B 3).
The density fluctuations at a constant-z cross section S(q) of Eq.1 can be written as n 0 S (q) = |δn (q, z)| 2 (10) assuming a bulk system with statistical invariance along the z-direction. S(q) can thus be computed from S 3d (q, q z ) by integration over q z , where s sets an ultraviolet cut-off coming from the discretization in the z-direction or, more physically, from the superconducting layer spacing [10].

Liquid and Solid phases without disorder
We start by reviewing the ideal case without disorder at thermal equilibrium, V D = 0 in Eq. 5. In the absence of anisotropies in the x − y plane, we expect to have roughly the same three-dimensional structure factor for the liquid, S 3d liq (q, q z ), and the solid Abrikosov lattice, S 3d sol (q, q z ), vortex phases at small enough wave-vectors, In this expression obtained from Eq.(5) [7,10,14,16], k B is the Boltzmann constant and T the bath temperature. We are again assuming that the elastic constants in the vortex liquid can be well approximated by those of the solid phase at the hydrodynamic scales (q·a 0 < 1). Eq. (12) can be obtained by relating the spontaneous density fluctuations measured by S 3d liq (q, q z ) with the associated linear response function n 0 q 2 [q 2 c 11 (q, q z )+ q 2 z c 44 (q, q z )] −1 via the static fluctuation-dissipation theorem. Therefore, the most important difference between the response of the two phases, a finite c 66 (q, q z ) in the solid phase, does not play any role here. We also note that at small wave-vectors compared to the wave-vector associated to the Bragg peaks, q 0 = 2π/a 0 , the solid phase has isotropic density fluctuations since the broken rotational and translational symmetries of this phase are irrelevant.
Neglecting for the moment the dispersivity along the z-direction (see Sec. IV for a specific discussion for our experimental case), and approximating c 11 and c 44 by their values at q z = 0, we obtain following Ref. 10 for a constant-z cross section of the liquid and solid vortex phases in a clean sample (i.e. V D = 0). If we assume constant c 11 (q, 0) and c 44 (q, 0) as q → 0 we find that Therefore, according to the classification of Ref. 2, in the ideal clean case the point pattern at any constant-z cross section is a class-II hyperuniform system (α = 1) for the liquid and solid Abrikosov lattice phases. Some examples of systems belonging to this hyperuniform universality class are quasicrystals, classical disordered ground states, zeros of the Riemann zeta function, eigenvalues of random matrices, fermionic point processes, superfluid helium, maximally random jammed packings, perturbed lattices, perfect glasses and density fluctuations in the early Universe. [2].

Weak uncorrelated disorder: Liquid phase
The presence of disorder introduces corrections to the prediction of Eq. (13). Weak bulk point disorder, characterized by the correlator V D (r, z)V D (r , z ) = ∆δ(r − r)δ(z − z ), produces an additive correction to the vortex liquid three-dimensional structure factor (Eq. 12) such that with ∆ measuring the strength of the bulk disorder and where we have again approximated the non-negligible elastic constants of the liquid by those of the Abrikosov lattice. This calculation neglects the dispersivity of the elastic constants along the z-direction by taking their values for q z = 0. The correction to the two-dimensional structure factor is thus obtained by integrating over q z , If we again assume that c 11 (q, 0) and c 44 (q, 0) tend to constants as q → 0, we get as also found for the clean case. Note however, that the prefactor is different: Only when n 0 ∆/c 11 (q, 0) < k B T , the system crossovers to the line-liquid in clean samples. This case of α = 1 or class-II hyperuniformity is different to the α = 2 class-I hyperuniformity that is found for instance for pancake vortices interacting logarithmically in a thin-film superconductor or the two-dimensional onecomponent plasma system. The two-dimensional density fluctuations of a vortex structure in a constant-z cross section of a thick superconductor, described by Eq. (17), are instead equivalent to the ones expected for an effective two-dimensional particle system interacting with the long-range Coulomb ∼ 1/r repulsion. Such case is indeed known to have a two-dimensional compression modulus diverging as ∼ 1/q, which complies with the fluctuationcompressibility theorem. Hyperuniformity in the liquid phase is rather robust since arises from the continuity of flux lines, see Eq. (9). This constraint implies that we need to deform vortex lines along the z-direction in order to compress point vortices at a given cross section, inducing a divergent effective two-dimensional compression modulus at equilibrium. It is thus an interesting example of hyperuniformity emerging in the three-dimensional density-contour levels of an extended system composed by objects with short-range repulsive interactions.
Summarizing, for constant-z cross sections the vortex liquid phase is expected to be an α = 1 class-II hyperuniform system, either in the absence or presence of weak uncorrelated disorder.

The Bragg glass phase
As mentioned, the Abrikosov lattice at thermal equilibrium is predicted to be a class-II hyperuniform system at constant-z cross sections. Disorder-driven corrections to the two-dimensional structure factor in the solid or glassy phases are subtle, and must be considered. For weak uncorrelated disorder the corrections in the solid phases are quite different than in the liquid. This is due to the strong relevance of periodicity in the coupling of the vortex density with disorder [4,17]. Fortunately enough, if disorder is weak and the temperature low, the elastic theory can be applied and the vortex structure can be described by a hydrodynamic two-component displacement field u α (r) measuring the distance of a vortex with respect to the corresponding perfect lattice position, such that u α (R i , z) = u i,α (z) is the displacement of the i-th vortex line at the slice z. In the case of weak uncorrelated disorder the Bragg glass phase with quasi long-range positional order is predicted [4,17,18]. The periodic coarse-grained vortex density describing the topologically ordered phase can be expressed as where the ". . . " denote rapidly oscillating terms that do not contribute to the hydrodynamic density modes we are interested in. The hydrodynamic three-dimensional structure factor is where u L (q, q z ) is the component of the displacement field in the direction parallel to q. This expression is rather general, see Ref. 19. If we consider the Abrikosov lattice corresponding to the clean system below the melting temperature, the thermal roughening of such displacement component is The structure factor of the lattice is thus essentially the same as for the liquid, since the structure factor in both phases is controlled by the longitudinal modes only. At finite temperatures, the Abrikosov lattice at constant-z hence displays the same α = 1 hyperuniformity as the liquid, as already discussed in the previous section. Weak disorder destroys the perfect long-range positional order, but not the topological order, and a particular quasi-long range order emerges. In this Bragg glass phase [4,17] the displacement field grows as a power-law at intermediate distances, and logarithmically at large distances, in contrast with the temperaturedependent saturation of displacements found in the Abrikosov lattice (associated to the Debye-Waller factor). The disorder-induced roughening of the displacement field is expected to be isotropic in the x − y plane at large-distances with its corresponding structure factor, see Eq. 19, to reduce to the clean thermal case of Eq.(12) on vanishing disorder. Heuristically, we can thus approximate the thermal-and disorder-induced displacements as with d = 3 the space dimension and ζ a characteristic exponent. Note that we have used the same c 44 and c 11 as for the clean case, assuming no disorder-induced renormalization. In presence of thermal noise we have ζ = 1 − d/2, thus reducing to Eq. (12). When ζ > 0 this exponent is the so-called roughness exponent of the lattice, such that the mean square displacement grows as with R the observation scale or system linear size in the x−y plane. For ζ = 0 the lattice is logarithmically rough u 2 L ∼ log R, while for ζ < 0 has a macroscopically flat displacement field.
By integrating over q z and assuming consistently that the elastic constants c 11 and c 44 tend to finite values at small q z , we obtain the constant-z cross section two-dimensional structure factor of the pinned solid, Therefore, under these assumptions, the elastic system is hyperuniform only if ζ < 0 (displacement field macroscopically flat) while for ζ = 0 is marginally hyperuniform, i.e. α = 0. This is particularly important for the Bragg glass and for systems with quasi-long-range order (such as two-dimensional crystals with short-range interactions), since ζ → 0 asymptotically. This so-called random-periodic regime (RP) of the Bragg glass phase exists for qR a 1, where R a is the scale at which the displacement field has fluctuations of order a 0 . For R c < R < R a , where R c is the Larkin radius on the x − y plane, there is a crossover to the so-called randommanifold regime (RM) where ζ ≡ ζ RM ≈ 0.2. For scales shorter than R c , the system crossovers to the Larkin regime, where ζ ≡ ζ L = (4 − d)/2 = 1/2. In none of these regimes the Bragg glass phase presents suppressed density fluctuations.
Summarizing, any arbitrarily weak disorder kills the class-II hyperuniformity of the Abrikosov lattice, and the Bragg glass phase is expected to be marginally hyperuniform at constant-z cross sections. Rather surprisingly, from the given arguments, the Bragg glass is expected to jump from marginal to a class-II hyperuniformity at melting, adding a new signature to this first-order phase transition [4,17].

CD correlated disorder: Liquid and Bose glass
The case of the correlated disorder generated by randomly-distributed CD, such that V D (r, z)V D (r , z ) = ∆ 1 δ(r − r), is special since the long-range correlation along z favors the localization of vortex lines into CD. A Bose glass transition is expected by lowering the temperature from the liquid phase [13], with the concomitant divergence ofc 44 approaching the transition temperature, whilec 11 remains finite. There is an exception for the putative Mott glass phase at the matching field B = B Φ ≡ n col Φ 0 wherec 11 is also expected to diverge. The predicted three-dimensional structure factor for the vortex liquid in presence of CD disorder is , wherec 11 andc 44 are the compression and tilt modulii in the presence of correlated CD disorder, and ∆ 1 measures the strength of the columnar disorder. Integrating over q z we obtain the two dimensional structure factor This expression yields a crossover wave-vector such that for scales q < q CD , the density fluctuations in the liquid are dominated by the pinning introduced by CD (second term in Eq. 27), while for q > q CD they are similar to those found in the liquid phase without correlated CD disorder (first term in Eq. 27).
On approaching the Bose glass transition on cooling, c 44 increases rapidly and diverges at the transition, con-comitantly with q CD → ∞. We hence expect . (29) in the Bose Glass phase where the localized vortex lines repel each other with short-range interactions. The compression modulus remains finite at the transition [13] and therefore a non-hyperuniform system is expected at equilibrium. Summarizing, correlated disorder generated by CD destroys hyperuniformity at constant-z cross sections, both in the Bose glass and in the liquid phases. However, since on increasing temperature q CD vanishes, a crossover to a nearly class-II hyperuniform liquid can occur within the liquid phase. This crossover can also occur on decreasing the magnitude and/or density of correlated disorder ∆ 1 .

Finite-size effects
In the previous subsections, we assumed infinite samples and showed that the ideal clean liquid and solid phases, as well as the liquid phase with weak disorder, are class-II hyperuniform systems. In samples with thickness L, finite effects in the z-direction will affect hyperuniform behavior inducing a crossover at a characteristic q F S (L). Indeed, by a simple dimensional analysis of the threedimensional structure factors in these hyperuniform cases (Eqs. 12 and 15), we get Physically, q F S is the characteristic wave-vector at which the correlation length of vortices in the z-direction becomes of order L [10]. In order to illustrate this, for simplicity we perform the calculation from Eq. 12 for the clean system. By Fourier inverting in the z-direction, with ξ (q) = q −1 c 44 (q, 2π/L)/c 11 (q, 2π/L).
Therefore, Eq. 30 is equivalent to the physical condition ξ (q F S ) = L. For q < q F S the three-dimensional system essentially behaves as a two-dimensional system of rigid lines with the structure factor obtained by setting z 1 = z 2 , and ξ = L, Provided L > λ ab at the crossover, rigid vortices repel with a short-range interaction yielding [14,15] c 11 (q) = 2 0 2πλ 2 with 0 = (Φ 0 /4πλ ab ) 2 the line tension. Therefore, since the length-scale in the x − y plane probed is q −1 λ ab , we have S(q) → const as q → 0. In other words, finitesize effects kill the class-II hyperuniformity of the threedimensional vortex liquid or solid phases in clean samples at thermal equilibrium. Although this result was obtained for a particular case, this finite-size effect is expected to be present in all vortex phases, even with disorder. Indeed, the existence of q F S (L) comes from first, a dimensional analysis of the competition between tilting and compression elastic responses, and second, from the non-dispersive behavior of the three-dimensional elastic constants.

Two-dimensional systems
It is interesting to discuss Eq. 33 in the limit of very thin superconductors such that L λ ab . In this limit the effective penetration length becomes Λ = 2λ 2 ab /L. If the observation length-scale in the x − y plane is qΛ > 1, two different power-law behaviors can be expected. At very small length-scales, qΛ 1, interactions are logarithmic and in Eq. 33, so the system at equilibrium is class-I hyperuniform with α = 2, for log(1/x) interactions. Recent simulations with this type of interaction show that the value α = 2 is robust under the presence of quenched uncorrelated disorder. [5] For larger length-scales, qΛ > 1, the interaction becomes Coulomb-like and yielding, at equilibrium, a class-II hyperuniform vortex pattern with α = 1 and thus for 1/x interactions. Interestingly, this is the constant-z cross section behavior predicted for the three-dimensional line array. In other words, the effective two-dimensional interaction at a constant-z cross section, mediated by the short-range three-dimensional interaction between vortex lines, mimics Coulomb interactions (as in Wigner crystals). This is analogous to the long-range elasticity of the triple contact line of a liquid meniscus, so called "fringe elasticity of the line of contact", reflecting the energetics of deformations of the liquid/gas interface [20,21]. Although we are assuming clean systems, since for the two-dimensional elastic lattice in the presence of weak disorder c 11 is not renormalized [15,22], we can expect the result α = 1 to hold even in the presence of weak uncorrelated pinning. This would be the case in the solid as well as in the liquid two-dimensional vortex structure, due to their similar uniaxial compression properties at small wave-vectors.

Summary of analytical predictions
In Table I we summarize our predictions regarding the type of hyperuniformity of different vortex phases at constant-z cross sections of a three-dimensional superconductor (disregarding the finite-size effects discussed above). Liquid vortex phases are expected to be class-II hyperuniform for weak uncorrelated pinning, while correlated disorder destroys hyperuniformity. As discussed, a crossover of the liquid with CD correlated disorder towards an apparent class-II hyperuniformity at large temperatures or weak disorder is possible. At low temperatures, only the Abrikosov lattice with thermal fluctuations is class-II hyperuniform. The Bragg glass phase, expected for weak uncorrelated disorder, is marginally hyperuniform due to the disorder-induced roughening of the displacements field, while the Bose glass phase is not hyperuniform due to the enhanced correlations along the z-axis. The putative Mott glass phase is expected to be hyperuniform due to the divergence of both, the tilt and the compression modulii. To the best of our knowledge, the wayc 11 diverges with q has not been reported for the Mott glass, and hence we can not guess its hyperuniformity class yet. All these predictions assume three-dimensional vortex line phases, non-dispersive elastic constants in the low q limit, and neglect finite-size effects.

C. Numerical Simulations
The hydrodynamic predictions presented in the previous sections are based on a coarse-grained continuous model and the assumption of thermal equilibrium. In order to test these analytical predictions, and to match two different levels of description of the problem, we have performed molecular dynamic simulations of a simple microscopic model. This helps us to understand finite-size and discretization effects in the z-direction, as well as the effect of the particular shape of the vortex-vortex interaction potential, without relying in the approximated elastic modulii of the continuum description.
We model vortices as elastic lines discretized in the z-direction, such that r i (z) ≡ (x i (z), y i (z)) describe their two-dimensional coordinates at the layer z, with z = 0, ..., N z − 1 and N z the total number of layers. Periodic boundary conditions are considered in all directions. The total energy of the elastic lines array is E[{r i (z)}] = E 1 + E int such that each line has an elastic tension energy given by Hook coupling with strength k, and the repulsive interaction energy between vortex-lines is modeled as with K 0 (x) the zero-order modified Bessel function of the second kind. Eq. 40 is derived from the London model, while Eq. 39 with k ∝ 0 is a harmonic local approximation for the single-vortex elastic tension. In layered superconductors as the one experimentally study here, the elastic tension energy arises from the attractive electromagnetic and Josephson couplings between pancakes of different superconducting layers. We consider an overdamped Langevin dynamics at a temperature T where η is the Bardeen-Stephen friction. At long enough times this system equilibrates in the canonical ensemble at temperature T . In order to approach experimental conditions we have simulated a range of a 0 close to the typically accessed in magnetic decorations. Figure 2 shows typical two-dimensional structure factors S(q) averaged over the N z layers. The liquid and solid without disorder in thermal equilibrium, have both an isotropic structure factor at low wave-vectors, although at large wave-vectors the solid displays the peaks corresponding to the triangular Abrikosov lattice. This isotropy at low q is consistent with the hydrodynamic assumption made in Sec. II B considering axiallysymmetric elastic constants. Thus, for the q-range we work with the angular average of the structure factor, S(q) ≡ (2π) −1 2π 0 dΨS(q cos(Ψ), q sin(Ψ)), in order to examine large length-scale density fluctuations.
The inserts of Fig. 3 show S(q) for a line liquid above the melting temperature and for an ordered line solid, both at equilibrium and without disorder. in the z-direction. There are clear finite size effects, i.e. N z dependence, in the z-direction, and they are amplified when a 0 increases. We find that these size effects can be quantified by appropriately normalizing the inset curves into a master curve: we consider a characteristic finite-size crossover scale q F S such that S(q)/S(q F S ) ∼ G(q/q F S ), with G(x) ∼ 1 for x 1 and a unique G(x) for x 1. For the range of a 0 analyzed, we find a good collapse in the low-q region using the crossover wavevector q F S = a 0 /N z , and S(q F S ) ∝ T a 2 0 /L, both for the liquid and solid phases. Therefore, the saturation of S(q) for q < q F S is actually a finite-size effect, and it should disappear in the thermodynamic N z → ∞ limit. The dependence of q F S with N z ∝ L confirms the size effects predicted in Sec. II B 5. In addittion, the shape of the master curve for q > q F S , G(x) ∼ x, confirms an α ≈ 1 typical of class-II hyperuniformity (see Eq 13).
The evolution of q F S with a 0 is less universal than the dependence with N z since relies on the precise functionality of the elastic modulii c 44 and c 11 that depend on the type of inter-vortex interaction and the field regime explored in the simulations. However, if we use the rigid vortex compression modulus of Eq. 34, c 11 (q, q z = 0) ∝ B 2 /(1 + q 2 λ 2 ab ), the single-vortex tilt modulus c 44 ∝ 0 ∼ B 2 a 2 0 (a constant in this particular model),

FIG. 3.
Normalized two-dimensional structure factor obtained from Molecular Dynamics simulations of a threedimensional elastic-line array in the ideal clean case, for various lattice spacings a0 and number of layers Nz. We use a finite-size crossover scale qF S ∝ a0/Nz and a normalization S(qF S ) ∝ T a 2 0 /Nz. (a) Liquid phase at equilibrium and (b) solid Abrikosov phase. Both phases display the same finitesize crossover at qF S . For q < qF S the system effectively becomes a non-hyperuniform two-dimensional system, while for q > qF S it shows a fair S(q) ∼ q hyperuniform behaviour at low q. Upward corrections at large wave-vectors, qλ ab 1, can be explained by the dispersivity of the compression modulus emerging from the particular interaction potential. Lines are a guide to the eye for S(q) ∼ q α with α = 2 (dashed line) and α = 1 (dashed-dotted line). The insert shows nonnormalized data. and approximate Eq. (30) as q F S ≈ c 44 /c 11 (0, 0)/L, we get q F S ∝ a 0 /L. Also, from Eq. 33 we get S(q F S ) ∼ n 0 k B T /c 11 L ∝ T a 2 0 /N z . Both results are in agreement with the outcome of simulations shown in Fig.3 confirming that the c 11 for rigid vortices is adequate for describing the S(q) saturation for q < q F S .
For q > q F S the master curve presents deviations from S(q) ∼ q for q q F S , crossing over to a more rapid increase of S(q) with q. This behavior can be also qualitatively explained in terms of the hydrodynamic prediction of Eq. 13 if we consider that S(q) ∝ q/ c 44 c 11 (q, 0) ∼ q 1 + λ 2 ab q 2 ≈ q(1 + q 2 λ 2 ab /2 + . . . ). Therefore, this microscopic model confirms the existence of finite-size ef-fects, and the class-II hyperuniformity of the liquid and solid vortex-line phases in ideal clean samples. We also explain the upward deviations from the α = 1 behavior observed at q > q F S controlled by the dispersive behavior of c 11 , emerging when qλ ab ∼ 1. As discussed later, dispersive effects may be important for the range of parameters given by the magnetic decoration experiments in the layered superconductor Bi 2 Sr 2 CaCu 2 O 8+δ . The study presented in this section bridges the hydrodynamic scale with the molecular dynamics scale, using a simple elastic-line array model.

III. DENSITY FLUCTUATIONS: EXPERIMENTAL RESULTS
In order to test the previous theoretical predictions we have performed magnetic decorations in various superconducting samples giving us direct access to twodimensional point patterns at the top surface. This layer is the particular z = L cross section where surface interactions may play an important role. However, it has been argued that the structure at the surface is representative of that at the bulk of the sample for the typical in-plane length-scales probed by decorations ∼ a 0 [10]. An alternative experimental method to the one described here would be to study small-angle neutron scattering data, yielding directly the structure factor. We followed this path but found that the experimental resolution of this technique at low q is not good enough as to ascertain whether hyperuniformity is present in the vortex system.
The choice of the extremely layered Bi 2 Sr 2 CaCu 2 O 8+δ system is based in the easiness to obtain freshly cleaved surfaces to perform several decoration experiments in the same crystal, and on the availability of samples with different types of disorder (pinning potentials). We studied vortex lattices nucleated on a large set of singlecrystals with the natural disorder coming from crystalline defects (pristine samples), and with additional disorder introduced by irradiation with electron and heavy ions. Electron-irradiated samples present a dense pointlike disorder whereas heavy ion-irradiated ones present a Poisson-like distribution of CD, columns of crystallographic defects aligned along the z-direction. Electron irradiation was performed with 2.3 MeV accelerated electrons in a van de Graaff accelerator coupled to a closedcycle hydrogen liquifier at theÉcole Polytechnique of Palaiseau, France [23]. This process is performed at low temperatures (20 K) as to guarantee the stability of Frenkel pairs created in the irradiation process. The data presented here corresponds to a sample irradiated with an electron density of 1.7 10 19 e/cm 2 [23].  We image individual vortex positions in a typical fieldof-view of thousands of vortices by performing magnetic decoration experiments at 4.2 K after field cooling at different applied fields in the range 5 < H < 150 Oe. Vortices are decorated with Fe particles attracted by the local field gradient generated around the vortex cores [24], observed as black dots in the inverted scanning-electronmicroscopy images of Fig. 4. Further details in the field-  Fig. 4. The intensity is shown in a logarithmic scale and the color level is the same for the three images, see color bar at the bottom. The gray crosses indicate the q-window affected by spurious effects arising from the field-of-view edges. Data in these crosses are not considered for the calculation of the angularly-averaged structure factor S(q).
cooling decoration protocol followed in this case can be found in Ref. 24. For every studied sample, several (∼ 10) magnetic decoration experiments were performed in freshly cleaved surfaces, eventually at different applied fields. We studied ∼ 30 pristine, 1 electron-irradiated, and 10 heavy-ion irradiated single crystals. Figure 4 show examples of magnetic decoration images at various applied fields in the three types of studied samples. Panels (a) and (b) correspond to snapshots of vortex positions taken at a field well within the quasi-crystalline Bragg glass phase for pristine and electron-irradiated samples. For these two types of samples, at fields larger than 15 Gauss, the vortex structure is single-crystalline, presents quasi-long-range positional order and very few topological defects associated with non-sixfold coordinated vortices. This is observed in the Delaunay triangulations superimposed to the pictures with blue lines joining first neighbors. Topological defects are highlighted with red dots. For both images only ∼ 2 % of vortices are involved in defects, mainly edge dislocations. For fields smaller than 15 Gauss, the structure breaks into small crystallites for pristine [24] as well as electron-irradiated samples. This polycrystalline structure results from vortex-vortex interaction weakening and disorder becoming more relevant on the viscous freezing dynamics [24]. Further details on the field-evolution of the structural properties (images, density of defects, S(q), displacement correlator and correlation lengths) of the Bragg glass phase in pristine and electron-irradiated samples can be found in our previous work Ref. 25.
The structural properties of the vortex matter nucleated in samples with CD are qualitatively different than for the case of point disorder: small misaligned crystallites with less than 20 vortices are observed at densities up to 100 Gauss, see Fig. 4, and for fields smaller than ∼ 15 Gauss the structure is amorphous. In the latter case the density of topological defects can be as high as 50 %, whereas at high fields the amount of non-sixfold coordinated vortices is always larger than 40 %. The nucleation of these structures with at best 10 lattice spacing shortrange positional order comes from the strong effect of pinning introduced by CD that are randomly distributed in the sample [26]. Nevertheless, the positional correlation of these vortex structures is not that of a random distribution of points mimicking the Poisson-like distribution of CD: the pair correlation function of the vortex structure presents one very good defined principal peak at r = a 0 and secondary weak peaks in some cases.
The decorated structures were frozen, at length scales of order a 0 , at temperatures at which the pinning generated by disorder sets in [27,28]. This freezing temperature is some K below the temperature at which magnetic response becomes irreversible [29]. This lenghtscale is shorter than the relevant for the hyperuniformity analysis. Due to their different relaxation rates, density modes corresponding to length-scales larger than a 0 are expected to freeze at higher temperatures than the modes with q = 2π/a 0 . We will further discuss this issue later.
The images of Fig. 4 are just some examples, but typically we study a panoramic images of the vortex structure with no less than 1500 and up to 15000 vortices. We digitalize vortex positions in these large fields-of-view and then we calculate the structure factor S(q) and the number variance σ 2 N . The larger the field-of-view, the smaller the value of q ≡ q 2 x + q 2 y that we can access to calculate the two-dimensional structure factor at the top surface. Figure 5 show the corresponding two-dimensional struc-FIG. 6. Angular average of the structure factor for the magnetically decorated vortex structure at the surface of (a) pristine, (b) electron-irradiated, and (c) Xe-irradiated (correlated CD disorder) Bi2Sr2CaCu2O 8+δ samples. The vortex density in every case is indicated. Data are shown as a function of q/q0 with q0 = 2π/a0 the Bragg wave-vector. The red line is the best power-law fit S(q) ∼ q α eff for the low wave-vector range a0/6 < q/q0 < a0/2.  Fig. 4. The S(q x , q y ) values are angularly averaged as to obtain the S(q) data shown in fig. 6.
We would like to point out some important technical difficulties that are quite specific to the study of the lowq density modes. On one hand, the borders and shape of our field-of-view hinders the study of the structure factor in the relevant low-q range due to the annoying window-ing effect. In rectangular fields-of-view the artifact shows up as an excess in S(q x , q y ) localized in a "+"-shaped region centered at q x = q y = 0. This artifact, associated to the Fourier transform of the field-of-view, is oriented along the principal directions of the rectangle and has an oscillatory decay on increasing q. This border effect is avoided in simulations with in-plane periodic boundary conditions. To get rid of this spurious effect, for the smallest q we perform a partial average over Ψ-values outside the '+"-shaped region, instead of averaging over all Ψ. This is possible due to the high angular localization of the artifact. Nevertheless, since the cross has a finite width, at the end we are forced to consider a safety minimal wave-vector q min > 2π/W f ov , with W f ov a characteristic linear size of the field-of-view. On the other hand, as we do not have an ensemble of many magneticallydecorated vortex configurations at exactly the same position of the samples to average over, inevitably, the lower the q, the bigger the statistical fluctuations in S(q). Henceforth, in practice, although our images span distances as large as W f ov ∼ 50a 0 , we end up analyzing density fluctuations with q > q min q 0 /6 = 1/6(2π/a 0 ).
Typical results for S(q) are shown in Fig. 6 (top panels) for the pristine (a), electron-irradiated (b), and heavyion irradiated (c) samples. The dashed and dasheddotted lines are guides to the eye, showing q 2 and q evolutions, respectively. For the lowest q/q 0 the behavior is consistent with class-II uniformity with α = 1 for the three samples. Indeed, red lines in the figure shows a fit to the data yielding α eff ∼ 1. Another way to analyze the density fluctuations is to study the distance-evolution of the vortex-number variance σ 2 N , (defined in Eq. 4). To this end, we considered circular regions with radius R and centers located at random within the panoramic image of the vortex structure. We pay attention to the circles not crossing or touching the edge of the field-of-view. We make statistics on the number of vortices contained over a large amount of circles of size R and the number variance σ 2 N (R) is then computed. Fig. 6 (bottom panels) shows the results of the R-evolution of σ 2 N corresponding exactly to the structures considered in the top panels. We observe a rough agreement with S(q), in the sense that a class II hyperuniform scaling σ 2 N ∼ R/a 0 fairly describes the data (see the dashed-dotted line). Indeed, the red lines in Figs. 6 (d), (e) and (f) correspond to functions ∼ R 2−α eff (see Eq. 4) with the α eff obtained from the fits of the S(q) ∼ q α eff data of the top panels.
In order to perform a comprehensive study of the occurrence of hyperuniformity in vortex matter with different types of disorder and for different vortex densities, we have systematically fitted S(q) ∼ q α eff for all a our studied samples and vortex densities (magnetic field). The effective power-law exponent α eff is obtained by fitting in the range q/q 0 ∈ [1/2, 1/6]. The effective exponents as a function of field for every type of disorder are shown in the three panels of Fig. 7, compiling data from over roughly 40 statistically independent cases. Although S(q) and σ 2 N qualitatively agree, we have found that a systematic fit of α eff using the expected scaling of σ 2 N is empirically more difficult due to the strong nonasymptotic corrections to the number variance [2]. As we can observe, α eff ≈ 1 ± 0.3 is rather robust for all the studied cases, independently of the type of disorder, correlated or uncorrelated, present in the samples.

IV. DISCUSSION AND PERSPECTIVES
A strict hyperuniformity analysis requires data in an asymptotic regime which is rather difficult to reach and assure experimentally. Nevertheless, our experimental results in an extensive data-set display a clear suppression of the amplitude of the density fluctuations with S(q) ∼ 10 −2 for the lowest q accessed (to be compared with S(q) = 1 for the Poisson or ideal gas particle distribution). Hence, the vortex structures frozen during the decoration field-cooling protocol are nearly hyperuniform two-dimensional point patterns at the superconductor surface. Furthermore, these structures display S(q) ∼ q α eff with α eff ≈ 1 ± 0.3 in the low-q range accessed in our experiments. This agrees with the equilibrium hydrodynamic predictions for the large-scale fluctuations of the liquid phase with weak uncorrelated disorder, as well as for the solid and liquid vortex phases in the ideal clean case. However, this exponent contrasts with the predictions for the Bragg glass phase (α = 0) and for the liquid and Bose glass phases with CD (S(q) → const as q → 0). We will argue that these discrepancies can be explained by considering the unavoidable relevance of memory effects during the decoration field-cooling protocol that affect the observed vortex structures at different length-scales. These non equilibrium effects are ignored in the hydrodynamic description of Sec. II. We will also argue that dispersion effects and finite-size effects are experimentally relevant for a systematic study of density fluctuations in vortex matter.
The two-dimensional point pattern obtained at the surface of the sample after a field-cooling magnetic decoration down to 4.2 K corresponds to an out-of-equilibrium structure. This is due to the strong dependence of the relaxation rate of the vortex structure with the wavelength of the density modes. In field-cooling experiments with a fixed cooling rate, large length-scale density fluctuations have a slower relaxation rate than smalllength-scale ones, and then the former fall out of equilibrium at larger temperatures. This out-of-equilibrium effect is expected to be enhanced by disorder that dramatically slows down the thermally-activated dynamics. In other words, a magnetic decoration image of the vortex structure at low temperatures is not just a snapshot but also a photograph with memory of its field cooling history. Then, the freezing temperature of the vortex structure is not a unique magnitude but rather depends on the length-scale of the density fluctuations, T freez = T freez (q). In a field-cooling experiment, at T ∼ T freez (q) the density fluctuations of mode q have a relaxation rate of the order of the experimental cooling rate. Namely, at T < T freez (q) all modes with wavelength ≤ 1/q are out of equilibrium.
In particular, the local ordering at the scale a 0 , observed in samples with weak uncorrelated disorder indicates that T freez (q ∼ a −1 0 ) ≤ T m , with T m the melting temperature of the Bragg glass phase. Density fluctuations associated to larger wavelengths are then expected to fall out of equilibrium at much higher temperatures, T freez (q a −1 0 ) > T m . These slow modes with q → 0 retain memory of the liquid phase. This can explain our experimental observation of a near class-II hyperuniform vortex structure in pristine and electron-irradiated samples (predicted in the equilibrium hydrodynamic approximation for the liquid phase with weak disorder). This argument can be further justified in the framework of disordered elastic systems without topological defects. Indeed, even a simple elastic string relaxing in a random medium after a temperature quench displays a logarithmically growing correlation length [14,30] separating equilibrated from non-equilibrated length-scales. Glassy vortex dynamics thus prevents the system to reach the marginal hyperuniformity predicted for the Bragg glass phase at equilibrium. This argument also explains why it has been in general difficult to observe the different crossover regimes predicted for the Bragg glass in magnetic decoration experiments. Indeed, decoration experiments have been so far only capable of revealing the random manifold, but not the random-periodic regime of lattice roughness in samples with weak uncorrelated disorder [29,31].
The above out-of-equilibrium qualitative explanation is more subtle to apply for the samples with CDs correlated disorder since in this case the line-liquid as well as the Bose glass are expected to be non-hyperuniform at equilibrium in an hydrodynamic approximation. How is then possible to obtain α eff ≈ 1 from the memory of a liquid with correlated disorder? A possible answer can come from the crossover of density fluctuations at q CD from a class-II hyperuniform liquid to a non-hyperuniform liquid on decreasing q. Observing the equilibrium structure at q → 0 implies measuring at very large length scales ∼ 1/q CD that might be hard to access experimentally in finite fields-of-view. To evaluate this hypothesis we need to compute where we have used in Eq. 28 that ∆ 1 ≈ U 2 0 b 4 0 /d 2 1 + O b 2 0 /d 2 is the disorder correlator strength for CD with radius b 0 and mean separation d b 0 . It is difficult to make a quantitative assessment of q CD due to the many microscopic parameters involved. However, correlated disorder by CD is expected to slow down the dynamics with respect to the weak disorder case, and then to increase T freez (q) for all q. Since q CD ∝ 1/T freez (q CD ), then q CD can become too small to be experimentally accessed for our typical fields-of-view. Thus, the observed structures effectively display a class-II hyperuniformity due to the memory of the liquid phase, with the first term of Eq. 27 dominating for q > q CD at very high temperatures. A systematic study is desirable to check this hypothesis, for instance comparing different densities of CD. In particular, increasing ∆ 1 will produce an enhancement of q CD , and then the predicted saturation of S(q) due to the dominance of the correlated disorder term might be observable in the same field-of-view. Infinitely thick sample L/λ ab 1 for various a0/λ ab corresponding to different applied magnetic fields. At high fields, corresponding to the smallest a0λ ab values, a crossover from the type II hyperuniform scaling S(q) ∼ q to S(q) ∼ q 2 is observed at qλc ≈ 1 (dotted lines are a guide to the eye for the two regimes). At low fields, the crossover shows up at larger q. The largest value of a0/λ ab = 2 (black curve) roughly corresponds to the typical field of 40 G applied in magnetic decorations for Bi2Sr2CaCu2O 8+δ with λ ab (T freez ) = 0.4 µm. Bottom Panel: Finite-thickness effect in the structure factor S(q) for various L/λ ab ratios at a fixed field such that a0/λ ab = 2. This effect produces the saturation of S(q) in a region of small q whose extension depends on L/λ ab . The case of L/λ ab = 1000 (black points) corresponds to a Bi2Sr2CaCu2O 8+δ sample with a thickness of 400 µm, a typical value in real samples. The gray shaded area corresponds to the q · λ ab range accessed experimentally for our typical field-of-view in decorations (obtained considered again λ ab (T freez ) = 0.4 µm). Finite size effects in S(q) are not dominant for the experimental window accessed in our experiments in Bi2Sr2CaCu2O 8+δ . However, dispersivity effects do play a role in our experimental window.
Another effect that we have so far not considered in the discussion is the anisotropy-enhanced dispersivity of the elastic constants for the lowest wave-vectors accessed in our fields-of-view. This is particularly important in layered materials as Bi 2 Sr 2 CaCu 2 O 8+δ since dispersivity can induce, as discussed, a crossover from class-II hyperuniformity to a different behavior on increasing q. Since we access scales q · λ ab 1, in-plane dispersivity is not important. However, dispersion in the z-direction is important when q·λ c = Γq·λ ab ∼ 1. The anisotropy parameter of our samples is Γ ≈ 170, and then the condition q ·λ c ∼ 1 can easily be reached in our experimental fieldsof-view since λ c (T freez ) = Γλ ab (T freez ) ≈ 70 µm. For an anisotropic superconductor the compression modulus is [11] c 11 (q, q z ) = B 2 4π 1 + λ 2 c (q 2 + q 2 z ) (1 + λ 2 ab (q 2 + q 2 z )) (1 + λ 2 c q 2 + λ 2 ab q 2 z ) (44) while the tilt modulus reads with c 44 (q z ) the isolated-vortex contribution. This last term, important at low fields, is dispersive in q z for small wave-length fluctuations only, such that q z λ ab ∼ 1. This dispersivity is due to the electromagnetic coupling between pancake vortices [14]. For the long wave-length fluctuations we are interested in, we can take the constant non-dispersive limit for the isolated vortex contribution, If we assume L → ∞ and integrate Eq. 12 over q z using the expressions of Eqs. 44 and 45, we obtain the twodimensional structure factor S(q) including dispersion effects. The result for our Bi 2 Sr 2 CaCu 2 O 8+δ samples is shown in Fig. 8, for different vortex densities in the case of an infinite sample, see top panel.
The bottom panel of Fig. 8 shows S(q) for different sample thicknesses at a fix vortex density a 0 /λ ab (T freez ) = 2 corresponding to a magnetic induction of ∼ 40 G. We considered λ ab (T freez ) ∼ 0.4 µm at the irreversibility temperature T irr (B ∼ 40 Gauss). This is the temperature at which pinning sets in, and the vortex structure is frozen at length-scales of a 0 [29]. The gray area indicates the q · λ ab (T freez ) range used for fitting α eff in our experimental S(q) data.
In the case of an infinite sample, S(q) presents a crossover, from S(q) ∼ q to S(q) ∼ q 2 at larger q · λ ab . This large-q deviation from the asymptotic hyperuniform class-II behavior for clean samples is due to the dispersivity in the elastic constants. The q-location of the crossover decreases on increasing field. In particular, in the limit of high fields, where the isolated vortex contribution c 44 can be neglected, the crossover scale tends to λ c (T freez ), see top panel of Fig. 8. Then, in order to ascertain whether a structure is hyperuniform, it is crucial to obtain α eff by fitting S(q) data in a q-range located below this dispersivity-induced crossover.
In addition, this fit can neither be performed in the low-q limit where finite size effects due to the finite sample thickness L destroy hyperuniformity. The bottom panel of Fig. 8 shows that this finite size effect, observed as a saturation of S(q) at low q values, extends over a larger q-range on decreasing L. This was already discussed in Sec. II B 5, though neglecting the dispersivity and strong anisotropy of the elastic constants for vortex matter in Bi 2 Sr 2 CaCu 2 O 8+δ . The figure highlights a black curve corresponding to the prediction for our approximate experimental situation, namely vortex density ∼ 40 Gauss and L ∼ 400 µm. In our experiments, the saturation associated to the finite size effect is well below the q-range taken into account in the α eff fits (gray-shaded area). Therefore our hyperuniformity study in vortex matter from magnetic decoration images is not affected by this finite-size effect. However, for our experimental length-scales, the saturation should be clearly visible for samples with L/λ ab (T freez ) ≈ 10, namely for L ≈ 4 µm. This is a rather thin sample to be obtained just by cleaving a thicker crystal. Magnetic decoration data in 1 µm thick Bi 2 Sr 2 CaCu 2 O 8+δ micron-sized samples are available, but the field-of-view is of only a few hundred of vortices, making difficult to perform a careful hyperuniformity analysis of such vortex nanocrystals [29,32].
Going back to the relevance of dispersivity effects for our hyperuniformity analysis, the departure from the S(q) ∼ q behavior for our typical experimental situation, see black curve in Fig. 8, starts at the upper half of our experimental fitting window. Therefore, although the easily cleavable Bi 2 Sr 2 CaCu 2 O 8+δ system makes the experimental realization of this magnetic decoration study possible, the hyperuniformity analysis is messed up by the dispersivity effects inherent to its extremely-layered nature. Then the α eff values obtained in our study might be slightly-overestimated due to this effect.
In order to study up to what point our obtained α eff values are overestimated due to out-of-equilibrium and/or dispersivity effects, further experiments should be performed. First, out-of-equilibrium effects could be quantified by altering the field-cooling process, either by significantly slowing down the cooling rate, or by adding an in-plane dithering field allowing the vortex system to more efficiently relax towards the free energy minimum [33]. Second, dispersivity effects can be reduced by choosing a less anisotropic host superconducting material. In a more general perspective, complementary experimental studies of vortex density fluctuations by applying other different techniques, such as neutron scattering, would be of great interest. From the theoretical point of view, extending the hydrodynamic theory to the realistic non-equilibrium situation of field-cooling experiments, as well as performing non-equilibrium simulations under the same conditions, would be very insightful.
In summary, we have shown that a directed threedimensional elastic line array with short-ranged interactions in all directions can generate, at equilibrium, a twodimensional hyperuniform point pattern at their cross sections. This behavior can be explained by the threedimensional bulk energetics involved in the compression of a single layer, and ultimately in the continuity of lines.
We have shown that disorder can play an important role on modifying in several interesting ways such hyperuniform behavior. In order to test the above physical picture we have chosen the case-study system of vortex structures in type-II superconductors, a paradigmatic experimental realization of the three-dimensional elastic line array. We have studied this issue both, theoretically via analytic calculations and numerical simulations, and experimentally in the anisotropic vortex matter nucleated in Bi 2 Sr 2 CaCu 2 O 8+δ samples with correlated and uncorrelated disorder. By magnetically decorating the frozen two-dimensional vortex point pattern at the top surface of different superconducting samples we find a systematic suppression of vortex-density fluctuations at long wavelengths. We have argued that this is roughly consistent with the phenomenological hydrodynamic predictions if we assume that the density modes of the observed configurations keep memory of the liquid-state at the largest length-scales probed by the structure factor. We have also discussed the effect of non-asymptotic behaviors, such as the ones arising from the dispersivity of elastic constants and non-equilibrium effects, and finite-size effects.
The kind of study we have done may offer a different viewpoint to analyze vortex phases, their transitions, and also to characterize the somewhat unavoidable nonequilibrium relaxation effects in disordered systems. This includes not only static vortex phases as the ones we have analyzed in this work, but also dynamical vortex phases such as current-driven lattices [4,15,[34][35][36][37], where genuine non-equilibrium stationary effects may add interesting extra ingredients to the density fluctuations. In all cases, it would be interesting to analyze whether new hyperuniform patterns emerge, either at the level of the whole system or in a subsystem (as the constant-z cross section we study here), and to check if they share some of the interesting response properties discussed in Ref. 2 for general hyperuniform states of matter.