Observation of a gel of quantum vortices in a superconductor at very low magnetic fields

A gel consists of a network of particles or molecules formed for example using the sol-gel process, by which a solution transforms into a porous solid. Particles or molecules in a gel are mainly organized on a scaffold that makes up a porous system. Quantized vortices in type II superconductors mostly form spatially homogeneous ordered or amorphous solids. Here we present high-resolution imaging of the vortex lattice displaying dense vortex clusters separated by sparse or entirely vortex-free regions in $\beta$-Bi$_2$Pd superconductor. We find that the intervortex distance diverges upon decreasing the magnetic field and that vortex lattice images follow a multifractal behavior. These properties, characteristic of gels, establish the presence of a novel vortex distribution, distinctly different from the well-studied disordered and glassy phases observed in high-temperature and conventional superconductors. The observed behavior is caused by a scaffold of one-dimensional structural defects with enhanced stress close to the defects. The vortex gel might often occur in type-II superconductors at low magnetic fields. Such vortex distributions should allow to considerably simplify control over vortex positions and manipulation of quantum vortex states.


INTRODUCTION
Quantized vortices in superconductors arrange spatially in structures that bear some similarities with atomic or molecular arrangements. The main difference is that vortices consist simply of single quantized fluxes with repulsive interactions, whereas atoms and molecules have numerous degrees of freedom and allow for bonding. Nevertheless, interactions among vortices can also be varied because these are related to the properties of the superconducting material hosting them 1,2 . Thus, solid, liquid or disordered glass, have been observed in vortex matter [3][4][5] . However, a gel with inhomogeneous vortex density distribution caused by a network has not yet been reported. Experiments made in superconductors at low magnetic fields have until now unveiled vortex clustering in hexagonal lattices due to attractive interactions or the intermediate mixed state in single crystals or polycrystalline vortex arrangements and glassy phases in presence of disorder [5][6][7][8][9][10][11][12][13][14][15][16] . The vortex density remains spatially ho-mogeneus in all these phases. Here we image the vortex configurations at very low magnetic fields in the type-II superconductor β−Bi 2 Pd. We find vortex clusters whose distribution has characteristics specific to gels, such as a wide distribution of intervortex separation, covering widely different distances that diverges when decreasing the magnetic field and is characterized by multiple fractal exponents.
The new vortex configuration is dominated by vortices arranging along linear defects, leaving isolated vortices well separated from their neighbors. Such a configuration should be helpful to manipulate topological quantum vortex states, such as those that might arise in the framework of the recently proposed topological superconductivity in β−Bi 2 Pd and in iron pnictide superconductors [17][18][19][20] .
β-Bi 2 Pd crystallizes in a tetragonal structure and becomes superconducting below T c ≈ 5 K 21,22 . Zerotemperature critical magnetic fields are H c1 = 225 G and H c2 = 6000 G (yielding coherence length ξ(0) ≈ 23.5 nm and penetration depth λ(0) ≈ 132 nm). The upper critical field anisotropy is small, of order of 10% 20,22,23 . Single crystals of this material can be easily cleaved, leaving large stable and atomically flat terraces, where a fully formed s-wave superconducting gap is observed over the entire surface 22 . Here we undertook a thorough examination of vortex states using SQUID-on-tip (SOT) microscopy 24,25 and magnetic force microscopy (MFM) 26 in a large range of magnetic fields, from 1 to 600 G.

I. EXPERIMENTAL
Sample preparation β-Bi 2 Pd single crystals were grown as described in Ref. 22. The β-Bi 2 Pd phase is thermodynamically unstable at room temperature. Crystals need to be quenched from the high temperature β-Bi 2 Pd phase by immersion in cold water. Structural characterization gives high quality single crystalline samples. The superconducting transition (T c = 5 K) is sharp in all thermodynamic measurements 21,23 . However, the residual resistivity, of 18 µΩcm is quite high, pointing out that, during the quench, a sizeable amount of defects has been left in the sample. The electronic mean free path is of ≈ 15 nm (from residual resistivity 22 ) and the thickness of the sample is of about 0.1 mm.
We prepare the sample surface by cleaving with scotch tape and obtain atomically flat surfaces suitable for scanning probe imaging. We provide a detailed description of defects and of the surface topography in Appendix A. The magnetic field is applied perpendicular to the surface using superconducting coils. The plate like single crystals used here give demagnetizing factors close to one and we usually measure in field-cooled conditions, heating the sample above T c after each field change.

Magnetic Force Microscopy
Magnetic force microscopy (MFM) measurements have been performed in a commercial Low-Temperature SPM equipment from Nanomagnetics Instruments, working in the 300 K -1.8 K temperature range. The microscope is inserted in a superconducting coil 27 . Simultaneous atomic force microscopy (AFM) and magnetic MFM images are obtained in dynamic mode, oscillating the cantilever at its resonance frequency, and using the so-called retrace mode for magnetic detection. In this mode, the tip scans twice over the sample surface at two different distances. A first scan is performed at distances of few nm to extract the AFM profile. Then, the tip is retracted by a selected distance (∼ 100 nm in this work) and the system repeats the same profile obtained in the first scan. Phase shift of the cantilever oscillation during this retrace scan, caused by the long range magnetostatic interaction, is used to build the MFM image. We use com-mercial MFM tips from Nanosensors (PPP-MFMR) and magnetize these prior to the measurement by applying a magnetic field of 1500 G at 10 K.

Scanning SQUID-on-tip Microscopy
Scanning SOT images provides high spatial resolution magnetic imaging 28,29 reaching single-spin sensitivity 24,25 and enabling detection of sub-nanometer and ultrafast vortex displacements 30,31 . The SOT used in this work had an effective diameter of 260 nm, 160 µA critical current and white flux noise of around 800 nΦ 0 Hz −1/2 above a few hundred Hz. The SOT was mounted in a home-built scanning probe microscope with a scanning range of 30 × 30 µm 2 and read out using a series SQUID array amplifier 32 . The SOT images in this manuscript show an area of 20 × 20 µm 2 with a pixel size of 100 × 100 nm 2 . The acquisition time was 100 s per image. All measurements were performed at 4.2 K in an open loop mode, at constant height of ∼300 nm above the sample.

