# Hyperuniform vortex patterns at the surface of type-II superconductors

Gonzalo Rumi, Jazmín Aragón Sánchez, Federico Elías, Raúl Cortés

Maldonado, Joaquín Puig, Néstor René Cejas Bolecek, Gladys Nieva, Marcin

Konczykowski, Yanina Fasano, Alejandro B Kolton

► To cite this version:

Gonzalo Rumi, Jazmín Aragón Sánchez, Federico Elías, Raúl Cortés Maldonado, Joaquín Puig, et al. Hyperuniform vortex patterns at the surface of type-II superconductors. Physical Review Research, 2019, 1 (3), pp.033057. 10.1103/PhysRevResearch.1.033057. hal-03800530

HAL Id: hal-03800530

https://hal.science/hal-03800530

Submitted on 6 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Hyperuniform vortex patterns at the surface of type-II superconductors

Gonzalo Rumi,¹ Jazmín Aragón Sánchez,¹ Federico Elías,¹ Raúl Cortés Maldonado,¹ Joaquín Puig,¹

Néstor René Cejas Bolecek,¹ Gladys Nieva,¹ Marcin Konczykowski,² Yanina Fasano,¹,\* and Alejandro B. Kolton¹

¹Centro Atómico Bariloche and Instituto Balseiro, CNEA, CONICET, and Universidad Nacional de Cuyo, Bariloche, Argentina

²Laboratoire des Solides Irradiés, CEA/DRF/IRAMIS, Ecole Polytechnique, CNRS, Institut Polytechnique de Paris, Palaiseau, France

<figure><img src="image_2.png" /></figure>

(Received 17 July 2019; revised manuscript received 6 September 2019; published 29 October 2019)

A many-particle system must possess 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 ~$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 nonhyperuniform 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 $Bi_2Sr_2CaCu_2O_{8+\delta}$ samples in the mixed state. For most cases, we find hyperuniform two-dimensional point patterns at the superconductor surface with an effective exponent $\alpha_{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. We argue these predictions may be observable and propose further imaging experiments to test them independently.

DOI: 10.1103/PhysRevResearch.1.033057

I. INTRODUCTION

Hyperuniform point patterns, defined by a complete sup-
pression of density fluctuations in the large-wavelength limit
[1], have attracted great interest in recent years. Such behavior
can spontaneously emerge, following either equilibrium or
nonequilibrium protocols, in disordered ground states, glass
formation, jamming, Coulomb systems, spin systems, pho-
tonic 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–3]. Hy-
peruniform 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 a 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
among elasticity, quenched disorder, thermal fluctuations,
driving forces, anisotropic effects, and finite-size effects, ei-
ther 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 nonequilibrium 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
nonhyperuniform distribution of point pins. The work of
Ref. [5] also predicts that vortex matter may exhibit disor-
dered 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 the presence of a Poisson distri-
bution 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 hype-
run uniformity [2] in a system with a nondivergent compression
modulus for nearly uniform deformations. On the other hand,

<footnote>*Corresponding author: yanina.fasano@gmail.com</footnote>
hyperuniformity is expected for clean thin-film superconductors with the penetration length much larger than the sample thickness and logarithmic vortex-vortex repulsion. Reference [5] shows that this expectation is satisfied also in the presence of Poisson-distributed point pins.

Interestingly, as we show in this work, three-dimensional systems made of many nonrigid interacting lines directed along the z direction are generically expected to follow a hyperuniform point pattern for any constant-z cross section of the embedding space. This behavior arises mainly from the constraint that elastic lines cannot start or 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₂Sr₂CaCu₂O₈+δ 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 finite-size 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 Sec. II A, we define the main observables that we will use to study hyperuniformity in magnetically decorated vortex structures. In Sec. 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 Sec. III, we study the large-wavelength density fluctuations in magnetically decorated vortex structures nucleated in Bi₂Sr₂CaCu₂O₈+δ samples with different types of disorder.

<figure><img src="image_3.png" /><figcaption>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 z direction. 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 two-dimensional point pattern at each constant-z cross section (full circles), in spite of repulsion between lines being short-ranged. The pattern frozen at the top layer can be accessed experimentally via bitter magnetic decoration experiments (see Sec. III).</figcaption></figure>

In Sec. IV, we discuss the theoretical interpretation of these results and report predictions expected for vortex matter in ad hoc experiments that might shed more light on 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