Ginzburg-Landau numerical simulations
The numerical simulations were performed within the Ginzburg-Landau (GL) formalism, taking into account the spatially inhomogeneous microscopic parameters. In their stationary form, the dimensionless GL equations for the superconducting order parameter Ψ and the magnetic vector potential A read: where f (r) = 1−t(r) 2 1+t(r) 2 and g(r) = 1 (1+t(r) 2 ) 2 contain the temperature dependence (proposed empirically in Ref. 33, and proven to describe well the standard experimental observations) and the spatially-dependent critical temperature through t(r) = T /T c (r). In Eqs. (1-2), κ(0) = λ(0)/ξ(0) is the GL parameter at zero temperature, the distance is measured in units of the coherence length ξ(0), the vector potential A in c /2eξ(0), and the order parameter Ψ is scaled to its value at zero field and temperature. We apply a finite-difference representation for the order parameter and the vector potential on a uniform 2D (x,y) Cartesian space grid, with spatial resolution of 0.1ξ(0). We use periodic boundary conditions on the simulated rectangular cell W x × W y , in the form A(r + b i ) = A(r) + ∇η i (r), and Ψ(r + b i ) = Ψ exp(2πiη i (r)/Φ 0 ), where b i=x,y are the supercell lattice vectors, and η i is the gauge potential. Since the sample is exposed to a homogeneous perpendicular magnetic field H = He z ) we employ the Landau gauge A ext = Hxe y for the external vector potential and η x = HW x y while η y = 0. The value of applied magnetic field must always match the integer number of flux quanta (vortices) in the simulated cell, as stipulated by the virial theorem 34 . In most simulations we deliberately choose W x = √ 3W y that should favor a periodic triangular lattice of vortices, so that every departure from an Abrikosov lattice pattern is easily seen.

RESULTS
In Fig. 1(a) we show a 7 × 7 µm 2 AFM image taken at 2 K. We observe flat surfaces interspersed with steps of 10 nm height or larger. The vortex lattice shown in Fig. 1(b) is measured simultaneously with MFM at 300 G. Fig. 1(c) shows an optical image of the SOT near the surface. In Fig. 1(d-i) we show image of the local magnetic field B(x, y) acquired using the SOT between 2 and 50 G taken in field-cooled conditions at 4.2 K. At low magnetic fields, Fig. 1(d,e), many vortices are located along one dominant line, while, in the rest of the frame, we observe large vortex-free regions. When increasing the magnetic field, vortices cluster along lines and the vortex lattice is formed in between ( Fig. 1(f-i)).
Note that for vortices located along lines, the value of the magnetic field at the vortex center is smaller than the value we find for vortices located far from lines. This is nicely visible at low fields when vortices are well separated and do not overlap. For example, in Fig. 1(d) we see that vortices located along the main linear defect visible in the image (from upper left to middle right) are not as bright as those located elsewhere in the image. This remains in the whole range of magnetic fields (in Fig. 1(f-i) there is strong vortex overlap, which produces linear vortex clusters that appear very bright). When we fit the spatial profile of vortices along lines in the images at low fields we find values for λ close to 340 ± 40 nm, which is about twice the penetration depth found in experiments measuring bulk field penetration with Hall probes 23 . Further details on the fitting procedure are given in the Appendix B.
We have made numerical simulations of vortex behavior close to linear defects using Ginzburg-Landau theory. We parametrize the simulations according to the values of λ measured, that are translated into small changes of T c close to the defects and weak disorder (i.e. a randomized value of T c ). The behavior close to the linear defect is modeled through a parabolic recovery of T c at a distance of 2 µm away from the step. This simulates the pinning potential of the linear defect. The obtained vortex configurations over a 10 × 20 µm 2 area are shown in Fig. 2(a-c), for three values of applied magnetic field. This indeed captures the evolution seen in the images. Namely, vortices first occupy locations along the defect, where superconductivity is suppressed, and their magnetic field is weaker than that of vortices found away from the defect. The accumulation of vortices at the defect strengthens repulsion of other vortices in the sample, so that a noticeable vortex-free band is formed between vortices of different brightness in the images. The vortex free band diminishes with increasing magnetic field. Similar vortex free bands are observed in the experiment ( Fig. 1(f-i)).
The more striking result appears at the smallest magnetic fields. There, vortices are located along lines separated by large vortex free areas ( Fig. 1(d)). We have Delaunay triangulated images and calculated intervor- , where LD is the lateral size of the defect. Note that we allow for a slow decay when leaving the linear defect. The position of the defect is marked by a red dashed line. Panels (a-c) show a field of view of 10 × 20 µm 2 , for applied magnetic fields of 2 G, 5 G, and 15 G, respectively. We plot the magnetic induction at a height of 300 nm above the sample surface. This gives a span of ≈ 15 G (a), 16 G (b) and 18 G (c) in each image. Inset in (a) shows comparative magnetic profiles of two selected vortices, marked in the main panel by green and magenta lines. tex distances. In Fig. 3(a) we show intervortex distance histograms and in Fig. 3(b) the standard deviation of intervortex distances. First of all, let us remark that all distance histograms ( Fig. 3(a)) show just a single peak located at a 0 = ( 4 3 ) 1/4 Φ0 This means that the average intervortex distance follows the Abrikosov rule for a triangular lattice. The magnetic field behavior of the position of this peak is shown in the inset of Fig. 3(b) and we can see that its value coincides with a 0 also at the weakest magnetic fields. This includes situations where vortices are strongly clustered with large voids in between clusters. At the weakest magnetic fields, the histograms of Fig. 3(a) are extremely broad and the distances between vortices become widespread. We observe intervortex distances d up to twice a 0 and many vortices located much closer to each other than a 0 . The standard deviation of the histograms as a function of the magnetic field ( Fig. 3(b)) diverges as 1/ √ µ 0 H below about 80 G. We compare our results with distance histograms obtained in the same range of magnetic fields in a cuprate superconductor, Bi 2 Sr 2 CaCu 2 O 8 in the cases of pristine samples and crystals with a strong and dense distribution of pinning centers (columnar defects) 35 . In both cases the vortex density is much more homogeneous than in β-Bi 2 Pd and the standard deviation remains constant when decreasing the magnetic field (black and violet points in Fig. 3(b)).

DISCUSSION
All our experiments are in field-cooled conditions, so during cooling, vortices are nucleated inside the sample. Our data show that there are locations with vortex pinning determined by linear defects, which provide a scaffold on which vortices arrange and areas in between with just a few vortices. Other mechanisms that can lead to a modified vortex distribution, different than pinning, are discussed in Appendices C and D and can be definitely excluded in β−Bi 2 Pd.
Vortices located on the scaffold have a different λ. We consider stress as a possible mechanism to explain the spatial variation of λ. Stress builds up close to defects and modifies the superconducting length scales (and the Ginzburg-Landau parameter), producing an effective interaction between the crystal lattice and the vortex lattice 36,37 . Recently, it has been shown that stress induced intervortex interaction can lead to square vortex lattice in tetragonal superconductors 38 . Such a coupling between crystalline elasticity and superconductivity can be treated using the dependence of T c with pressure, dT c /dP 36,[39][40][41] . Generally, with dT c /dP > 0, vortices are repelled from places with internal stress and the opposite occurs with dT c /dP < 0. The pressure dependence of T c in β-Bi 2 Pd was measured in bulk samples 42,43 , giving dT c /dP = −0.025 K/kbar. T c,bulk = 5 K for the unstressed region and our SOT measurements are performed at T = 4.2 K. Thus, a variation in T c of a fraction of a degree is sufficient to change penetration depth by a factor ∼2. If we consider that λ(T ) = λ(0)/ 1 − (T /T c ) 4 , we estimate T c,def ect = 4.38 K. According to the study of Ref. 42, this T c corresponds to a local pressure of ∼ 20 kbar. This subtle local change in T c is very difficult to detect (we provide further details in the Appendix E).
To analyse further the vortex distributions, we calculate the elastic energy associated to pairs of vortices, F , at different locations in our images. We compare the re-sult for vortices located at a step with the elastic energy for pairs of vortices far from the steps. To this end, we use F = φ 2 0 4πµ0λ 2 log(κ) + φ 2 0 4πµ0λ 2 K 0 (d/λ) for the free energy per unit length of two vortices interacting with each other at a distance d 44 . The first term comes from the energy of superfluid currents, giving the line tension of the vortex, and the second term denotes the interaction energy between vortices. K 0 is the order zero modified Bessel function of the second kind. We then calculate F for vortices far from defects using the bulk λ. For vortices at the steps we use λ measured close to defects. This is an approximation since, during a field cooling procedure, the vortex lattice can form closer to T c where λ is significantly larger 45 . However, as we take SOT images at quite high temperature and with an increased λ, we can still gain useful insight with these approximations. We find that the difference in free energy between the two situations is δF ≈ 7 × 10 −12 Jm −1 . This result is independent of the intervortex distance, for fields of order of a few tens of G or less. Below ≈ 50 G the intervortex distances vary from 0.5 to 4 µm and the second term of the interaction energy remains negligible with respect to the first term. Thus, at low magnetic fields, the intervortex interaction essentially vanishes and vortices behave as nearly isolated entities.
For fields above ≈ 50 G, the vortex lattice density increases and the previous two-vortex interaction approximation is no longer valid. We take into account the interaction with the first six neighbours arranged in a hexagonal lattice. We can use the same equation as before, but with an interaction term of 6× The difference in energy between six vortices close to a defect obtained using the two different values of λ considered here changes with the intervortex distance. We find that when vortices at the linear defect are closer than about 250 nm, it is no longer energetically favorable to add new vortices in there, so that linear defects saturate and a homogeneous vortex lattice is formed over the whole sample. The field at which we change regime is consistent with the experiment.
Thus, our direct imaging of the evolutionary patterns of superconducting vortex configurations in β-Bi 2 Pd in a broad range of magnetic fields shows that linear pinning centers produce a scaffold that cluster vortices along lines and leave large vortex free areas. At higher fields, when the intervortex distance is much smaller than the separation between defects, a triangular lattice forms within the vortex free areas observed at low magnetic fields.
Another striking insight from our measurements is that the standard deviation of the distribution of vortices diverges for very small vortex densities as the inverse of the square root of the magnetic field. The flux is quantized at each vortex, so that we can assume that B = N A Φ 0 within the fields of view studied here (N is the total amount of vortices and A the area) and in agreement with the magnetic field dependence of the intervortex distance following the triangular lattice expression (in- In the inset we show intervortex distance vs the magnetic field for the same images. set in Fig. 3(b)). In case of a purely one-dimensional row of vortices of length L, the average magnetic field along the row goes like B ∝ N L Φ 0 and the intervortex distance diverges as d 1D ∝ Φ0 B , which gives intervortex distances d 1D much smaller than a 0 if we consider vortex rows in two-dimensional areas. A scaffold of onedimensional lines distributed within a two-dimensional area A presents a standard deviation of intervortex distances normalized by a 0 which increases when decreasing the magnetic field and diverges as 1/ √ µ 0 H, as we observe experimentally.
On the contrary, if we consider a homogeneous random distribution of point pinning centers, the vortices are pinned randomly. As the intervortex distance varies with the vortex density, the deviation of intervortex distances over the average value remains constant when decreasing the magnetic field. This explains the difference between the data in the vortex configurations observed in cuprates and the divergence observed in β−Bi 2 Pd in Fig. 3(b). Actually, we can understand the behavior in β-Bi 2 Pd as 1D dominated pinning at very low magnetic fields, which becomes usual 2D or 3D dominated pinning at higher magnetic fields.
In addition, the difference between the vortex distributions observed in β−Bi 2 Pd and in the cuprates is that in the former there are strong spatial variations in the vortex density. Thus, there is a parallel with structural glasses and gels that is straightforward. In a similar way that gels are different from glasses because they are created from spatially highly inhomogenoeus amorphous particle arrangements, the vortex gel differs from other vortex configurations, such as the Bose glass or the Bragg glass, by a highly inhomogeneous vortex distribution. There are multiple forms of gels and all share the presence of a scaffold that holds a network of particles. Accordingly, statistical properties of particle distribution and sizes are varied and mimic the properties of the scaffold. Here we have very simple particles, all equal to each other and characterized by carrying exactly a quantized flux, whereas most gels are formed by polymers or other complex structures. The scaffold is also relatively simple, linear defects on a crystalline lattice. Therefore, In the inset we show the generalized dimension Dq as a function of the set of exponents q. Note that the curves strongly vary for low magnetic fields, but remain centered around 2 for high magnetic fields and for the results in the cuprate superconductor. Details on the calculation of multifractal parameters are provided in Appendix F.
the statistics of the vortex gel reveals the distribution of these linear defects in the sample. Often, gels form self-similar or fractal structures. It has been shown that the fractal dimension remains constant with density, unless the density is significantly decreased, in which case the rugosity of constituents increases with decreasing density and the corresponding fractal dimension too 46 . It is thus interesting to ask whether the distribution of vortices at low magnetic fields also shows similar features. We have calculated the generalized dimension D q and the multifractal spectrum f (α) (see Appendix F for details on the calculation). The results are shown in Fig. 4. Images showing triangular or disordered lattices provide distribution of fractal subsets centered at α = D q = 2. When we start observing strong variations in the vortex density, the lacunarity increases, which leads to f (α) that are much broader and whose maximum deviates from 2. D q also increases for small values of the multiscaling exponent q. Thus, the vortex distributions at small magnetic fields are multifractal, with a probability of fractal subsets that strongly increases when decreasing the magnetic field, leading to a widening of f (α). On the other hand, the maximum of f (α) and the value of D q for small q is larger than two and increases when the magnetic field is decreased. As we discuss in Appendix A, we find similar fractal properties in images of the linear defects of β−Bi 2 Pd. We further show in Appendix G the Voronoi tesselation of vortex images, showing that vortices cluster along the lines and that vortices in between lines are practically isolated.
The outstanding feature of the vortex gel is that the standard deviation of vortex positions diverges for decreasing magnetic fields, suggesting that the magnetic field penetration is neither a Meissner state interspersed with normal areas (the intermediate state 7 ), nor an intermediate mixed state with hexagonal vortex clusters (as in pure Nb 16,47 ), nor the disordered vortex lattices observed in high T c cuprate superconductors 1,2 . Other arrangements with strong variations in vortex densities are found in conformal crystal arrays of pinning centers and when periodic and random pinning potentials are formed [48][49][50] . In the latter case simulations show that the presence of square crystal and hexagonal vortex lattices can lead to one-dimensional fractal structures, in a small range somewhat below the pinning strength range where the combination of random pinning and the crystal lattice is sufficiently strong to form vortex chains 48 . Our work shows that this actually occurs in a wide range of magnetic fields in presence of one-dimensional pinning centers, leading to vortex distributions with structural properties typical of a gel.
Vortex positions are determined by pinning, and we do not see indications of thermal fluctuations playing a role in the vortex location. Thus, the transition from the normal Abrikosov lattice into the vortex liquid (which can only appear in an extremely small temperature range in β−Bi 2 Pd) and the vortex gel should be continuous and be produced by dynamical arrest. Thus, it depends rather on the diffusion process of vortices along the pinning centers than on equilibrium properties of the superconductor. As such, it should appear in all crystalline superconductors with pinning centers having a spatially inhomogeneous distribution at sufficiently low magnetic fields. Vortex lattices showing voids have been indeed observed a few times in some Fe based materials and in the cuprate superconductors [51][52][53] . The inhomogeneous vortex distributions occur either due to twin boundaries in orthorhombic superconductors 52 , or to highly inhomogeneous pinning centers 53 . Data were taken at single values of the magnetic field, sometimes two orders of magnitude below the ones we discuss here 51 . Thus, the divergence of the standard deviation and the multifractal distribution could not be followed with magnetic field nor identified. Nevertheless, the highly inhomogeneous vortex distributions found there suggests that the vortex gel might occur in many superconductors at low magnetic fields.