We will be interested in the two-dimensional large-wavelength density fluctuations at a constant-z cross section of a three-dimensional array of vortex lines nucleated in a real sample of thickness L. The average magnetic induction at equilibrium, $B = B\hat{z}$, determines the number of vortices per unit area in any slice, $n_0 = B/\Phi_0$, with $\Phi_0 = 2.07 \times 10^{-7} \text{ 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 approximation, that vortex lines do not present overhangs or 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 \gg 1$ points at a constant-z cross section is defined as

$$S(\mathbf{q}, z) = N^{-1} \left| \sum_{n=1}^{N} e^{-iq \cdot r_n(z)} \right|^2 - N \delta_{\mathbf{q}, 0} \quad (1)$$

with $\mathbf{q} = (q_x, q_y)$. For thick systems in the z direction, $L \to \infty$, and the bulk value of $S(\mathbf{q}, z)$ is expected to be independent of z, namely $S(\mathbf{q}, z) = S(\mathbf{q})$ [9].

A constant-z cross-section pattern is considered hyperuniform if its two-dimensional large-wavelength (small-wavevector) density fluctuations are suppressed:

$$\lim_{q \to 0} S(\mathbf{q}) = 0. \quad (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,

$$S(\mathbf{q}) \sim q^{\alpha}, \quad (q \to 0), \qquad (3)$$

with $\alpha \ge 0$ being a characteristic exponent. In the asymptotic limit $N \gg 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$, $\sigma_N^2(R) = \langle N^2(R) \rangle - \langle N(R) \rangle^2$. In the latter, $\langle \cdots \rangle$ denotes average over randomly distributed windows and the expression is valid in the large-$R$ limit. The variance is predicted to scale as [1]

$$\sigma_N^2(R) \sim \begin{cases} R^{d-1}, & \alpha > 1, \\ R^{d-1} \ln R, & \alpha = 1 \\ R^{d-\alpha}, & 0 < \alpha < 1. \end{cases} \quad (R \to \infty), \quad (4)$$

For instance, a Poisson distribution of points is not hyperuniform since the number variance is extensive with the window volume, i.e., $\sigma_N^2 \sim R^d$. A hyperuniform system is considered to be class I if $\alpha > 1$, class II if $\alpha = 1$, and class III if $0 < \alpha < 1$. A large and diverse list of physical and mathematical systems for each class is given in Ref. [2].

## 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, $\lambda_{ab} \ll L$ (see Fig. 1), vortex lines repel each other with a roughly exponential dependence on distance for scales of order $\lambda_{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 \to 0$, $q_z \to 0$ limit of quasiuniform compressions and tilting deformations.

The starting point for the elastic line hydrodynamics of line liquids is the Landau free-energy functional,

$$F = \frac{1}{2n_0^2} \int d^2\mathbf{r}d\mathbf{z} \int d^2\mathbf{r}' \int d\mathbf{z}' \\ \times [c_{44}(\mathbf{r}-\mathbf{r}', \mathbf{z}-\mathbf{z}')\mathbf{t}(\mathbf{r}, \mathbf{z}) \cdot \mathbf{t}(\mathbf{r}', \mathbf{z}') \\ + c_{11}(\mathbf{r}-\mathbf{r}', \mathbf{z}-\mathbf{z}')\delta n(\mathbf{r}, \mathbf{z})\delta n(\mathbf{r}', \mathbf{z}')] \\ + \int d^2\mathbf{r} \int d\mathbf{z} V_D(\mathbf{r}, \mathbf{z})\delta n(\mathbf{r}, \mathbf{z}), \qquad (5)$$

where $\mathbf{r} = (x, y)$, $c_{44}(\mathbf{r}-\mathbf{r}', \mathbf{z}-\mathbf{z}')$, and $c_{11}(\mathbf{r}-\mathbf{r}', \mathbf{z}-\mathbf{z}')$ are the nonlocal tilt and compression modulii [12] with Fourier transforms $c_{11}(q, q_z)$ and $c_{44}(q, q_z)$,

$$\delta n(\mathbf{r}, \mathbf{z}) = \sum_{j=1}^{N} \delta[\mathbf{r} - \mathbf{r}_j(\mathbf{z})] - n_0 \qquad (6)$$

are the two-dimensional vortex density fluctuations at layer $z$ around its mean value $n_0$, and

$$\mathbf{t}(\mathbf{r}, \mathbf{z}) = \sum_{j=1}^{N} \frac{d\mathbf{r}_j}{dz} \delta[\mathbf{r} - \mathbf{r}_j(\mathbf{z})] \qquad (7)$$

is the two-dimensional tangent field density for a collection of $N \gg 1$ vortex-lines positioned at $\mathbf{r}_j(\mathbf{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 \equiv B/\Phi_0$ (see Fig. 1). The last term in $F$ describes the coupling of the vortex density with the pinning potential $V_D(\mathbf{r}, \mathbf{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 $\delta n(\mathbf{q}, \mathbf{z}) = n(\mathbf{q}, \mathbf{z}) - n_0 A \delta_{q,0}$, where $\mathbf{q}$ and $\mathbf{z}$ are the wave vectors in the in-plane directions and in the average-$z$ direction of the lines. Fluctuations are conveniently quantified by the full three-dimensional structure factor

$$n_0 S^{3d}(\mathbf{q}, \mathbf{q}_z) = \overline{\langle |\delta n(\mathbf{q}, \mathbf{q}_z)|^2 \rangle}, \qquad (8)$$

where $\cdots$ and $\langle \cdots \rangle$ denote averages over disorder and thermal fluctuations, respectively. These averages can be calculated, at equilibrium, from Eq. (5) with the constraint

$$\partial_z \delta n + \nabla_{\perp} \cdot \mathbf{t} = 0 \qquad (9)$$

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 as $\tilde{c}_{44}$ and $\tilde{c}_{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) = \langle |\delta n(q, z)|^2 \rangle \quad (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$ (see Appendix A)

$$S(q) \propto \int_{0}^{2\pi/s} dq_z S^{3d}(q, q_z), \quad (11)$$

where $s$ sets an ultraviolet cutoff coming from the discretization in the z direction or, more physically, from the superconducting layer spacing [10].

## 1. 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_{\text{liq}}^{3d}(q, q_z)$, and the solid Abrikosov lattice, $S_{\text{sol}}^{3d}(q, q_z)$, vortex phases at small enough wave vectors,

$$S_{\text{liq}}^{3d}(q, q_z) \approx S_{\text{sol}}^{3d}(q, q_z) = \frac{n_0 k_B T q^2}{q^2 c_{11}(q, q_z) + q_z^2 c_{44}(q, q_z)}, \quad (12)$$

In this expression obtained from Eq. (5) [7,10,14,16], $k_B$ is the Boltzmann constant and $T$ is 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 ($qa_0 < 1$). Equation (12) can be obtained by relating the spontaneous density fluctuations measured by $S_{\text{liq}}^{3d}(q, q_z)$ with the associated linear response function $n_0 q^2 [q^2 c_{11}(q, q_z) + q_z^2 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\pi/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 by following Ref. [10]

$$S_{\text{liq}}(q) \approx S_{\text{sol}}(q) = \frac{n_0 k_B T}{\sqrt{c_{44}(q, 0)c_{11}(q, 0)}} q \quad (13)$$

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 \to 0$, we find that

$$S_{\text{liq}}(q) \approx S_{\text{sol}}(q) \propto q. \quad (14)$$

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 form ($\alpha = 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].

## 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 $\overline{V_D(r, z)V_D(r', z')} = \Delta\delta(r' - r)\delta(z - z')$, produces an additive correction to the vortex liquid three-dimensional structure factor [Eq. (12)] such that

$$S_{\text{pin-liq}}^{3d}(q, q_z) = n_0 \Delta \left[ \frac{n_0 q^2}{q^2 c_{11}(q, 0) + q_z^2 c_{44}(q, 0)} \right]^2 + S_{\text{liq}}^{3d}(q, q_z) \quad (15)$$

with $\Delta$ 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$,

$$S_{\text{pin-liq}}(q) = \frac{n_0^3 \Delta}{\sqrt{c_{44}(q, 0)c_{11}^3(q, 0)}} q + S_{\text{liq}}(q). \quad (16)$$

If we again assume that $c_{11}(q, 0)$ and $c_{44}(q, 0)$ tend to constants as $q \to 0$, we get

$$S_{\text{pin-liq}}(q) \propto q, \quad (17)$$

as also found for the clean case. Note, however, that the prefactor is different: Only when $n_0 \Delta/c_{11}(q, 0) < k_B T$ will the system cross over to the line liquid in clean samples.

This case of $\alpha = 1$ or class-II hyperuniformity is different to the $\alpha = 2$ class-I hyperuniformity that is found, for instance, for pancake vortices interacting logarithmically in a thin-film superconductor or the two-dimensional one-component 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 $\sim 1/r$ repulsion. Such case is indeed known to
have a two-dimensional compression modulus diverging as ~1/q, which complies with the fluctuation-compressibility theorem. Hyperuniformity in the liquid phase is rather robust since it 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 $\alpha = 1$ class-II hyperuniform system, either in the absence or presence of weak uncorrelated disorder.

### 3. 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_\alpha(\mathbf{r})$ measuring the distance of a vortex with respect to the corresponding perfect lattice position, such that $u_\alpha(R_i, z) = u_{i,\alpha}(z)$ is the displacement of the ith 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

$$n(\mathbf{r}, z) = n_0(1 - \partial_{\mathbf{r}} \mathbf{u}(\mathbf{r}, z) + \dots), \quad (18)$$

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

$$S_{\text{sol}}^{3d}(q, q_z) \approx n_0 q^2 \overline{\langle |u_L(\mathbf{q}, q_z)|^2 \rangle}, \quad (19)$$

where $u_L(\mathbf{q}, q_z)$ is the component of the displacement field in the direction parallel to $\mathbf{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

$$\overline{\langle |u_L(\mathbf{q}, q_z)|^2 \rangle} \sim \frac{k_B T}{q^2 c_{11}(q, q_z) + q_z^2 c_{44}(q, q_z)}. \quad (20)$$

The structure factor of the lattice is thus essentially the same as for the liquid,

$$S_{\text{sol}}^{3d}(q, q_z) \approx \frac{n_0 k_B T q^2}{q^2 c_{11}(q, q_z) + q_z^2 c_{44}(q, q_z)} \approx S_{\text{liq}}^{3d}(q, q_z), \quad (21)$$

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 $\alpha = 1$ hyperuniformity as the liquid,

$$S_{\text{sol}}(q) \propto q, \quad (22)$$

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 temperature-dependent 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

$$\overline{\langle |u_L(q, q_z)|^2 \rangle} \sim [q^2 c_{11}(q, q_z) + q_z^2 c_{44}(q, q_z)]^{-\frac{d+2\zeta}{2}}, \quad (23)$$

with $d = 3$ being the space dimension and $\zeta$ being 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 the presence of thermal noise, we have $\zeta = 1 - d/2$, thus reducing to Eq. (12). When $\zeta > 0$, this exponent is the so-called roughness exponent of the lattice, such that the mean square displacement grows as

$$\langle u_L^2 \rangle \sim \int dq \, dq_z \, q^{d-2} \overline{\langle |u_L(\mathbf{q}, q_z)|^2 \rangle} \sim R^{2\zeta} \quad (24)$$

with $R$ being the observation scale or system linear size in the x-y plane. For $\zeta = 0$, the lattice is logarithmically rough $\langle u_L^2 \rangle \sim \log R$, while for $\zeta < 0$ it 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, $S_{\text{pin-sol}}(q) \sim q^{3-d-2\zeta}[c_{44}c_{11}^{d+2\zeta-1}]^{-1/2}$, and thus for $d = 3$

$$S_{\text{pin-sol}}(q) \sim \frac{q^{-2\zeta}}{\sqrt{c_{44}(q, 0)c_{11}^{2+2\zeta}(q, 0)}}. \quad (25)$$

Therefore, under these assumptions, the elastic system is hyperuniform only if $\zeta < 0$ (displacement field macroscopically flat) while for $\zeta = 0$ is marginally hyperuniform, i.e., $\alpha = 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 $\zeta \to 0$ asymptotically. This so-called random-periodic regime (RP) of the Bragg glass phase exists for $qR_a \ll 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 random-manifold regime (RM), where $\zeta \equiv \zeta_{RM} \approx 0.2$. For scales shorter than $R_c$, the system crossovers to the Larkin regime, where $\zeta \equiv \zeta_L = (4-d)/2 = 1/2$. In none of these regimes will the Bragg glass phase present suppressed density fluctuations.
In summary, 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].

### 4. CD correlated disorder: Liquid and Bose glass

The case of the correlated disorder generated by randomly distributed CD, such that $V_D(\mathbf{r}, z)V_D(\mathbf{r}', z') = \Delta_1\delta(\mathbf{r}' - \mathbf{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 of $\tilde{c}_{44}$ approaching the transition temperature, while $\tilde{c}_{11}$ remains finite. There is an exception for the putative Mott glass phase at the matching field $B = B_\Phi \equiv n_{\text{col}}\Phi_0$, where $\tilde{c}_{11}$ is also expected to diverge. The predicted three-dimensional structure factor for the vortex liquid in the presence of CD disorder is

$$S_{\text{col-liq}}^{3d}(\mathbf{q}, q_z) \approx \frac{n_0 k_B T q^2}{q^2 \tilde{c}_{11}(q, q_z) + q_z^2 \tilde{c}_{44}(q, q_z)} + \Delta_1 \frac{n_0^3}{\tilde{c}_{11}^2(q, 0)} \delta(q_z), \quad (26)$$

where $\tilde{c}_{11}$ and $\tilde{c}_{44}$ are the compression and tilt modulii in the presence of correlated CD disorder and $\Delta_1$ measures the strength of the columnar disorder. Integrating over $q_z$, we obtain the two-dimensional structure factor

$$S_{\text{col-liq}}(q) \approx \frac{n_0 k_B T}{\sqrt{\tilde{c}_{44}(q, 0)\tilde{c}_{11}(q, 0)}} q + \Delta_1 \frac{n_0^3}{\tilde{c}_{11}^2(q, 0)}. \quad (27)$$

This expression yields a crossover wave vector

$$q_{\text{CD}} = \frac{n_0^2 \Delta_1}{\tilde{c}_{11} k_B T} \sqrt{\frac{\tilde{c}_{44}}{\tilde{c}_{11}}} \quad (28)$$

such that for scales $q < q_{\text{CD}}$, the density fluctuations in the liquid are dominated by the pinning introduced by CD [second term in Eq. (27)], while for $q > q_{\text{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, $\tilde{c}_{44}$ increases rapidly and diverges at the transition, concomitantly with $q_{\text{CD}} \to \infty$. We hence expect

$$S_{\text{col-sol}}(\mathbf{q}) \approx \Delta_1 \frac{n_0^3}{\tilde{c}_{11}^2(q, 0)} \quad (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 nonhyperuniform 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_{\text{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 $\Delta_1$.

### 5. 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_{\text{FS}}(L)$. Indeed, by a simple dimensional analysis of the three-dimensional structure factors in these hyperuniform cases [Eqs. (12) and (15)], we get

$$q_{\text{FS}} \sim \sqrt{\frac{c_{44}(q_{\text{FS}}, 2\pi/L) 2\pi}{c_{11}(q_{\text{FS}}, 2\pi/L) L}} \quad (30)$$

Physically, $q_{\text{FS}}$ 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,

$$S^{3d}(\mathbf{q}, z_1 - z_2) = \frac{n_0 k_B T}{c_{11}(q, 0) \xi_{\parallel}(q)} e^{-|z_1 - z_2| / \xi_{\parallel}(q)} \quad (31)$$

with

$$\xi_{\parallel}(q) = q^{-1} \sqrt{c_{44}(q, 2\pi/L)/c_{11}(q, 2\pi/L)}. \quad (32)$$

Therefore, Eq. (30) is equivalent to the physical condition $\xi_{\parallel}(q_{\text{FS}}) = L$. For $q < q_{\text{FS}}$, 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 $\xi_{\parallel} = L$,

$$S(\mathbf{q}) = \frac{n_0 k_B T}{c_{11}(q, 0)L}. \quad (33)$$

Provided $L > \lambda_{ab}$ at the crossover, rigid vortices repel with a short-range interaction, yielding [14,15]

$$c_{11}(q) = 2\epsilon_0 \frac{2\pi \lambda_{ab}^2}{1 + \lambda_{ab}^2 q^2} \frac{B^2}{\Phi_0^2} \quad (34)$$

with $\epsilon_0 = (\Phi_0/4\pi\lambda_{ab})^2$ being the line tension. Therefore, since the length scale in the x-y plane probed is $q^{-1} \gg \lambda_{ab}$, we have $S(\mathbf{q}) \to \text{const}$ as $q \to 0$. In other words, finite-size effects kill the class-II hyperuniformity of the three-dimensional 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_{\text{FS}}(L)$ comes from first, a dimensional analysis of the competition between tilting and compression elastic responses, and second, from the nondispersive behavior of the three-dimensional elastic constants.

### 6. Two-dimensional systems

It is interesting to discuss Eq. (33) in the limit of very thin superconductors such that $L \ll \lambda_{ab}$. In this limit, the effective penetration length becomes $\Lambda = 2\lambda_{ab}^2/L$. If the observation length scale in the x-y plane is $q\Lambda > 1$, two different power-law behaviors can be expected. At very small length scales,
<table><caption>TABLE I. Analytical predictions for the type of hyperuniformity for different three-dimensional vortex-line phases at constant-z cross sections in the thermodynamic limit.</caption><thead><tr><th>Phases</th><th>Hyperuniformity</th></tr></thead><tbody><tr><td>Liquid without disorder</td><td><b>Yes</b>: class II α = 1</td></tr><tr><td>Liquid with weak uncorrelated disorder</td><td><b>Yes</b>: class II α = 1</td></tr><tr><td>Liquid with correlated disorder by CDs</td><td><b>No</b></td></tr><tr><td>Abrikosov crystal</td><td><b>Yes</b>: class II α = 1</td></tr><tr><td>Bragg glass</td><td><b>Marginal</b>: α = 0</td></tr><tr><td>Bose glass</td><td><b>No</b></td></tr><tr><td>Mott glass</td><td><b>Yes</b>: unknown class</td></tr></tbody></table>

$q\Lambda \gg 1$, interactions are logarithmic and

$$c_{11}(\mathbf{q}) = 2\epsilon_0 \frac{2\pi}{q^2} \frac{B^2}{\Phi_0^2} \qquad (35)$$

in Eq. (33), so the system at equilibrium is class-I hyperuniform with $\alpha = 2$,

$$S(q) \rightarrow q^2 \qquad (36)$$

for $\log(1/x)$ interactions. Recent simulations with this type of interaction show that the value $\alpha = 2$ is robust under the presence of quenched uncorrelated disorder [5].

For larger length scales, $q\Lambda > 1$, the interaction becomes Coulomb-like and

$$c_{11}(\mathbf{q}) = 2\epsilon_0\Lambda \frac{2\pi}{q} \left(\frac{B}{\Phi_0}\right)^2, \qquad (37)$$

yielding, at equilibrium, a class-II hyperuniform vortex pattern with $\alpha = 1$ and thus

$$S(\mathbf{q}) \rightarrow q \qquad (38)$$

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, the 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 $\alpha = 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, as a result of their similar uniaxial compression properties at small wave vectors.

## 7. 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 toward 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 way $\tilde{c}_{11}$ diverges with $q$ has not been reported for the Mott glass, and hence we cannot guess its hyperuniformity class yet. All these predictions assume three-dimensional vortex line phases and nondispersive 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 on the approximated elastic modulii of the continuum description.

We model vortices as elastic lines discretized in the z direction, such that $\mathbf{r}_i(z) = (x_i(z), y_i(z))$ describe their two-dimensional coordinates at the layer z, with $z = 0, \dots, 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[\{\mathbf{r}_i(z)\}] = E_1 + E_{\text{int}}$ such that each line has an elastic tension energy given by Hooke coupling with strength $k$,

$$E_1 = \sum_{i=1}^{N} \sum_{z=0}^{N_z-1} \frac{k}{2} |\mathbf{r}_i(z+1) - \mathbf{r}_i(z)|^2 \qquad (39)$$

and the repulsive interaction energy between vortex lines is modeled as

$$E_{\text{int}} = \sum_{i \neq j} \sum_{z=0}^{N_z-1} \epsilon_0 K_0 \left( \frac{|\mathbf{r}_j(z) - \mathbf{r}_i(z)|}{\lambda_{ab}} \right) \qquad (40)$$

with $K_0(x)$ being the zero-order modified Bessel function of the second kind. Equation (40) is derived from the London model, while Eq. (39) with $k \propto \epsilon_0$ is a harmonic local approximation for the single-vortex elastic tension. In layered superconductors such as the one experimentally studied here, the elastic tension energy arises from the attractive electromagnetic and Josephson couplings between pancakes of different superconducting layers. We consider an overdamped
<figure><img src="image_4.png" /><figcaption>FIG. 2. Typical structure factors for a liquid and a solid at equilibrium at finite T for an ideal clean superconductor, as obtained from molecular dynamics simulations of an elastic-line array with short-ranged elasticity and repulsive short-ranged interactions. Axial symmetry in both phases is observed in the large-wavelength (low-q) limit.</figcaption></figure>

Langevin dynamics at a temperature T

$$ \eta \partial_t \mathbf{r}_i(z) = -\frac{\delta E}{\delta \mathbf{r}_i(z)} + \xi(\mathbf{r}_i(z), t), \quad (41) $$

$$ \langle \xi(\mathbf{r}_i, t) \xi(\mathbf{r}'_j, t') \rangle = 2\eta k_B T \delta_{ij} \delta(t-t'), \quad (42) $$

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₀ 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 both have 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 axially symmetric elastic constants. Thus, for the q-range we work with the angular average of the structure factor, S(q) = (2π)⁻¹ ∫₀²π dΨS(q cos(Ψ), q sin(Ψ)), with Ψ being the azimuthal angle, 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. The figure shows results for different lattice spacings a₀/λab = 1, 3, 4, where A = Na₀²√3/2 is the area of the computational box, and different number of layers N_z = 16, 36, 50 in the z direction.

<figure><img src="image_5.png" /><figcaption>FIG. 3. Normalized two-dimensional structure factor obtained from molecular dynamics simulations of a three-dimensional elastic-line array in the ideal clean case, for various lattice spacings a₀ and number of layers N_z. We use a finite-size crossover scale q FS ∝ a₀/N_z and a normalization S(q FS ) ∝ T a 2 $_{0}$/N_z. (a) Liquid phase at equilibrium and (b) solid Abrikosov phase. Both phases display the same finite-size crossover at q FS . For q < q FS , the system effectively becomes a nonhyperuniform two-dimensional system, while for q > q FS it shows a fair S ( q ) ∼ q hyperuniform behavior 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 (dash-dotted line). The insert shows non-normalized data.</figcaption></figure>

There are clear finite-size effects, i.e., N_z dependence, in the z direction, and they are amplified when a₀ 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 FS such that S(q)/S(q FS ) ∼ G(q/q FS ), with G(x) ∼ 1 for x ≪ 1 and a unique G(x) for x ≫ 1. For the range of a₀ analyzed, we find a good collapse in the low-q region using the crossover wave vector q FS = a₀/N_z and S(q FS ) ∝ T a 2 $_{0}$/L, both for the liquid and solid phases. Therefore, the saturation of S(q) for q < q FS is actually a finite-size effect, and it should disappear in the thermodynamic N_z → ∞ limit. The dependence of q FS with N_z ∝ L confirms the size effects predicted in Sec. II B 5. In addition, the shape of the master curve for q > q FS , G(x) ∼ x, confirms an α ≈ 1 typical of class-II hyperuniformity [see Eq. (13)].
The evolution of $q_{FS}$ with $a_0$ is less universal than the dependence with $N_z$ since it relies on the precise functionality of the elastic modulii $c_{44}$ and $c_{11}$ that depend on the type of intervortex 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) \propto B^2/(1 + q^2\lambda_{ab}^2)$, the single-vortex tilt modulus $c_{44} \propto \epsilon_0 \sim B^2a_0^2$ (a constant in this particular model), and approximate Eq. (30) as $q_{FS} \approx \sqrt{c_{44}/c_{11}(0,0)}/L$, we get $q_{FS} \propto a_0/L$. Also, from Eq. (33) we get $S(q_{FS}) \sim n_0k_B T/c_{11}L \propto Ta_0^2/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_{FS}$.

For $q > q_{FS}$, the master curve presents deviations from $S(q) \sim q$ for $q \gg q_{FS}$, 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) \propto q/\sqrt{c_{44}c_{11}(q,0)} \sim q\sqrt{1+\lambda_{ab}^2q^2} \approx q(1+q^2\lambda_{ab}^2/2+\cdots)$. Therefore, this microscopic model confirms the existence of finite-size effects 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 $\alpha = 1$ behavior observed at $q > q_{FS}$ controlled by the dispersive behavior of $c_{11}$, emerging when $q\lambda_{ab} \sim 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_2Sr_2CaCu_2O_{8+\delta}$. 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 two-dimensional 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 $\sim 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 the experimental parameters in our measurements were set to maximize the scattered signal, and we found that in this case the experimental resolution at low $q$ is not good enough as to ascertain whether hyperuniformity is present in the vortex system. Future small-angle neutron experiments optimizing the experimental configuration as to improve the resolution for the low $q$ range will certainly shed more light on this issue.

The choice of the extremely layered $Bi_2Sr_2CaCu_2O_{8+\delta}$ system is based on 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 single-crystals 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 closed-cycle 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 correspond to a sample irradiated with an electron density of $1.7 \times 10^{19} \text{ e/cm}^2$ [23]. The irradiation of samples with energetic ($\approx 1$ GeV) heavy ions of Xe and Pb resulted in samples with CD correlated disorder characterized by the matching field $B_\Phi$. The studied samples have a low density of CD: $B_\Phi = 30$ and $60$ G for the Xe-irradiated samples and $B_\Phi = 45$ and $100$ G for the Pb-irradiated ones. The whole set of studied samples are single crystals with typical dimensions of $1 \times 1 \text{ mm}^2$ in-plane area and thickness ranging from 300 to 600 $\mu$m. The size-effect mentioned in previous sections comes from the finite thickness of the sample. Details on the vortex phase diagrams for these samples can be found in Ref. [24].

We image individual vortex positions in a typical field of 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 [25], observed as black dots in the inverted scanning-electron-microscopy images of Fig. 4. Further details in the field-cooling decoration protocol followed in this case can be found in Ref. [25]. For every studied sample, several ($\approx 10$) magnetic decoration experiments were performed in freshly cleaved surfaces, eventually at different applied fields. We studied $\approx 30$ pristine, 1 electron-irradiated, and 10 heavy-ion irradiated single crystals.

Figure 4 shows examples of magnetic decoration images at various applied fields in the three types of studied samples. Figures 4(a) and 4(b) correspond to snapshots of vortex positions taken at a field well within the quasicrystalline Bragg glass phase for pristine and electron-irradiated samples. For these two types of samples, at fields larger than 15 G, the vortex structure is single crystalline, presents quasi-long-range positional order, and has 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 $\approx 2\%$ of vortices are involved in defects, mainly edge dislocations. For fields smaller than 15 G, the structure breaks into small crystallites for pristine [25] 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 [25]. 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 [24].
<figure><img src="image_6.png" /><figcaption>FIG. 4. Magnetic decoration images of the vortex structure (black dots) nucleated in Bi$_2$Sr$_2$CaCu$_2$O$_{8+δ}$ samples with different types of disorder: (a) pristine sample with point pins; (b) electron-irradiated sample with extra point pins; (c) heavy-ion-irradiated sample with a low density of columnar defects ($B_Φ = 30$ G, Xe-irradiated). The magnetic induction *B* controlling vortex density is indicated in every panel. In all cases, the white bar corresponds to 5 µm. Each decoration image has overimposed Delaunay triangulations joining near-neighbor vortices with blue lines. Non-sixfold-coordinated vortices are highlighted in red.</figcaption></figure>

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 G, see Fig. 4,

and for fields smaller than ≈15 G 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 short-range 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 Kelvins below the temperature at which magnetic response becomes irreversible [29]. This length scale is shorter than the relevant for the hyperuniformity analysis. Because of 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\pi/a_0$. We will further discuss this issue later.

The images of Fig. 4 are just some examples, but typically we study panoramic images of the vortex structure with no less than 1500 and up to 15 000 vortices. We digitalize vortex positions in these large fields of view and then we calculate the structure factor $S(q)$ and the number variance $\sigma_N^2$. As the field of view becomes larger, the value of $q = \sqrt{q_x^2 + q_y^2}$ decreases, so that we can access to calculate the two-dimensional structure factor at the top surface. Figure 5 shows the corresponding two-dimensional structure factors $S(q_x, q_y)$ for the examples of vortex structures shown in Fig. 4. The $S(q_x, q_y)$ values are angularly averaged as to obtain the $S(q)$ data shown in Fig. 6. We find that $S(q)$ is independent of the image resolutions considered in our experimental data and grid spacings considered for the calculation of $S(q)$.

We would like to point out some important technical difficulties that are quite specific to the study of the low-$q$ 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 windowing 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 the azimuthal angle $\Psi$ values outside the “+”-shaped region, instead of averaging over all $\Psi$. 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\pi/W_{fov}$, with $W_{fov}$ being a characteristic linear size of the field of view. On the other hand, as we do not have an ensemble of many magnetically decorated vortex configurations at exactly the same position of the samples to average over, inevitably, as the $q$ is lower, the statistical fluctuations in $S(q)$ are bigger.
<figure><img src="image_7.png" /><figcaption>FIG. 5. Two-dimensional structure factors $S(q_x, q_y)$ for the examples of vortex structures shown in 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)$.</figcaption></figure>

Henceforth, in practice, although our images span distances as large as $W_{\text{fov}} \sim 50a_0$, we end up analyzing density fluctuations with $q > q_{\text{min}} \gtrsim q_0/6 = 1/6(2\pi/a_0)$.

Typical results for $S(q)$ are shown in Fig. 6 (top panels) for the pristine [Fig. 6(a)], electron-irradiated [Fig. 6(b)], and heavy-ion-irradiated [Fig. 6(c)] samples. The dashed and dash-dotted lines are guides to the eye, showing $q^2$ and $q$ evolutions, respectively. For the lowest $q/q_0$, the behavior is consistent with hyperuniformity with $\alpha = 1$ for the three samples. Indeed, red lines in the figure shows a fit to the data yielding $\alpha_{\text{eff}} \sim 1$. Another way to analyze the density fluctuations is to study the distance evolution of the vortex-number variance $\sigma_N^2$ [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 $\sigma_N^2(R)$ is then computed. Figure 6 (bottom panels) shows the results of the $R$ evolution of $\sigma_N^2$ 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 $\sigma_N^2 \sim R/a_0$ fairly describes the data (see the dash-dotted line). Indeed, the red lines in Figs. 6(d)–6(f) correspond to functions $\approx R^{2-\alpha_{\text{eff}}}$ [see Eq. (4)] with the $\alpha_{\text{eff}}$ obtained from the fits of the $S(q) \sim q^{\alpha_{\text{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) \sim q^{\alpha_{\text{eff}}}$ for all our studied samples and vortex densities (magnetic field). The effective power-law exponent $\alpha_{\text{eff}}$ is obtained by fitting in the range $q/q_0 \in [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 (see Appendix C for the experimental error bar estimate). Although $S(q)$ and $\sigma_N^2$ qualitatively agree, we have found that a systematic fit of $\alpha_{\text{eff}}$ using the expected scaling of $\sigma_N^2$ is empirically more difficult due to the strong nonasymptotic corrections to the number variance [2]. As we can observe, $\alpha_{\text{eff}} \approx 1 \pm 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) \sim 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) \sim q^{\alpha_{\text{eff}}}$ with $\alpha_{\text{eff}} \approx 1 \pm 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.
<figure><img src="image_8.png" /><figcaption>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 (CD correlated disorder) Bi$_2$Sr$_2$CaCu$_2$O$_{8+δ}$ samples. The vortex density in every case is indicated. Data are shown as a function of $q/q_0$ with $q_0 = 2\pi/a_0$ the Bragg wave vector. The red line is the best power-law fit $S(q) \sim q_\text{eff}^\alpha$ for the low wave-vector range $a_0/6 < q/q_0 < a_0/2$. Number variance $\sigma_N^2$ as a function of $R/a_0$ corresponding to the same samples of the top panel: (d) pristine, (e) electron-irradiated, and (f) Xe-irradiated samples. The red line shows functions $\sim R^{2-\alpha_\text{eff}}$ with the $\alpha_\text{eff}$ obtained from fitting the $S(q)$ data of the top panels. The dashed line with slope $\alpha = 2$ and the dash-dotted line with slope $\alpha = 1$ are guides to the eye.</figcaption></figure>

However, this exponent contrasts with the predictions for the Bragg glass phase ($\alpha = 0$) and for the liquid and Bose glass phases with CD [$S(q) \to \text{const}$ as $q \to 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 nonequilibrium 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 small length-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_{\text{freez}} = T_{\text{freez}}(q)$. In a field-cooling experiment, at $T \sim T_{\text{freez}}(q)$ the density fluctuations of mode $q$ have a relaxation rate of the order of the experimental cooling rate.

<figure><img src="image_9.png" /><figcaption>FIG. 7. Effective exponents obtained fitting the structure factor as $S(q) \sim q^{\alpha_{\text{eff}}}$ for the magnetically decorated vortex structures in (a) pristine, (b) electron-irradiated, and (c) heavy-ion irradiated samples as a function of vortex densities $B$.</figcaption></figure>
Namely, at $T < T_{freez}(q)$ all modes with wavelength $\le 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 \sim a_0^{-1}) \le T_m$, with $T_m$ being 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 \ll a_0^{-1}) > T_m$. These slow modes with $q \to 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 nonequilibrated length scales. Glassy vortex dynamics thus prevents the system from reaching 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 capable of only 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 nonhyperuniform at equilibrium in an hydrodynamic approximation. How is it possible to obtain $\alpha_{eff} \approx 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 nonhyperuniform liquid on decreasing $q$. Observing the equilibrium structure at $q \to 0$ implies measuring at very large length scales $\sim 1/q_{CD}$ that might be hard to access experimentally in finite fields of view. To evaluate this hypothesis, we need to compute

$$q_{CD} \approx \frac{n_0^2(U_0^2 b_0^4/d^2)}{\bar{c}_{11} k_B T_{freez}(q_{CD})} \sqrt{\frac{\bar{c}_{44}}{\bar{c}_{11}}}, \quad (43)$$

where we have used in Eq. (28) that $\Delta_1 \approx (U_0^2 b_0^4/d^2)[1 + \mathcal{O}(b_0^2/d^2)]$ is the disorder correlator strength for CD with radius $b_0$ and mean separation $d \gg 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} \propto 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 $\Delta_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.

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+\delta}$ since dispersivity can induce, as discussed, a crossover from class-II hyperuniformity to a different behavior on increasing $q$. Since we access scales $q \cdot \lambda_{ab} \ll 1$, in-plane dispersivity is not important. However, dispersion in the $z$-direction is important when $q \cdot \lambda_c = \Gamma q \cdot \lambda_{ab} \sim 1$. The anisotropy parameter of our samples is $\Gamma \approx 170$, and then the condition $q \cdot \lambda_c \sim 1$ can easily be reached in our experimental fields-of-view since $\lambda_c(T_{freez}) = \Gamma \lambda_{ab}(T_{freez}) \approx 70 \mu m$. For an anisotropic superconductor the compression modulus is [11]

$$c_{11}(q, q_z) = \frac{B^2}{4\pi} \frac{1 + \lambda_c^2(q^2 + q_z^2)}{(1 + \lambda_{ab}^2(q^2 + q_z^2))(1 + \lambda_c^2 q^2 + \lambda_{ab}^2 q_z^2)} \quad (44)$$

while the tilt modulus reads

$$c_{44}(q, q_z) = \frac{B^2}{4\pi} \left[ \frac{1}{1 + \lambda_c^2 q^2 + \lambda_{ab}^2 q_z^2} \right] + c'_{44}(q_z) \quad (45)$$

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 \lambda_{ab} \sim 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,

$$c'_{44} \approx \frac{B^2}{4\pi} \left( \frac{a_0}{\lambda_{ab}} \right)^2 \quad (46)$$

If we assume $L \to \infty$ and integrate Eq. (12) over $q_z$ using the expressions of Eqs. (44) and (45), we obtain the two-dimensional structure factor $S(q)$ including dispersion effects. The result for our Bi$_2$Sr$_2$CaCu$_2$O$_{8+\delta}$ 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/\lambda_{ab}(T_{freez}) = 2$ corresponding to a magnetic induction of $\sim 40$ G. We considered $\lambda_{ab}(T_{freez}) \sim 0.4 \mu m$ at the irreversibility temperature $T_{irr}(B \sim 40$ G). 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 \lambda_{ab}(T_{freez})$ range used for fitting $\alpha_{eff}$ in our experimental $S(q)$ data.

In the case of an infinite sample, $S(q)$ presents a crossover, from $S(q) \sim q$ to $S(q) \sim q^2$ at larger $q \lambda_{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 with 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 $\lambda_c(T_{freez})$; see the top panel of Fig. 8. Then, in order to ascertain whether a structure is hyperuniform, it is crucial to obtain $\alpha_{eff}$ by fitting $S(q)$ data in a $q$ range located below this dispersivity-induced crossover.
<figure><img src="image_10.png" /><figcaption>FIG. 8. Constant-z cross-section structure factor $S(q)$ for the three-dimensional line liquid without disorder. Top panel: Infinitely thick sample $L/\lambda_{ab} \gg 1$ for various $a_0/\lambda_{ab}$ corresponding to different applied magnetic fields. At high fields, corresponding to the smallest $a_0\lambda_{ab}$ values, a crossover from the type II hyperuniform scaling $S(q) \sim q$ to $S(q) \sim q^2$ is observed at $q\lambda_c \approx 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 $a_0/\lambda_{ab} = 2$ (black curve) roughly corresponds to the typical field of 40 G applied in magnetic decorations for $Bi_2Sr_2CaCu_2O_{8+\delta}$ with $\lambda_{ab}(T_{freez}) = 0.4$ µm. Bottom panel: Finite-thickness effect in the structure factor $S(q)$ for various $L/\lambda_{ab}$ ratios at a fixed field such that $a_0/\lambda_{ab} = 2$. This effect produces the saturation of $S(q)$ in a region of small $q$ whose extension depends on $L/\lambda_{ab}$. The case of $L/\lambda_{ab} = 1000$ (black points) corresponds to a $Bi_2Sr_2CaCu_2O_{8+\delta}$ sample with a thickness of 400 µm, a typical value in real samples. The gray shaded area corresponds to the $q\lambda_{ab}$ range accessed experimentally for our typical field of view in decorations (obtained considered again $\lambda_{ab}(T_{freez}) = 0.4$ µm). Finite-size effects in $S(q)$ are not dominant for the experimental window accessed in our experiments in $Bi_2Sr_2CaCu_2O_{8+\delta}$. However, dispersivity effects do play a role in our experimental window.</figcaption></figure>

In addition, this fit cannot 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 with neglecting the dispersivity and strong anisotropy of the elastic constants for vortex matter in $Bi_2Sr_2CaCu_2O_{8+\delta}$. The figure highlights a black curve corresponding to the prediction for our approximate experimental situation, namely vortex density $\approx 40$ G and $L \sim 400$ µm. In our experiments, the saturation associated to the finite-size effect is well below the $q$ range taken into account in the $\alpha_{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/\lambda_{ab}(T_{freez}) \approx 10$, namely for $L \approx 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_2Sr_2CaCu_2O_{8+\delta}$ micron-sized samples are available, but the field of view is of only a few hundred of vortices, making it 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) \sim q$ behavior for our typical experimental situation (see the black curve in Fig. 8) starts at the upper half of our experimental fitting window. Therefore, although the easily cleavable $Bi_2Sr_2CaCu_2O_{8+\delta}$ system makes the experimental realization of this magnetic decoration study possible, the hyperuniformity analysis is muddied by the dispersivity effects inherent to its extremely layered nature. Then the $\alpha_{eff}$ values obtained in our study might be slightly overestimated due to this effect.

In order to study up to what point our obtained $\alpha_{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 toward 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 nonequilibrium situation of field-cooling experiments, as well as performing nonequilibrium simulations under the same conditions, would be very insightful.

In summary, we have shown that a directed three-dimensional elastic line array with short-ranged interactions in all directions can generate, at equilibrium, a two-dimensional hyperuniform point pattern at their cross sections. Strikingly, this means that, at thermal equilibrium, every constant-z cross section becomes incompressible in the limit $q \to 0$ in response to a weak external field applied only in such cross section, while the response of the three-dimensional system to an external field applied exerting a uniform pressure in all layers is finite; i.e., it is compressible in the $q \to 0$ limit. This can be understood considering the fluctuation-dissipation theorem that implies that the
considered responses are proportional to the vanishing $S(q)$ and the nonvanishing $S^{3d}(q, q_z = 0)$ respectively (see Appendix B). This interesting behavior can be explained by the three-dimensional 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 in modifying such hyperuniform behavior in several interesting ways. 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_2Sr_2CaCu_2O_{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 nonasymptotic behaviors, such as the ones arising from the dispersivity of elastic constants and non equilibrium effects, and finite-size effects.

The spontaneously formed hyperuniform magnetic patterns we find may be, for instance, used to affect the properties of a second system which mainly couples to the surface of the superconductor. This kind of surface interactions is indeed possible, and the Fe clumps resulting from magnetic decoration is just probably the simplest example of a system affected by such a two-dimensional magnetic pattern. One could imagine a different second system, for instance, a film with particular physical properties (optical, mechanical, magnetic, or electrical) on top of the superconductor and affected by the presence of the hyperuniform magnetic pattern. This magnetic substrate may also affect transport properties of a second system. If the second system is another superconducting system with vortices and an applied current, an enhanced pinning effect, analogous to the one described in Ref. [5], could be envisaged. Proposing concrete applications of our results is highly speculative at this point and out of the scope of our paper.

On the other hand, the kind of study we have done may offer a different viewpoint for analyzing vortex phases and their transitions, and also to characterizing 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–37], where genuine nonequilibrium 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.

## ACKNOWLEDGMENT

We acknowledge illuminating discussions with C. Olson-Reichhardt, P. Le Doussal, T. Giamarchi, and H. Suderow.

## APPENDIX A: STRUCTURE FACTORS

Here, we demonstrate the validity of Eq. (11), relating the three-dimensional structure factor $S^{3d}(\mathbf{q}, q_z)$ to the two-dimensional structure factor of a single-layer $S(\mathbf{q})$, and make explicit their relation with different susceptibilities at equilibrium. Let us consider a layered medium of linear sizes $L_x$, $L_y$, $L_z = N_z s$ containing $N_v$ vortices, with layer spacing $s$ in the $z$ direction. The vortex density $n(\mathbf{r}, z)$ is considered a continuous function of $\mathbf{r} \equiv (x, y)$ and a discrete function of $z$, such that $\mathbf{r} \equiv (x, y)$ with $0 < x < L_x$ and $0 < y < L_y$, and $z \equiv z_m = ms$, with $m = 0, \dots, N_z - 1$.

The two-dimensional structure factor of a single layer at a given $z$ can be written as

$$S(\mathbf{q}, z) = \frac{1}{n_0} \langle |n(\mathbf{q}, z)|^2 \rangle, \quad (A1)$$

where $\langle \cdots \rangle$ denotes ensemble average, $n_0 = A^2/N_v$ is the average areal density, and $n(\mathbf{q}, z)$ is the Fourier transform of the vortex density in the $x$ and $y$ coordinates,

$$n(\mathbf{q}, z) = \frac{1}{L_x L_y} \int d\mathbf{x} e^{-i\mathbf{q} \cdot \mathbf{x}} n(\mathbf{x}, z). \quad (A2)$$

The discrete transform over the $z$ coordinate yields

$$n(\mathbf{q}, q_z) = \frac{1}{N_z L_x L_y} \sum_z \int d\mathbf{x} e^{-i\mathbf{q} \cdot \mathbf{x} - iq_z z} n(\mathbf{x}, z). \quad (A3)$$

When $N_z \gg 1$, every layer in the bulk is expected to be statistically equivalent, so we can write

$$S(\mathbf{q}, z) \approx \frac{1}{N_z} \sum_{n=1}^{N_z} S(\mathbf{q}, z_m) \equiv S(\mathbf{q}). \quad (A4)$$

Then, using Eq. (A1) we get

$$S(\mathbf{q}) = \frac{1}{N_z n_0} \sum_z \langle |n(\mathbf{q}, z)|^2 \rangle = \frac{1}{N_z n_0} \sum_z \langle n(\mathbf{q}, z) n^*(\mathbf{q}, z) \rangle. \quad (A5)$$

Writing $n(\mathbf{q}, z)$ in terms of its inverse discrete Fourier transform in the $z$ variable, we can write

$$ \begin{aligned} S(\mathbf{q}) &= \frac{1}{n_0} \sum_{q_z, q'_z} \langle n(\mathbf{q}, q_z) n^*(\mathbf{q}, q'_z) \rangle \frac{1}{N_z} \sum_z e^{i(q_z - q'_z)z} \\ &= \sum_{q_z, q'_z} \frac{1}{n_0} \langle n(\mathbf{q}, q_z) n^*(\mathbf{q}, q'_z) \rangle \delta_{q_z, q'_z} \\ &= \sum_{q_z} \frac{1}{n_0} \langle |n(\mathbf{q}, q_z)|^2 \rangle \equiv \sum_{q_z} S^{3d}(\mathbf{q}, q_z). \end{aligned} \quad (A6) $$

In the continuum limit, we have the general relation

$$S(q) \propto \int dq_z S^{3d}(\mathbf{q}, q_z) \quad (A7)$$

mentioned in Eq. (11) in the main text.
Marchetti and Nelson [10] have shown that border effects and physical surface interactions are not relevant for the typical experimentally accessible range of $q$ in magnetic decorations, i.e., $S(q, z = 0) \approx S(q)$. Therefore, decoration patterns are expected to be representative of the behavior in the bulk layers, with $0 \ll z \ll L_z$.

## APPENDIX B: LINEAR RESPONSE

Let us consider a weak periodically modulated external field $\phi(\mathbf{r}, z)$ coupled to the density,

$$ \phi(\mathbf{r}, z) = \frac{1}{L_x L_y N_z} \delta\phi(\mathbf{q}, q_z) \exp(-i[\mathbf{q} \cdot \mathbf{r} + q_z z]), \quad (B1) $$

such that $\delta\phi(\mathbf{q}, q_z)$ is small compared to $k_B T$. The total potential energy is then

$$ \Phi = \sum_z \int_{x,y} dr \phi(\mathbf{r}, z) n(\mathbf{r}, z) = \frac{1}{L_x L_y N_z} \delta\phi(\mathbf{q}, q_z) n(\mathbf{q}, q_z). \quad (B2) $$

This external field will induce a modulation of the local density $n(\mathbf{r}, z)$. In the linear response regime, the equilibrium fluctuation-dissipation theorem states that the local susceptibility $\chi(\mathbf{q}, q_z) \equiv \langle n(\mathbf{q}, q_z) \rangle / \delta\phi(\mathbf{q}, q_z)$ is proportional to the spontaneous density fluctuations [38]

$$ \chi(\mathbf{q}, q_z) = -\frac{n_0}{k_B T} S^{3d}(\mathbf{q}, q_z) \quad (B3) $$

with $S^{3d}(\mathbf{q}, q_z) \equiv \langle n(\mathbf{q}, q_z) n(-\mathbf{q}, -q_z) \rangle / n_0$ being the three-dimensional structure factor defined in Eq. (8).

If the external field is uniform in the $z$ direction (but oscillatory in the plane), the susceptibility

$$ \chi(\mathbf{q}, q_z = 0) = -\frac{n_0}{k_B T} S^{3d}(\mathbf{q}, q_z = 0) \quad (B4) $$

measures, in the $q \to 0$ limit, the compressibility of the vortex lattice to an isotropic pressure from the sides, in the $x-y$ direction. This quantity is analogous to the isothermal compressibility $\kappa_T$ of fluids under hydrostatic pressure.

If instead of a uniform external field we apply the external field only in a given $z = z_0$ layer, $\phi(\mathbf{r}, z) \equiv \phi(\mathbf{r})\delta_{z,z_0}$. If this layer is statistically equivalent to any other layer, the corresponding response $\chi(\mathbf{q}, q_0)$ equals the superposition of all $z$ sum $\chi(\mathbf{q}, q_z)$ over all $q_z$,

$$ \chi(\mathbf{q}, z_0) \equiv \sum_{q_z} \chi(\mathbf{q}, q_z) = -\frac{n_0}{k_B T} S(\mathbf{q}). \quad (B5) $$

In the $q \to 0$ limit, this response function measures the compressibility of a single representative layer to a localized isotropic pressure in the $x-y$ plane. This quantity is thus analogous to the effective two-dimensional isothermal compressibility $\kappa_T$ of the $z = z_0$ subsystem.

## APPENDIX C: ERROR ESTIMATE FOR $\alpha_{eff}$

As described in the main text, by fitting the experimental data in a range [$q_{min}$, $q_{max}$] with the model $S(q) \sim q^\alpha$ we obtain by a least squares procedure the $\alpha_{eff}$ shown in Fig. 7. In order to estimate the error bar, $\Delta\alpha_{eff}$, we consider different sources of uncertainty by using the simple ad hoc formula

$$ \Delta\alpha_{eff} = \sqrt{\sigma_{fit}^2 + (\alpha_{high} - \alpha_{low})^2/4 + \sigma_{low}^2 + \sigma_{high}^2}, \quad (C1) $$

where the following are true:

(i) $\sigma_{fit}^2$ is the least squares error of $\alpha$ in a fit of $S(q) \sim Aq^\alpha$ in the range [$q_{min}$, $q_{max}$].

(ii) By defining $q_M = (q_{max} + q_{min})/2$, $\alpha_{low}$ and $\alpha_{high}$ are the values of $\alpha$ fitted in the ranges [$q_{min}$, $q_M$] and [$q_M$, $q_{max}$], respectively, while $\sigma_{low}^2$ and $\sigma_{high}^2$ are their corresponding least squares errors.

In this way, we incorporate, in a conservative way, uncertainties associated to the fitting range extent [$q_{min}$, $q_{max}$], as well as those given by the least squares fits in the whole range.

[1] S. Torquato and F. H. Stillinger, Local density fluctuations, hyperuniformity, and order metrics, Phys. Rev. E 68, 041113 (2003).

[2] S. Torquato, Hyperuniform states of matter, Phys. Rep. 745, 1 (2018).

[3] M. Florescu, S. Torquato, and P. J. Steinhardt, Designer disordered materials with large, complete photonic band gaps, Proc. Natl. Acad. Sci. USA 106, 20658 (2009).

[4] P. Le Doussal, Novel phases of vortices in superconductors, Int. J. Mod. Phys. B 24, 3855 (2010).

[5] Q. Le Thien, D. McDermott, C. J. O. Reichhardt, and C. Reichhardt, Enhanced pinning for vortices in hyperuniform pinning arrays and emergent hyperuniform vortex configurations with quenched disorder, Phys. Rev. B 96, 094516 (2017).

[6] I. A. Sadovskyy, A. Koshelev, W-K. Kwok, U. Welp, and A. Glatz, Targeted evolution of pinning landscapes for large superconducting critical currents, Proc. Nat. Acad. Sci. USA 116, 10291 (2019).

[7] D. R. Nelson and P. Le Doussal, Correlations in flux liquids with weak disorder, Phys. Rev. B 42, 10113 (1990).

[8] M. C. Marchetti and D. R. Nelson, Dynamics of flux-line liquids in high-$T_c$ superconductors, Phys. C (Amsterdam, Neth.) 174, 40 (1991).

[9] In a realistic finite system, the $z$ dependence must be kept in order to capture the surface-modified vortex interaction due to the stray field interaction of vortex tips near the top $z = L$ and bottom $z = 0$ layers. This effect was shown to appear only at very small $q \ll 1/L$, currently inaccessible in the typical fields of view of decoration experiments. We will thus drop the $z$ dependence, $S(\mathbf{q}, z) \equiv S(\mathbf{q})$, and further assume that the behavior of the top layer is representative of that of the inner layers [10].

[10] M. C. Marchetti and D. R. Nelson, Translational correlations in the vortex array at the surface of a type-II superconductor, Phys. Rev. B 47, 12214 (1993).

[11] E. H. Brandt, The flux-line lattice in superconductors, Rep. Prog. Phys. 58, 1465 (1995).
[12] We are abusing the notation by using the labels $c_{11}$ and $c_{44}$ for the elastic modulii in the solid as well as in the liquid phases. Nevertheless, in the small wave-vector limit, the corresponding elastic modulii are expected to be close [11].

[13] D. R. Nelson and V. M. Vinokur, Boson localization and cor-
related pinning of superconducting vortex arrays, Phys. Rev. B
**48**, 13060 (1993).

[14] G. Blatter, M. V. Feigel'man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Vortices in high-temperature superconductors, Rev. Mod. Phys. **66**, 1125 (1994).

[15] T. Nattermann and S. Scheidl, Vortex-glass phases in type-II superconductors, Adv. Phys. **49**, 607 (2000).

[16] L. Radzihovsky and E. Frey, Kinetic theory of flux-line hydro-dynamics: Liquid phase with disorder, Phys. Rev. B **48**, 10357 (1993).

[17] T. Giamarchi and P. Le Doussal, Elastic theory of flux lattices in the presence of weak disorder, Phys. Rev. B **52**, 1242 (1995).

[18] T. Klein, I. Joumard, S. Blanchard, J. Marcus, R. Cubitt, T. Giamarchi, and P. Le Doussal, A Bragg glass phase in the vortex lattice of a type II superconductor, Nature (London) **413**, 404 (2001).

[19] J. Kim and S. Torquato, Effect of imperfections on the hyperuniformity of many-body systems, Phys. Rev. B **97**, 054105 (2018).

[20] J. F. Joanny and P. G. de Gennes, A model for contact angle hysteresis, J. Chem. Phys. **81**, 552 (1984).

[21] P.-G. de Gennes, *Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves* (Springer, Berlin, 2003).

[22] D. Carpentier and P. Le Doussal, Glass phase of two-dimensional triangular elastic lattices with disorder, Phys. Rev. B **55**, 12128 (1997).

[23] M. Konczykowski, C. J. van der Beek, A. E. Koshelev, V. Mosser, M. Li, and P. H. Kes, Vortex matter in $Bi_2Sr_2CaCu_2O_{8+\delta}$ with pointlike disorder, J. Phys. Conf. Ser. **150**, 052119 (2009).

[24] J. Aragón Sánchez, R. Cortés Maldonado, N. R. Cejas Bolecek, G. Rumi, P. Pedrazzini, M. I. Dolz, G. Nieva, C. J. van der Beek, M. Konczykowski, C. D. Dewhurst, R. Cubitt, A. B. Kolton, A. Pautrat, and Y. Fasano, Unveiling the vortex glass phase in the surface and volume of a type-II superconductor, arXiv:1905.06908 [Commun. Phys. Nat. (to be published)].

[25] Y. Fasano, J. Herbsommer, and F. de la Cruz, Superficial periodic pinning induced by bitter decoration applied to the study of vortex structure nucleation and growth, Phys. Status Solidi B **215**, 563 (1999).

[26] M. Menghini, Y. Fasano, F. de la Cruz, S. S. Banerjee, Y. Myasoedov, E. Zeldov, C. J. van der Beek, M. Konczykowski, and T. Tamegai, First-Order Phase Transition from the Vortex Liquid to an Amorphous Solid, Phys. Rev. Lett. **90**, 147001 (2003).

[27] F. Pardo, A. P. Mackenzie, F. de la Cruz, and J. Guimpel, Effect of the reversibility region on the low-temperature vortex structure imaged by bitter magnetic decoration, Phys. Rev. B **55**, 14610 (1997).

[28] Y. Fasano, M. De Seta, M. Menghini, H. Pastoriza, and F. de la Cruz, Commensurability and stability in nonperiodic systems, Proc. Natl. Acad. Sci. USA **102**, 3898 (2005).

[29] N. R. Cejas Bolecek, M. I. Dolz, H. Pastoriza, M. Konczykowski, C. J. van der Beek, A. B. Kolton, and Y. Fasano, Excess of topological defects induced by confinement in vortex nanocrystals, Phys. Rev. B **96**, 024507 (2017).

[30] A. B. Kolton, A. Rosso, and T. Giamarchi, Nonequilibrium Relaxation of an Elastic String in a Random Potential, Phys. Rev. Lett. **95**, 180604 (2005).

[31] P. Kim, Z. Yao, C. A. Bolle, and C. M. Lieber, Structure of flux line lattices with weak disorder at large length scales, Phys. Rev. B **60**, R12589 (1999).

[32] M. I. Dolz, Y. Fasano, N. R. Cejas Bolecek, H. Pastoriza, V. Mosser, M. Li, and M. Konczykowski, Size-Induced Depression of First-Order Transition Lines and Entropy Jump in Extremely Layered Nanocrystalline Vortex Matter, Phys. Rev. Lett. **115**, 137003 (2015).

[33] J. N. Avraham, B. Khaykovich, Y. Myasoedov, M. Rappaport, D. E. Shtrikman, H. C. J. Feldman, T. Tamegai, P. H. Kes, M. Li, M. Konczykowski, C. J. van der Beek, and E. Zeldov, Inverse melting of a vortex lattice, Nature (London) **411**, 451 (2001).

[34] F. Pardo, F. de la Cruz, P. L. Gammel, E. Bucher, and D. J. Bishop, Observation of smectic and moving-Bragg-glass phases in flowing vortex lattices, Nature (London) **396**, 348 (1998).

[35] P. Le Doussal and T. Giamarchi, Moving glass theory of driven lattices with disorder, Phys. Rev. B **57**, 11356 (1998).

[36] L. Balents, M. C. Marchetti, and L. Radzihovsky, Nonequilibrium steady states of driven periodic media, Phys. Rev. B **57**, 7705 (1998).

[37] S. Scheidl and V. M. Vinokur, Driven dynamics of periodic elastic media in disorder, Phys. Rev. E **57**, 2574 (1998).

[38] J. L. Barrat and J. P. Hansen, *Basic Concepts for Simple and Complex Liquids* (Cambridge University Press, Cambridge, UK, 2003).