CONCLUSIONS AND OUTLOOK
Our measurements illuminate the interplay of geometric defects and crystalline stress at very low magnetic fields in superconductors, which turns out to be much more varied than previously thought. We reveal novel vortex arrangements at low magnetic fields, governed by organizing principles that combine pinning centers as scaffolds, screening and intervortex interactions.
Screening has been considered to design mechanical resonators that are isolated from the environment for quantum circuits, to produce a scaffold of traps for cold atoms close to the surface of a superconductor or to design magnetic cloaks [54][55][56] . Recent work shows that Majorana modes might be present in superconducting vortices 57 . This calls for methods to manipulate vortices and entangle their states to search for non-Abelian statistics 18,58 . Superconductors considered are, among others, Fe based materials whose structural properties and superconducting parameters are similar to β−Bi 2 Pd and β−Bi 2 Pd itself 17-20,59,60 .
As we discuss in Appendix H, there is mounting evidence for the presence of triplet correlations in β−Bi 2 Pd, particularly from experiments that are probing the surface. The question then arises of why we did not observe a half-integer flux quantum in our experiment. The answer might be that triplet correlations could form only close to the surface. Magnetic vortices, however, are threads through the bulk of the material and the surface properties might be masked by the flux quantization from the bulk. It would be very interesting to repeat our experiments in thin films or very thin layers of β−Bi 2 Pd. The appearance of half-integer flux quanta is, in view of the recent reports 19,61,62 , quite likely. Vortices might then carry a Majorana fermion, as proposed in Refs. 63-67. In that case, the large intervortex distances found in the range of magnetic fields we study here should be very helpful to facilitate manipulation of vortices and braiding experiments. For example, a couple of vortices located in between lines could be easily moved around each other, as proposed for instance in Refs. [68][69][70] .
Vortex manipulation devices remain indeed very difficult to realize in dense vortex lattices at high magnetic fields. It is fundamental to have well separated and isolated vortices to be able to manipulate and entangle vortices (see Appendix B showing isolated vortices in between lines and the evidences for unconventional superconductivity in Appendix H). The vortex gel produces intrinsically areas with flux expulsion and flux concentration. Our results show that vortices are nearly independent to each other at very low magnetic fields and that their position is locked to the defect structure in the sample. This suggests that vortices in this field range are also highly manipulable, much more than in usual hexagonal or disordered vortex lattices.

ACKNOWLEDGMENTS
We acknowledge support, discussions and critical reading of the manuscript from Eli Zeldov, who also devised and set-up the SOT system. We also acknowledge critical reading and suggestions of Vladimir The bonds in the tetragonal structure of β-Bi 2 Pd are such that the surface is most likely formed by the square lattice of Bi atoms 71 . The cleaving plane is thus very well defined and lies perpendicular to the c-axis. There are no indications from van der Waals like bonds as in transition metal dichalcogenides 72 -this material has well established three-dimensional electronic properties. Nevertheless, it is a fact that it can be easily cleaved using scotch tape 22,23,61,62 . Cleavage occurs without any residues, as thick sheets of the material are removed when cleaving. The obtained surface is shown in a optical and scanning electron microscope (SEM) images in the Fig. 5. The surface is very shiny and has features which are important for the results discussed in the main text.
Cleaving or breaking of a surface occurs through the establishment of a fracture or crack at a few places close to the edges of the sample. The fracture then propagates as a crack front through the whole sample. The action on the sample during fracture consists of tear, shear and compressive forces (called mode I, II and III fracture, respectively, see Fig. 5(a) [73][74][75] ). If the action would occur only along the c-axis, just tear forces that separate layers would be active. However, the competition between elastic energy and surface energy is in-plane anisotropic, leading to crack behavior that depends on the in-plane properties of the material. This occurs irrespective of the anisotropy of the in-plane vs out-of-plane crystalline structure, and is even observed in two-dimensional van der Waals materials 76 . In three-dimensional crystals this issue is even more important. The sample is of course not perfect and the crack front encounters defects such as small angle grain boundaries 77,78 . In a cleaving process using scotch tape, shear forces appear easily, because scotch is highly deformable. Shear forces change the direction of crack propagation away from high symmetry crystalline lines. In addition, fracture produces a release of stress that might have been left over during crystal growth (for example by the presence of small temperature gradients). The process of releasing this load is influenced by defects and imperfections in the crystal. All this leads to regions with alternating compressive and tensile stress and causes a twist action (combining shear and tear, modes II and III, see Fig. 5(a)) that might well propagate far below the surface and influence large areas of the crystal 73 .
The features produced by this twist action are called twist hackles and must not necessarily follow a crystalline direction. They rather run parallel to the crack propagation direction. In β-Bi 2 Pd we observe features that can be associated to twist hackles. We mark a few of such features by yellow dashed lines in Fig. 5(b)-(d). Close to the sample edges, twist hackles have a strong tendency to start or arrive to the end of the sample at an angle to the surface, following the crack propagation direction, as we observe in the images (Fig. 5(d)). The crack propagation direction is also influenced by the direction of the crystalline axis. Thus, there are also a number of linear structures at 45 or so degrees to the twist hackles that can be associated to crystalline axes (we mark a few by red dashed lines in the Fig. 5(b), (e)). There are furthermore linear features perpendicular to all them. Thus, the images show resulting from twisting efforts produced during fracture. These efforts might produce enough stress to influence locally the superconducting properties.
Let us note that we also identify large wrinkles on the surface (Fig. 5(b)). The wrinkles appear close to very large defects (broken or open features in the Fig. 5 (a)). A closer analysis using SEM reveals a large number of separated layers close to wrinkles ( Fig. 5(g)). Generally, step edges appear strongly marked in SEM images (Fig. 5(fh)), suggesting that some of parts of the sample separate as layers (green arrows in Fig. 5(g-h)). All this supports the presence of large twisting efforts during fracture.
Our experiments are made close to the center of the sample, in locations showing no large wrinkles and the tip was carefully positioned away from optically visible defects. So that we are far from wrinkles produced during cleaving (red arrows, Fig. 5(b)). Instead, we perform our experiments close to linear structures due to twist hackles.
It is further interesting to search for similarities in the images showing linear features on the surface and vortex lattice images. To this end, we have taken SEM and optical images of the surface of a similar size as vortex lattice images and calculated their fractal properties (see Appendix F for details). The images of the surface show the linear features that act as a scaffold for vortex pinning, so that we can expect some relation among them. Unfortunately, we cannot analyze exactly the same field of view, because SOT images do not provide surface topographic scans. However, we have found fields of view which are very similar than the ones in vortex lattice imaging, both in terms of density of defects, orientation and distribution. Another caveat is that SEM and optical images provide a continous grey scale corresponding to steps or areas that have a large amount of defects, wheras vortex lattice images are a set of pixelized points (being one at a vortex position and zero elsewhere) that are not joined together. We have thus made a comparison for the image with the largest vortex density, which is the image shown in Fig. 1(i) and also increased the contrast of the SEM and optical image by squaring the grey scale of the image. The result is shown in Fig. 6. We see that we do obtain the same set of fractal dimensions in vortex lattice images and in the images made using SEM and an optical microscope. Thus, there is a clear relationship between the scaffold of defects and the vortex lattice distribution.

APPENDIX B: FITTING PROCEDURE OF SOT SCANS
There are two unknown parameters that are essential to obtain a meaningful fit to a magnetic field profile obtained with SOT.
One is the distance between the sample and the SOT d SOT and the other is the SOT transfer function that converts the measured current that flows in the SOT into a magnetic field (dI SOT /dB). The latter is normally determined by measuring the response of the SOT to a known applied magnetic field. Here, the presence of the superconducting sample partially screens the applied field. It is therefore somewhat hard to know which part of the applied field is screened by the sample. We thus assume that vortices far from defects, those that appear brightest in the images, have a flux of Φ 0 and λ = 186 nm 23 . Then, we model the field B(x, y) by a monopole located at λ below the surface i.e. at a distance λ+d SOT from the SOT. All vortices far from the defects that appear bright in Fig. 1(d)-(e) are fitted with fitting parameters d SOT and dI SOT /dB. Averaging results in different vortices we find that d SOT = 270 nm. We also obtain the magnetic field scale in the images by obtaining dI SOT /dB. We considered nine bright and ten less pronounced vortices from Fig. 1(d)-(e).
Let us now discuss two representative examples, one showing bright spot and another one a less pronounced spot in the SOT images. The magnetic field profile B(r) across each vortex is shown in the Fig. 7, together with two different fits. We extract two profiles in orthogonal directions for the vortices on the defect to show that the vortices are isotropic. We obtain similar fits in both direction when vortices are well separated.
We first assumed that the flux is smaller for a vortex located at a defect. To model this situation a fixed λ = 186 nm is considered while the magnetic flux is left as a free parameter. The result is shown in Fig. 7 (green lines). Bright vortices are well described by that value of λ and a value for the flux obtained from the fit very close to the flux quantum, at 1.06Φ 0 . The magnetic field profile of the vortex on the defect is however not well reproduced.
We then assumed that Φ 0 is the same for all vortices and that lambda differs at different locations. The results are shown in Fig. 7(black lines). Here, vortices far from defect yield 172 nm for λ, which is very close to the reported value 23 . The vortex close to the defect yields This leads to a constant Φ 0 and is much more consistent with the experiments. We thus conclude that λ is not homogeneous over the sample and is modified at defects in β−Bi 2 Pd. Vortex field profiles Bz(r) along vortices far from a defect (blue points) and close to a defect along the defect (red points) and along the direction perpendicular to the defect (black points). Profiles are taken from Fig. 1(d) along the white lines. The fit to the vortex profile leaving the magnetic flux φ as a free parameter is shown in blue. Vortices far from defects provide an excellent fit with 1.06 Φ0 where Φ0 is the flux quantum. However, close to the defects the fit does not provide a good account of the experimental data. By contrast, a fit of the same vortex profiles using the penetration depth λ as a free parameter leads to orange lines, which much better account for the experimental data in all situations. We find λ = 172 nm (which is of order of ∼ λ measured using Hall probes, as discussed in the main article) for vortices far from a defect and λ = 360 nm close to the average value found for vortices at a defect. extend to our observations. We also recall the discussions of possible two-gap superconductivity in β-Bi 2 Pd 17 and that vortex clustering, stripes and vortex free areas at low magnetic fields have all been observed in the archetypal two-gap superconductor MgB 2 8,9,15 . There, two characteristic distances in the vortex pattern have been reported, with two peaks in the histogram of nearest neighbor vortex distributions. One which is called an intergroup distance and follows a 0 when varying the magnetic field and an intragroup distance that remains almost constant with respect to the applied magnetic field. It has been argued that the vortex stripes are independent of the crystal lattice and therefore cannot be related to strong pinning. Authors of Refs. 6, 8, 9, and 15 relate instead their findings to the magnetic competition of two coexisting gap-condensates in superconducting MgB 2 . The vortex patterns we report here for β-Bi 2 Pd are different. At very low fields, the patterns do contain vortex stripes, clusters and vortex free regions. However, the intervortex distances do not cluster around two values as seen in MgB 2 . Instead, vortices arrange in lines along crystalline defects. As convinc-ingly shown in Refs. 22, 23, and 81, β-Bi 2 Pd is probably a single-gap superconductor, so vortex clustering cannot be associated to the hybridization of multiple gaps.
Vortex clustering has been seen in single-gap superconductors as well, with κ 1 2 . The intermediate mixed state consists of clusters of vortices with widely differing intervortex distances that are often smaller than d Abrikosov . Vortex lattices at fields close to or below H c1 imaged using magnetic decoration, scanning Hallprobe microscopy, and/or studied using small-angle neutron scattering often show a mixed intermediate state characterized by strong magnetic field gradients inside the sample, due to flux-free areas coexisting with areas having vortices inside 14,[82][83][84][85][86][87] . In high quality single crystals, such as e.g. Nb with κ = 1.1 1 2 , flux-free regions coexist with domains of vortex lattice 14 , where the shape of domains resembles the intermediate state found in type I superconductors 7 . Small-angle neutron scattering finds the intervortex distance inside those domains exactly as expected at H c1 16,47 . In our experiments, the magnetic flux integrated over large areas provides a magnetic induction of the same order in presence of an ordered vortex lattice and in presence of strong vortex clustering. Vortices are nucleated in the sample during the field-cooled procedure, and remain pinned even at very low magnetic fields, in spite of the strong field gradients produced when cooling below H c1 (T ). Such a feature has never been previously reported, to our knowledge, at low magnetic fields and in presence of strongly inhomogeneous vortex distributions.

APPENDIX D: FLUCTUATIONS VS PINNING
The elastic displacement of vortices u(0) caused by a pinning force F depends on the magnetic field and temperature following u(0) ≈ F ( 4π BΦ0 ) 1/2 µ 0 ) 3/2 , for an isotropic superconductor with penetration depth λ (µ 0 is the magnetic permeability) 2,85,88 . Close to steps, the increased λ thus yields an increased displacement of vortices u(0). In addition, B c1 ≈ Φ0 4πλ 2 (lnκ + 0.5) is decreased, whereas B c2 = Φ0 πξ 2 remains similar, provided that ξ is limited by the mean free path in the dirty limit regime.
Vortices suffer positional changes due to thermal activation, which can be discussed using elastic theory of the vortex lattice 2,85,89 .
The changes in the position are related to the temperature following < u 2 >≈ d 2 0 ( 3GiB π 2 Bc2(0) ) 1/2 T Tc λ 2 (T ) )) −1/2 (a 0 is the Abrikosov intervortex distance and G i the Levanyuk-Ginzburg number) 2,85,90-93 . The fluctuation range is very small, because in any event G i likely small in β−Bi 2 Pd. Anisotropy enhances fluctuations but there is practically no anisotropy in β − Bi 2 P d 22,23 . However, fluctuations are enhanced for large wavelengths, by taking into account nonlocal elastic interactions, and in presence of pinning induced disorder due to the mobility of dislocations and defects of the flux lattice 89,[94][95][96][97] . Vortex nucleation occurs here very close to the zero field T c for very low magnetic fields, which would suggest that fluctuations might play a role. Vortex entry occurs through thermal diffusion for long wavelengths k 2,85 . There are two processes of vortex diffusion, giving an exponential decay exp(−Γ 1 t) with rates Γ 1 and Γ 2 85,98 . First, through the lattice compression, giving exp(−Γ 1 t) with Γ 1 ≈ B µ0Bc2 ρ N k 2 with ρ N being the normal state resistivity. For very low magnetic fields of a few Gauss, and taking ρ N ≈ 10 −6 Ωm, we find Γ 1 ≈ 310 2 L 2 in s, with L the wavelength in meters. Even for relatively long wavelengths of order of the sample size, this process is fast, in the µs range or below. Second, through shear stress, with a rate Γ 2 ≈ c 66 k 2 . This process is much slower, because the shear modulus c 66 is small. If we take c 66 ≈ BΦ0 16πλ 2 µ0 we can find c 66 ∝ 1 L 2 and time scales three order of magnitude above estimation, but still of the order of a small fraction of a s. A typical cooling procedure after applying the magnetic field is certainly much slower, of the order of one s.
Therefore, we conclude that we do not identify a strong effect of fluctuations and that we are instead observing the sole effect of pinning in our experiments.

APPENDIX E: TEMPERATURE DEPENDENCE OF THE MFM IMAGES
The evolution of the flux line lattice with increasing temperature is shown in Fig. 8. In the MFM experiment, we observe the hexagonal vortex lattice up to temperatures of 4.5 K at a magnetic fields of 300 G. We mostly observe defect free vortex lattices and see no particular increase in the defect structure of the vortex lattice when increasing temperature. In Fig. 8(a) there are a few defects in the vortex lattice (vortices showing seven or five nearest neighbors). These are washed out when increasing temperature (Fig. 8(b),(c)). Thus, the range of vortex liquid is very small in β−Bi 2 Pd.
The diagonal line located at the bottom of the image and visible at all temperatures is probably caused by a crosstalk between charging effects between tip and sample and the magnetic signal close to a large step. The crosstalk between the two signals should not be temperature dependent and therefore we expect variations in this feature to be related to temperature induced variations in the magnetic properties. Part of this feature might thus be associated to a decrease in T c along this line. The crosstalk makes it however very difficult to make any quantitative estimation of T c . On the other hand, the image taken at the highest temperatures ( Fig. 8(b),(c)) show a whitish vertical line highlighted by a black arrow. That line is absent in other images as shown and might result from a change in T c along this line. With our present resolution in MFM, SOT and available STM experiments 22,23 we cannot be more quantitative. SOT suffers from problems with temperature control and the estimated change in T c corresponds to a change in the gap size of 50 µV which is impossible to detect in STM measurements taking into account the small but finite gap distribution close to 50 µV found in this compound 22,62 .

APPENDIX F: DETAILS ON THE CALCULATION OF THE FRACTAL PROPERTIES
To characterize a multifractal system, we mainly use two functions f (α) vs α and D q vs q.
f (α) vs α is the multifractal singularity spectrum and is typically a concave downward curve with a maximum and a certain width. The α value corresponding to the maximum in f (α) gives the fractal dimension. For example, if we consider pixels consisting of zeroes and ones distributed over a square, we can find different results, depending on the distribution of these pixels. A random distribution of pixels provides a small concave downward curve which is only defined very close to α = 2 and where f (α) is nearly constant. A monofractal distribution of pixels gives a maximum of f (α) at the fractal dimension. A multifractal distribution of pixels gives a concave downward f (α) defined over a wide range of α. D q vs q gives the generalized dimension for the set of scaling exponents q. The scaling exponents are used to highlight different regions of the image with more or less concentration of pixels. In random and in monofractal images, D q is a flat line because the dimension does not depend on q. In a multifractal, the dimension changes with the area, and thus with the scaling exponent q. We obtain a sigmoidal curve for D q vs q. q can be varied in the range of [∞, −∞] but for the implementation in the calculation, the limits depend on the convergence of the curve D q vs q for different ranges of q. In our case a convergence in D q was obtained in the range of [10, −10].
To go into the details of how to obtain f (α) and D q , we first remind a few aspects of the calculation we have made. We used the box counting method 99 . For a given binary matrix of points, as Figure 9(a), we calculate the number of points, m i ( ), in each box of length , and compute the probability of finding a white pixel in each box with: being N i( ) the number of boxes with length containing at least one point. Now we introduce the set of exponents q, which provide the dimensions in the multifractal spectra. We calculate: where I q represents how the pixels are distributed in space. The smaller I q , the larger the homogeneity in the number of pixels inside the boxes of same . µ qi( ) is equivalent to P i ( ) but taking into account the different behaviour of scaling exponents q. Then we calculate which gives the α q value by calculating the slope of log(A q ) vs log( ) as shown in Fig. 9(b). When a discrete set of points is given, as in the case of vortex image, only the box with a size of larger than minimum distance between first neighbours is taken into account. That is why the curves of Fig.9(b-c) flatten at small . We also compute shown in Fig. 9(c). From the slope of log(τ q ) vs log( ) we obtain τ q . Finally we can calculate To calculate D q and f (α) for the vortex lattice images (Fig. 1 ), we have first searched all vortex positions and The change of slope is due to the density of points. We only take in account the points before the change of slope. This gives us αq. In (c) we show log(τ q ) vs log( ) curves at the same set of q as in (b). The same behaviour with the size arises and we treat it similarly. The slope of this curve gives us τq.
calculated the images with one at a vortex and zero elsewhere. This leads to maps of pixels with values zero and one. The results are shown in Fig. 4. Images showing triangular or disordered lattices provide a distribution of fractal subsets centred at α = D q = 2. When we start observing variations in the vortex density, multifractality increases, which leads to f (α) that is much broader and whose maximum deviates from 2. D q also increases for small values of the multi-scaling exponent q. Thus, the vortex distributions at small magnetic fields are multifractal, with a probability of fractal subsets that strongly increases when decreasing the magnetic field, leading to a widening of f (α) and a considerable dependence of D q on q.

APPENDIX G: VORONOI TESSELATION OF THE IMAGES
Using the Delaunay triangulation, we have also calculated the Voronoi pattern of the vortex lattice images. Each vortex is then located inside a cell. We calculate the number of sides of each Voronoi cell (in the hexagonal lattice, each vortex cell has 6 sides) and the area, compared to the result in a perfect hexagonal lattice (a 2 0 ). The Voronoi patterns and the evolution of number of sides and area is shown in Fig. 10. At small magnetic fields, the number of sides of the cells significantly deviates from 6 ( Fig. 10(a),(c). More cells have less than 6 sides, although the distribution flattens out considerably. Cells with less than 6 sides corresponds to vortices aligned along a row, where vortices are located in cells with mostly 4 sides. On the other hand, vortices in between rows have six sides, or more when they do not have many neighbors. As can be expected (von Neumann's law), the area occupied by vortex cells increases with the number of sides (Fig. 10(d) for all magnetic fields). However, for low magnetic fields, the cells with more than 6 sides also have a larger area. This shows that vortices tend to be isolated in between rows and lines. Measurements of angular resolved photoemission in the normal phase above T c suggest non-trivial topological behavior at the surface of β-Bi 2 Pd 61,62 . The bandstructure observed by photoemission coincides with calculations and has mixed contributions from Pd 4d and Bi 6p orbitals that give three main sheets of convolved geometries that partially overlap with each other 71,100 . Photoemission reveals a Dirac cone well below the Fermi level 61 . Spin resolved measurements provide polarized bands close to the Dirac cone. The same authors suggest that topologically non-trivial spin polarized bands crossing the Fermi level might rise up to the surface. The STM experiment of Ref. 62 provides indications for a small triplet component appearing at the surface. In another STM experiment on epitaxially grown thin films, authors found superconducting properties that are different from the bulk behavior 17 -the critical temperature was somewhat larger and two gaps were detected in the tunneling conductance. Furthermore, a zero-bias peak appears in the center of the vortex cores, indicating the formation of vortex bound states [63][64][65] . Authors argue that these states could be topologically non-trivial, contrasting earlier results found in other superconducting materials. A recent report shows non-integer flux quantization in rings made of thin films of Bi 2 Pd, suggesting the formation of a π junction that authors associate with unconventional superconducting properties 19 .
Our results provide new insight in this debate. The flux carried by single vortices can vary with respect to the flux quantum Φ 0 in topological superconductors, due to enhanced stability of fractional vortices 101 . The observation of vortices with a weaker magnetic field profile as the ones located on a defect in our case could be explained by their flux below Φ 0 , instead of a varying penetration depth. However, the change in λ describes our results significantly better. Let us also remark that fractional flux quantization can also occur in long Josephson 0 − π junctions, where the phase difference varies along the junction 102,103 . Such junctions can form in superconductors with anisotropic gap structures or ferromagnetic inclusions. It is not clear how such a situation can be formed close to a linear defect in β-Bi 2 Pd and lead to the observed increase in λ.
Thus, there is mounting evidence for the presence of triplet correlations in β−Bi 2 Pd, particularly from experiments that are probing the surface. The question then arises of why we did not observe a half-integer flux quantum in our experiment. The answer might be that triplet correlations could form only close to the surface. Magnetic vortices, however, are threads through the bulk of the material and the surface properties might be masked by the flux quantization from the bulk. It would be very interesting to repeat our experiments in thin films or very thin layers of β−Bi 2 Pd. The appearance of half-integer flux quanta is, in view of the recent reports 19,61,62 , quite likely. Vortices might then carry a Majorana fermion, as proposed in Refs. 63-67. In that case, the large intervortex distances found in the range of magnetic fields we study here should be very helpful to facilitate manipulation of vortices and braiding experiments. For example, a couple of vortices located in between lines could be easily moved around each other.