# Unveiling the three-dimensional spin texture of skyrmion tubes

Daniel Wolf  Leibniz Institute for Solid State and Materials Research, IFW
Dresden, Helmholtzstr. 20, 01069 Dresden, Germany  Sebastian Schneider
Dresden Center for Nanoanalysis, cfaed, Technische Universität Dresden, 01069
Dresden, Germany  Leibniz Institute for Solid State and Materials Research,
IFW Dresden, Helmholtzstr. 20, 01069 Dresden, Germany  Ulrich K. Rößler
Leibniz Institute for Solid State and Materials Research, IFW Dresden,
Helmholtzstr. 20, 01069 Dresden, Germany  András Kovács  Ernst Ruska-Centre
for Microscopy and Spectroscopy with Electrons and Peter Grünberg Institute,
Forschungszentrum Jülich, 52425 Jülich, Germany  Marcus Schmidt  Department
Chemical Metal Science, Max Planck Institute for Chemical Physics of Solids,
Nöthnitzer Str. 40, 01187 Dresden, Germany  Rafal E. Dunin-Borkowski  Ernst
Ruska-Centre for Microscopy and Spectroscopy with Electrons and Peter Grünberg
Institute, Forschungszentrum Jülich, 52425 Jülich, Germany  Bernd Büchner
Leibniz Institute for Solid State and Materials Research, IFW Dresden,
Helmholtzstr. 20, 01069 Dresden, Germany  Institute of Solid State and
Materials Physics, TU Dresden, 01069 Dresden, Germany  Bernd Rellinghaus
Dresden Center for Nanoanalysis, cfaed, Technische Universität Dresden, 01069
Dresden, Germany  Axel Lubk  Leibniz Institute for Solid State and Materials
Research, IFW Dresden, Helmholtzstr. 20, 01069 Dresden, Germany  Institute of
Solid State and Materials Physics, TU Dresden, 01069 Dresden, Germany
Corresponding author.

Magnetic skyrmions 1 are stable topological solitons with complex non-coplanar
spin structures. Their nanoscopic size and the low electric currents required
to initiate and control their motion has opened a new field of research,
skyrmionics, that aims at using skyrmions as information carriers for data
storage and manipulation 2, 3, 4, 5. Recent advances in skyrmionics prompt for
a thorough understanding of the detailed three-dimensional spin texture of a
skyrmion 6, 7, 8, 9, 10, 11 including skyrmion-skyrmion interactions and their
coupling to surfaces and interfaces. These properties crucially affect
application-related aspects such as the stability and mobility of skyrmions in
confined structures. To date, however, experimental techniques to measure the
three-dimensional (3D) spin texture with nanometer resolution are largely
missing. We therefore adapt holographic vector field electron tomography12 to
the problem and report on the first quantitative reconstruction of the 3D spin
texture of skyrmions with sub-10 nanometer resolution. The reconstructed
textures reveal a variety of previously unseen local deviations from a
homogeneous Bloch character within the skyrmion tubes (SkTs), details of the
collapse of the skyrmion texture at surfaces, and a correlated modulation of
the SkT in FeGe along their tube axes. The quantitative 3D data of the
magnetic induction also allow to experimentally confirm some principles of
skyrmion formation by deriving spatially resolved maps of the magnetic energy
density across these magnetic solitons.

## Introduction

The unique features of magnetic skyrmions such as their competing magnetic
interactions, topological structure, or dynamics are of great interest in both
fundamental and applied physics. As multidimensional solitons, these particle-
like states are localized in two-dimensions which requires a definite
mechanism through additional frustrating magnetic couplings for their
stabilization 13. As consequence of their solitonic character, they can
condense into thermodynamically stable phases, in particular dense packed
lattices under applied fields 1. The stabilization mechanism of these phases
and their formation principles are ruled by effective skyrmion-skyrmion
interactions 13. However, the morphology of these phases in the phase-diagrams
of real materials is dictated by the problem of condensation of two-
dimensional periodic arrays, as in vortex-lattices of type-II superconductors
1. In particular, the field-temperature phase diagram may hold various
transitions between different condensed phases of skyrmions 14, 15. Very
recently, some studies address this problem for skyrmionic phases
theoretically 16 and in experiments 17. In three dimensional bulk materials or
thicker films the skyrmions are extended string-like objects; in the simplest
formation they are homogeneously continued as skyrmion tubes (SkT) preserving
translational invariance along their axis. In magnetic nanoobjects, however,
the influence of surfaces will affect formation, shape and interaction of
skyrmions and the stabilization of condensed skyrmionic phases. Already from
the earliest observations of skyrmionic phases in films of chiral helimagnets
18, 19, it is known that their phase diagram massively deviate form those of
bulk materials. It has been shown that 3D surface twists can stabilize SkTs in
thin films 6, 7, 20, 21 and that 3D modulations of SkTs embedded in a conical
host phase may introduce an attractive interaction between these tubes 22. 3D
SkT modulations also affect emergent electric and magnetic fields acting on
spin-polarized electrons and magnons 23, which results in unusual transport
phenomena 24 on top of the ”normal” topological Hall effect in static and
current-driven skyrmion crystals 25, 26, 27, 28.

Figure 1: a, Holographic vector-field electron tomography (VFET) of skyrmions
in FeGe. A needle-shaped sample is placed above an out-of-plane magnetized
$\mathrm{Sm_{2}Co_{17}}$ ring in a liquid nitrogen cooled TEM holder. The low
temperature and the remanent stray field of the ring stabilizes skyrmion tubes
and their orientation with respect to the holder. A tilt series of 2D phase
images of the transmitted electron wave is obtained from off-axis electron
holographic tilt series around the $x$ (left) and $y$ axes (right).
Subsequently, $x$ and $y$ components of the magnetic induction
$\boldsymbol{B}$ are tomographically reconstructed from the corresponding
phase tilt series. Solving $\mathrm{div}\boldsymbol{B}=0$ finally yields the
$z$ component, hence the full 3D vector of $\boldsymbol{B}(x,y,z)$. b, 3D map
of the resulting magnetic induction. For clarity, only the experimental
$B_{x}$ and $B_{y}$ components are shown here.

Similarly, e.g., in hybrid chiral ferromagnet/superconductor systems 29, the
functionalization of skyrmions is modified by 3D modulations of the SkT at the
interface. Finally, the observation of unusually strong topological quantum
Hall effects 30 may indicate the presence of abrupt magnetization changes such
as in magnetic bobbers (incorporating a Bloch point) in surface regions 31.

Notwithstanding the importance of 3D effects, neither exact 3D models of SkTs
in realistic confined geometries nor high-resolution experimental mappings of
their spin texture are currently available, although effects of confinement
19, 6, 20 and of anisotropies 32 have received attention. This lack of data
prevents a deeper understanding of skyrmion lattice defects 11, 33, influence
of surface anisotropies, curvatures 34, and real structure effects in the
modulation of 3D skyrmionic spin textures. Among the various high-resolution
magnetic imaging techniques, transmission electron microscopy (TEM) based
electron holography (EH)35, 36, 37, 12 and X-ray magnetic chiral dichrosim
(XMCD) 38, 39, 40 can be conducted in a tomographic way to determine the 3D
magnetic induction, $\boldsymbol{B}$, or magnetization, $\boldsymbol{M}$, of a
sample, respectively. In this work, we employ holographic vector-field
electron tomography (VFET) as it provides a higher spatial resolution (below
10 nm) than X-ray based methods, which is crucial for resolving the details of
magnetic textures in nanomagnetic structures such as vortices 12 or skyrmions.
The limited space in a high-resolution TEM instrument, however, has so far
prevented any in-situ applications of rotatable (out-of-plane) magnetic fields
to a cryogenically cooled sample, which is essential for the acquisition of
tomographic tilt series of electron holograms from a sample that needs to be
magnetically stabilized. This limitation impedes the measurement and 3D
reconstruction of spin textures for a large class of materials with a
metastable skyrmion phase at non-zero applied fields below room temperature
(e.g., many isotropic helimagnets). For the present experiments, we have
therefore devised a setup that overcomes these obstacles.

## High resolution vector-field tomography in an external magnetic field

Figure 2: Spin texture of skyrmion tubes in a FeGe needle: a, Volume rendering
(colored) of the in-plane ($x$, $y$) components of the reconstructed magnetic
induction $\boldsymbol{B}$ and iso-surface of the mean inner potential
highlighting the sample shape (grey, bottom half only). Rectangles indicate
cross-sections, whose details are shown in c and d. b, Planar $x$-$y$ cross-
section. Color and size of the arrows indicate the direction and magnitude of
the in-plane component of $\boldsymbol{B}$, respectively. c and d, direction
of $\boldsymbol{B}$ (arrow orientation) and magnitude (color) of $B_{x}$ and
$B_{y}$ in $y$-$z$ and $x$-$z$ cross-sections through SkT 3. Here, the SkT was
artificially aligned along its $z$ axis (see text for details).

The tomographic investigation of the magnetic texture of skyrmions was
conducted on a sample of the isotropic helimagnet FeGe with $\mathrm{P2_{1}3}$
structure (B20 phase). The material was chosen, since FeGe is an otherwise
well studied archetypical skyrmion host with a rather large skyrmion phase
pocket in the phase diagram spanned by temperature and external field 18, 41.
A needle-shaped sample (cf. Supplementary Sect. LABEL:suppl:note:geometry) was
cut from a FeGe single crystal by focused ion beam (FIB) including ion
polishing to restrict the ion beam damage to a surface layer of some
nanometers (see Methods). The dimensions and shape of the needle ensure that,
even at high tilt angles, the sample is fully electron-transparent and the
obtained holographic projections cover the same sample region. Additionally,
the elongated shape has some technological significance for anticipated
spintronic devices such as racetrack memories 42. In order to (i) adjust the
skyrmion phase below the Curie temperature and (ii) stabilize the orientation
of the skyrmion lattice with respect to the TEM holder, the FeGe needle was
steadily exposed to an out-of-plane magnetic field of
$\mu_{0}H_{\mathrm{ext}}\approx 170\,\mathrm{mT}$. The field was provided
through the remanent stray field of a ring-shaped $\rm Sm_{2}Co_{17}$ hard
magnet that was placed under the sample in a tomography-adapted liquid
nitrogen TEM cooling holder (cf. Fig. 1a). The field is virtually homogeneous
across the micron-sized sample (cf. Supplementary Sect.
LABEL:suppl:note:ring).

Using this special setup, we have recorded three holographic tilt series as
required for VFET 36, 12 (see Methods and Supplementary Sect.
LABEL:suppl:note:geometry for details of the imaging conditions). The first
series of holograms was acquired by tilting the sample around the $x$ axis at
room temperature, since above the Curie temperature of $\rm T_{C}=278.7\,K$
43, the phase $\varphi_{e}$ reconstructed from the holograms is of pure
electrostatic origin. The (scalar) electrostatic potential $\rm\Phi$ was then
determined by inverting the Radon transformation (i.e., linear projection law)
linking $\rm\Phi$ and $\varphi_{e}$. The resulting 3D mean inner potential
distribution is nearly homogeneous as discussed in Supplementary Sect.
LABEL:suppl:note:MIP (see Methods for the tomographic reconstruction details).

In the following two series, the sample was tilted around the $x$ and $y$ axes
(cf. Fig. 1a) at $T=95\,\mathrm{K}$. At this temperature below $T_{C}$, the
magnetic fields imposes an additional Aharonov-Bohm phase $\varphi_{m}$ on the
imaging electrons. After subtracting the pre-determined electrostatic
contribution from the total phase shift, the in-plane components of the
magnetic induction, $B_{x}(x,y,z)$ and $B_{y}(x,y,z)$ (see Fig. 1b), were
reconstructed from the remaining $\varphi_{m}$ in 3D by inverse Radon
transformation of another linear projection law linking the gradient of
$\varphi_{m}$ and $B_{x,y}$ (see Methods for details). The spatial resolution
of the reconstructed $B_{x}$ and $B_{y}$ components was below
$10\,\mathrm{nm}$ in directions outside of the missing tilt range (cf.
Supplementary Sect. LABEL:suppl:note:resolution).

In order to change the tilt axis from $x$ to $y$, the sample required to be
warmed up to room temperature and rotated in-plane by $90^{\circ}$ in the
sample holder. The associated heating and (field) cooling of the sample
followed precisely the same protocol as prior to tilting it around $x$. As a
result, at the here investigated and most confined tip region of the FeGe
needle, the skyrmion patterns obtained after cooling prior to acquiring the
$x$ and $y$ tilt series were almost perfectly identical, while at the less
confined broader end of the needle, the skyrmion pattern had changed
considerably (see Suppl. Fig. LABEL:suppl:Fig_Specimen for details).

Based on $B_{x,y}$, also the remaining third component $B_{z}(x,y,z)$ was
finally determined by solving $\mathrm{div}\,\boldsymbol{B}=0$, thereby
yielding the full 3D vector field of the magnetic induction
$\boldsymbol{B}(x,y,z)$ (see Methods for details and Supplementary Movie 1 for
3D animations of the tomograms). Since the calculation of $B_{z}$ is, however,
based on the differentiation of $B_{x}$ and $B_{y}$, it suffers from noise
amplification and corresponding artefacts, which needs to be taken into
account in the following analysis.

## Spin texture of Skyrmion tubes

In the following, we analyze this comprehensive 3D set of
$\boldsymbol{B}(x,y,z)$ data in order to extract characteristic magnetic
features and quantities of the SkTs in FeGe. Fig. 1b reveals that the tip of
the needle hosts a single row of SkTs that are elliptically distorted towards
the sideward surfaces of the needle, i.e., perpendicular to both the needle
axis and the stabilizing external field. Towards the broader back of the
needle (top region of Fig. 1b), these elongated SkTs develop into a zig-zag
chain of Bloch SkTs, when the width surpasses a critical value of roughly 150
nm. This width corresponds to about twice the characteristic helical
modulation length $L_{D}$ and the next-nearest neighbour distance in a close-
packed skyrmion lattice in FeGe 33. An evaluation of the out-of-plane
component of $\boldsymbol{B}$ (cf. Supplementary Sect.
LABEL:suppl:note:Lattice_Packing) reveals a ratio of areas with positive and
negative $B_{z}$ of 0.85, which also points to a close-packing of the SkTs in
this region 16.

Figure 3: Axial modulations of the SkTs. a, 3D rendering of the ($x$, $y$)
components of $\boldsymbol{B}$ for eight close-packed SkTs in the FeGe needle.
The red lines represent the SkT axes, and $\boldsymbol{q_{1}}$,
$\boldsymbol{q_{2}}$, $\boldsymbol{q_{3}}$ indicate the three NN directions of
the SkT lattice. b and c, $z$ dependent positions of the cores of SkTs 1-5 and
SkTs 6-8 (f and g) along $\boldsymbol{q_{2}}$ and $\boldsymbol{q_{1}}$. The
data points are colored according to the labels in a. d and e, $x$-$z$ slices
through SkTs 3 d and 8 e.

Fig. 2 and Supplementary Movies 2, 3 represent an in-depth view of the spin
texture within these SkTs. Fig. 2b shows a planar cross-section of the spin
texture through the vertical center of the needle. Here, the color of the
arrows indicate the direction of the in-plane ($x$, $y$) components of
$\boldsymbol{B}$ according to the color wheel in Fig. 2a. While all four SkTs
in this section feature radial Bloch walls, details of the spin texture
exhibit subtle discrepancies from that of an undisturbed, perfect Bloch
skyrmion: (i) The skyrmions may exhibit significant distortions and
(partially) lose their axial symmetry. Neither the direction of the in-plane
component of $\boldsymbol{B}$ remain tangential nor is its magnitude constant
for a given radius. (ii) Frequently, the maximum out-of-plane orientation (as
indicated by a virtually vanishing size of the arrows) is not located in the
center of the skyrmions. (iii) Unlike expected for isolated magnetic solitons,
some distortions of the skyrmionic spin textures seem to go along with
magnetic flux ”leaking” between neighboring SkTs as highlighted by the dashed
region ”I” between SkTs 4 and 7. The mere similarity of this appearance with a
confined helical ”band” resembles the evolution of a metastable isolated
skyrmion towards a helical modulation (”strip-out”) discussed in bi-layer thin
films 44. Both scenarios, however, violate the topology of a skyrmion and
necessitate the occurrence of singular Bloch points. Features indicating the
occurrence of Bloch points are in fact frequently observed: Figs. 2c, d show
two orthogonal vertical slices through SkT 3. Since the SkT axes are found to
be bent and (partially) twisted (see below and Fig. 3), SkT 3 was artificially
aligned along the $z$ axis for this presentation. To this end, each $x$-$y$
slice of the tube was laterally shifted such that the minima of the in-plane
component $B_{\|}=\sqrt{B_{x}^{2}+B_{y}^{2}}$ of all slices are aligned along
the $z$ axis. The dashed circles labeled ”II” and ”III” highlight sections
that are indeed reminiscent of a Bloch point. Besides, both cross-sections
confirm the lack of axial symmetry and substantiate the overall inhomogeneity
of the spin texture in the SkT already seen from the planar cross-section. In
contrast to ”pure” Bloch SkT, we frequently (but not systematically) observe a
(partial) radial (i.e., a small Néel character) orientation of the local
magnetic induction. These imperfections grow upon approaching the surface and
finally lead to a total collapse of the skyrmion structure. This becomes most
apparent in the $x$-$z$ cross-section in Fig. 2d, where the thickness of the
needle decreases. This region should be understood as result of surface
symmetry breaking and concomitant effects, such as pinning by surface
anisotropies, modified magnetic properties due to FIB surface damage and
demagnetization fields.

Figure 4: Planar maps of the predominant magnetic energy density contributions
$DB_{z}(\nabla\times\boldsymbol{B})_{z}$ (DM interaction, (a and b)) and
$A\|\left(\nabla\times\boldsymbol{B}\right)_{z}\|^{2}$ (exchange interaction,
(c and d)) and their sum (e and f), together with radial averages as function
of the radius (g, h). Left: Simulations. Right: Reconstructed from the central
part of SkT 3. Overlaid arrows indicate the in-plane magnetic induction.

In the cubic helimagnet FeGe, a twisting in the third direction, i.e. along
the axis of the SkT, could result in a gain of energy through the
Dzyaloshinskii-Moryia (DM) exchange. However, such an effect will not create
triply twisted structures of skyrmions, as this is geometrically impossible,
because the ferromagnetic vector can be rotated only in two directions in the
cutting plane perpendicular to the skyrmion axis. Hence, the chiral twist6
could only affect the shape of the SkTs. E.g., a modulation could arise as a
tertiary conformational deviation from straight cylindrical SkT shape.
Confinement or surface pinning may promote such morphology changes. Indeed,
Fig. 3a illustrates that the axis of SkTs (red lines) are axially bent and
twisted rather than extending as straight cylindrical objects along the $z$
axis in the close-packed region of the needle (similar rendering of $B_{\|}$
as in Fig. 2a). In order to study possible correlations between these
deformations, we have analyzed the in-plane positions of the SkT axes along
the nearest neighbour (NN) directions $\boldsymbol{q_{1}}$,
$\boldsymbol{q_{2}}$, $\boldsymbol{q_{3}}$ (indicated by white arrows). The
resulting dependencies of the deviations from an average axial position along
$\boldsymbol{q_{2}}$ and $\boldsymbol{q_{1}}$, i.e., in directions that are
largely affected by the lateral confinement, are shown in Fig. 3b, c for the
bottom row of SkTs (nos. 1-5) and in Fig. 3f, g for SkTs 6-8 in the top row.
Except for SkT 1 (small blue circles in b and c), which is least close-packed
and rather has two elliptically elongated SkT neighbours, and SkT 7 (small
pink circles in b and c), which is additionally distorted due to an unusual
magnetic coupling to SkT 4 (see above), all SkT axes exhibit pronounced
sideward deformations. As indicated by grey bands (guides to the eye only),
these lateral modulations are correlated among SkTs in the same row. These
deformations are harmonically modulated with a modulation length of
approximately $\rm 80\,nm$ that is close to the helical period $L_{D}\simeq
70\,\mathrm{nm}$ in FeGe 18, 45 pointing to the DM interaction as a possible
origin of the deformations. Note, however, that comparisons with $y$-$z$
cross-sections through SkTs 3 and 8 in Figs. 3d, e reveal that these
modulations correlate with the occurrence of uniformly polarized edge states
46. These edge states reside at the sidelong rims of the FeGe needle (cf. left
and right surfaces in Figs. 1b and 2a) and are separated by a very narrow
magnetic transition regions (resembling domain walls) of some 10 nm in width
from the SkTs. The correlation of the deformation of the SkTs with these edge
states is corroborated by the facts that (i) the central deformations are
directed towards the center on both sides of the needle and (ii) the magnetic
orientations of the edge states and the outer rims of the SkTs’ spin textures
are concurrently reversed between the right (SkTs 1-5) and left side (SkTs
6-8) of the needle, respectively. This results in qualitatively identical
interactions between the SkTs and the edge states on either side. In contrast,
the deformations of the SkTs along the largely unconfined direction
$\boldsymbol{q_{3}}$ do not exhibit any obvious correlations (not shown).

The availability of 3D vector data of the magnetic induction enables us for
the first time to experimentally derive from the volume of a sample spatial
maps of free energy density contributions from magnetic exchange and DM
interactions, respectively. These energetic contributions are most essential
for the formation and stabilization of skyrmions and SkTs, as they are
expected to reduce the free energy in the centers of the SkTs, while the
interstitial regions in a SkT lattice may be considered as domain walls of
increased energy 47. We have calculated from $\boldsymbol{B}(x,y,z)$ the
solenoidal part of the magnetic exchange energy density
$\mathit{w}_{\mathrm{ex}}=A\|(\nabla\times\boldsymbol{B})\|^{2}$ and the
volume contribution of the DM energy density
$\mathit{w}_{\mathrm{DM}}=D\boldsymbol{B}\cdot(\nabla\times\boldsymbol{B})$.
Here, $A=8.78\ \mathrm{\frac{pJ}{m}}$ and $D=1.58\ \mathrm{\frac{mJ}{m^{2}}}$
denote the exchange stiffness and the DM interaction strength of FeGe,
respectively 48. Due to the vanishing magnetic charge density $\rho_{m}\approx
0$, the conservative part of the exchange energy
$\|\nabla\cdot\boldsymbol{M}\|^{2}$ is small in Bloch skyrmions, and other
contributions are divergence contributions that can be collapsed to surface
terms, which is the reason why both are neglected here. As the spin texture of
the SkTs is highly disturbed in the near-surface region (cf. Fig. 2c, d), and
in order to account for the axial deformation of the SkTs, only magnetic
induction data from the central part of the SkT (cf. grey shaded boxes in
Figs. 2c, d) was used and projected in the $x$-$y$ plane to calculate the
planar distribution of energy densities. For comparison, such energy density
maps were also calculated for a simplified skyrmion lattice model employing
the circular cell approximation 49 taking into account the shape of the needle
(Supplementary Sect. LABEL:suppl:note:Magnetostatics).

Fig. 4 shows the resulting simulated (left column) and experimentally
determined energy density maps (right column) for the contributions arising
from the DM and exchange interactions, and their sum. Here, we only plot
energy densities dominated by the in-plane components of the magnetic
induction to suppress some of the artefacts afflicting the $B_{z}$ component,
which is, however, sufficient and consistent with calculations that take the
full $\boldsymbol{B}(x,y,z)$ into account (cf. Supplementary Sect.
LABEL:suppl:note:mag-energies and Fig. LABEL:Suppl_Fig_energetics_simulation).
Given the apparent real structure effects and noise in the data, the
experimental results are in excellent qualitative agreement with the
simulations. In particular, the course of the radially averaged contribution
(g, h) confirms experimentally the prediction that the reduction in the free
energy density due to the DM interaction overcompensates the energetic costs
of the exchange in the core of the skyrmion tube, which results in the overall
energetic stabilization of the SkT (lattice). In comparison with the
simulations, the experimental energy density landscapes are slightly
compressed in $x$ direction, which is attributed to the interaction with the
edge state and contributes, besides the noise, to the quantitative reduction
of the radial averages in Fig. 4h. Noteworthy, the energetically least stable
part of the SkT is the center of the circumferential Bloch wall, which is the
region with the highest in-plane orientation of the magnetic induction in the
SkT. It is remarkable that this is precisely that region, where we had found
indications of ”leaking flux” between SkTs 4 and 7 (see discussion above and
Fig. 2b). Apparently, the center of the Bloch wall is the most ”forgiving”
zone of the SkT, which in turn may explain its overall stability against the
variety of observed magnetic defects. All in all the energy maps confirm the
heterogeneous particle-like nature of skyrmions. Unlike the energetically
homogeneous spiral states of a chiral helimagnet like FeGe, the skyrmions have
a definite shape and size that is caused by the frustration between the
different exchange energies, which can be lifted only partially through the
doubly twisted core.

## Summary

Low temperature holographic vector-field electron tomography was used in
combination with the spatial stabilization of the specimen’s magnetic state by
an external magnetic field to reconstruct the full vector-field
$\boldsymbol{B}$ of the skyrmionic spin texture in FeGe in all three
dimensions at nanometer resolution. The unrivaled resolution of this 3D
magnetic microscopy of a volume sample reveals unprecedented insight into the
details of the 3D spin texture of skyrmions.

Besides a characterization of the complicated breakdown of the skyrmion
texture upon approaching surfaces in axial directions, we observe a variety of
imperfections in the spatial extension of skyrmion tubes. Among them are axial
and planar distortions of the SkTs, local losses of axial symmetry, and the
occurrence of unexpected radial (rather than purely tangential) tilts of the
magnetic induction in the circumferential Bloch walls. Even indications of in-
plane magnetic flux leaking among neighboring SkTs in close-packed regions and
abrupt changes of the magnetic induction that may be indicative of the
occurrence of Bloch points are found. Also, the 3D course of the SkT axes was
investigated in great detail. Here, we observe a substantial bending and
twisting of these axes that is locally correlated with the occurrence of
pronounced edge states, specifically in directions that are affected by
confinements. Noteworthy, these deformations appear at length scales, where
harmonic modulations are promoted by the DM interaction.

Planar energy density maps across the SkTs were derived from the volume data
of the magnetic induction and confirm for the first time experimentally the
anticipated formation and stabilization mechanism of skyrmions for a volume
sample. The results reveal a substantial energetic gain due to the DM
interaction that overcompensates the energetic effort associated with the
magnetic exchange interaction in the core of the SkT thereby stabilizing the
SkT lattice as a whole.

We anticipate that this novel experimental approach will pave the way to a
better understanding of spin textures in a large variety of complex
topologically protected and non-topological magnetization patterns, including
other members of the skyrmion family, thereby moving the fields of both
nanomagnetism and spintronics significantly forward.

## Methods

### Sample preparation.

Based on the results of crystal growth by chemical vapour transport in the
system Fe/Ge 50, 51 single crystals of FeGe in the B20 structure were grown
via chemical transport reaction using iodine as transport agent. Starting from
a homogeneous mixture of the element powders iron (Alfa Aesar 99,995 %) and
germanium (Alfa Aesar 99,999 %) the cubic modification of FeGe crystallized by
a chemical transport reaction very slowly in a temperature gradient from 850 K
(source) to 810 K (sink), and a transport agent concentration of
$0.2\,\mathrm{\frac{mg}{cm^{3}}}$ iodine (Alfa Aesar 99,998 %). The chemical
vapour transport was made perpendicular to the tube axis over a diffusion
distance of 38 mm. Selected crystals were characterized by EDXS, WDXS, and
especially X-ray single crystal diffraction to verify the present
modification.

The preparation of the FeGe needle was carried out via focused ion beam (FIB)
technique on a Thermo Scientific Helios 660. A rough cut of the needle
geometry $(700\times 700\,\mathrm{nm})$ was performed with currents of 790 and
$430\,\mathrm{pA}$ respectively. For further fine shaping $(300\times
300\,\mathrm{nm})$ the current was reduced to 80 and $40\,\mathrm{pA}$. The
final polishing was carried out at $24\,\mathrm{pA}$. In order to remove
preparation residue the needle was finally cleaned in a Fischione Model 1070
NanoClean for $1\,\mathrm{min}$.

High-resolution TEM images indicate, that the crystalline core of the needle
is surrounded by an amorphous surface layer of 4 nm. STEM-EDX measurements
reveal, that this layer consists of iron oxide (cf. Supplementary Sect.
LABEL:suppl:note:surface_layer).

The skyrmion phase within the needle was stabilized by the stray field of a
ring-shaped $\mathrm{Sm_{2}Co_{17}}$ permanent magnet fitted into a GATAN 636
double tilt liquid nitrogen sample holder. The ring was prepared from a bulk
magnet by sinker spark eroding and mechanical grinding and provided a magnetic
field in $z$-direction of approximately $170\,\mathrm{mT}$ (cf. Supplementary
Sect. LABEL:suppl:note:ring).

### Acquisition and reconstruction of the holographic tilt series.

Holographic tilt series were recorded at an FEI Titan G2 60-300 HOLO52 in
Lorentz-Mode (conventional objective lens switched off) operated at
$300\,\mathrm{kV}$. The voltage of the electrostatic Möllenstedt bisprism was
set to 120 V leading to a fringe spacing of $2.3\,\mathrm{nm}$ in the electron
hologram (cf. Supplementary Sect. LABEL:suppl:note:solitonic_skyrmions). For
the acquisition of the latter, a GATAN K2 Summit direct detection camera in
counting modes was used yielding a holographic fringe contrast of $40\,\%$.
The acquisition process was performed semi-automatically with an in-house
developed software package 53 to collect three holographic tilt series
consisting of object and object-free empty holograms, two at $95\,\mathrm{K}$
and one at room temperature. For the first tilt series at $95\,\mathrm{K}$,
the angle between the needle and tilt axis amounts to $30^{\circ}$. For the
second tilt series, the specimen was manually rotated outside the microscope
in-plane by $70^{\circ}$ (ideal is $90^{\circ}$) resulting in an angle between
the needle and tilt axis of $-40^{\circ}$ (cf. Supplementary Sect.
LABEL:suppl:note:geometry for the details). The tilt range of each tilt series
was from $-66^{\circ}$ to $+65^{\circ}$ in $3^{\circ}$ steps. To obtain the
full phase shift ($>2\pi$), the phase images were unwrapped automatically by
the Flynn algorithm 54 and manually at regions, where the phase signal is too
noisy or undersampled, by using prior knowledge of the phase shift (e.g., from
adjacent projections) 55. Potential phase wedges in vacuum caused by the
magnetic stray-field of the ring were corrected in all three tilt series. An
analysis of these stray-field contributions is presented in the Supplementary
Sect. LABEL:suppl:note:el_mag_fringing_fields.

### Tomographic reconstruction.

All three phase image tilt series were aligned, i.e., corrected for image
displacements with respect to their common tilt axis by cross-correlation,
center-of-mass method and common-line approach12. Thereby obtained aligned
datasets correspond to the following linear projection laws (Radon
transformations):

$\varphi_{e}\left(p,\theta,z\right)=C_{E}\iint_{\boldsymbol{e}\cdot\boldsymbol{r}}\Phi\left(x,y,z\right)\mathrm{d}x\mathrm{d}y$
(1)

and

$\frac{\partial\varphi_{m}(p,\theta,z)}{\partial
z}=\frac{e}{\hbar}\iint_{\boldsymbol{e}\cdot\boldsymbol{r}}B_{p=y,x}\left(x,y,z\right)\mathrm{d}x\mathrm{d}y.$
(2)

Here, $C_{E}$ is a kinetic constant depending solely on the acceleration
voltage, $p$ and $z$ are the 2D detector coordinates, $\theta$ the tilt angle,
$\boldsymbol{r}=(x,y)^{T}$, and
$\boldsymbol{e}=(\cos{\theta},\sin{\theta})^{T}$. The index to the integral
indicates a collapse of the 2D integral to the projection line defined by
$\boldsymbol{e}\cdot\boldsymbol{r}$. The subsequent tomographic 3D
reconstruction of the aligned phase tilt series (i.e., the inverse Radon
transformation) was numerically carried out using weighted simultaneous
iterative reconstruction technique (W-SIRT) 56.

The three resulting tomograms represent the incremental 3D phase shift per
voxel that we refer to as 3D phase maps. The two 3D phase maps obtained at
$95\,\mathrm{K}$ were released from their electrostatic (MIP) contribution by
superposition and subtraction of the 3D phase map obtained at room
temperature. Then, the derivation of each of the two resulting magnetic 3D
phase maps in direction perpendicular to the experimental tilt axis and
multiplication with the factor $\hbar/e$ leads to one component of the
magnetic induction in the respective direction. Since the specimen was rotated
only by $70^{\circ}$ in the underlying tomographic experiment for the
reconstruction of these two $\boldsymbol{B}$-field components, one of them was
projected on the orthogonal direction of the other to receive finally the 3D
$B_{x}$ and $B_{y}$ component. A verification of the experimental workflow
repeated on simulated data is provided in Supplementary Sect.
LABEL:suppl:note:Rec-of-Sim.

### Calculation of the third magnetic B field component.

The third $\boldsymbol{B}$ field component $B_{z}$ is obtained by solving
Gauss’s law for magnetism $\mathrm{div}\,\boldsymbol{B}=0$ with appropriate
boundary conditions on the surface of the reconstruction volume. Here, we
employed periodic boundary conditions for solving this differential equation
in Fourier space endowed with coordinates $\boldsymbol{k}$, i.e.,

$B_{z}\left(\mathbf{k}\right)=-\frac{k_{x}B_{x}\left(\mathbf{k}\right)+k_{y}B_{y}\left(\mathbf{k}\right)}{k_{z}}$
(3)

The zero frequency component (integration constant) was fixed by setting the
average of $B_{z}$ to zero on the boundary of the reconstruction volume. To
suppress noise amplification by this procedure, a butterworth-type low pass
filter was applied.

### Magnetic energy densities.

Following Ref. 12 the exchange energy density may be split into contributions
from magnetic charges, currents and surface terms
$A\left(\left(\nabla\cdot\boldsymbol{M}\right)^{2}+\left|\nabla\times\boldsymbol{M}\right|^{2}-\mathit{w}_{\mathrm{surf}}\right)$.
In the magnetostatic limit considered here, the magnetization in the second
term may be replaced by $\boldsymbol{B}/\mu_{0}$ and can be reconstructed from
the tomographic data. In the case of the DM interaction we have the following
identities

$\displaystyle E_{\text{DM}}\left[\boldsymbol{M}\right]$
$\displaystyle=D\int\boldsymbol{M}\cdot\left(\nabla\times\boldsymbol{M}\right)dV$
(4)
$\displaystyle=\frac{D}{\mu_{0}}\int\boldsymbol{B}\cdot\left(\nabla\times\boldsymbol{B}\right)dV+D\int\nabla\cdot\left(\Phi\boldsymbol{j}_{b}\right)dV$
$\displaystyle=\frac{D}{\mu_{0}}\int\boldsymbol{B}\cdot\left(\nabla\times\boldsymbol{B}\right)dV+D\varoiint\boldsymbol{w}_{\mathrm{surf}}\cdot
d\boldsymbol{S}$

Here $\boldsymbol{j}_{b}$ denotes the bound current and $\Phi$ the scalar
magnetic potential. The last line identifies that part of the DM energy
density, which can be derived solely from the $\boldsymbol{B}$-field, may be
identified as a volume contribution, which can be reconstructed from
tomographic data. The remainder can be collapsed to a surface term.

#### Acknowledgements

We thank D. Pohl for helpful discussions in the process of planning the
experimental setup and the data analysis. We furthermore acknowledge A. Tahn
and T. Walter for the preparation of the FIB needle and magnetic ring,
respectively. The authors are indebted to Vacuumschmelze GmbH & Co. KG for
providing the $\mathrm{Sm_{2}Co_{17}}$ magnet. AL, BR, and SS gratefully
acknowledge financial support through the Priority Program SPP2137 of the
German Research Foundation (DFG) within projects LU-2261/2-1 and RE-1164/6-1.
DW and AL have received funding from the European Research Council (ERC) under
the Horizon 2020 research and innovation program of the European Union (grant
agreement number 715620). AK and RD received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant No. 856538, project “3D MAGiC”), from the European
Union’s Horizon 2020 Research and Innovation Programme (Grant No. 823717,
project “ESTEEM3”) and from the European Union’s Horizon 2020 Research and
Innovation Programme (Grant No. 766970, project “Q-SORT”).

### Author Contributions

SS devised the experimental setup for the magnetic field stabilization. DW
conducted the holographic VFET experiments with active support of SS and AK
and performed the holographic and tomographic reconstructions. SS, DW, BR and
AL analyzed the data. AL performed the magnetic simulations. AL, SS, DW, BR
and UR wrote the manuscript. The FeGe single crystal was grown by MS. All
authors contributed to the critical discussion and revision of the manuscript.

### Supplementary Information

Supplementary Information is available for this paper.

### Author Information

Correspondence and requests for materials should be addressed to AL
(a.lubk@ifw-dresden.de).

## References

  * 1 Bogdanov, A. & Yablonskii, D.  Thermodynamically stable ”vortices” in magnetically ordered crystals. The mixed state of magnets.  _Zh. Eksp. Teor. Fiz_ 95, 178 (1989).
  * 2 Kiselev, N. S., Bogdanov, A. N., Schäfer, R. & Röler, U. K.  Chiral skyrmions in thin magnetic films: New objects for magnetic storage technologies?  _Journal of Physics D: Applied Physics_ 44, 392001 (2011).
  * 3 Sampaio, J., Cros, V., Rohart, S., Thiaville, A. & Fert, A.  Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures.  _Nature Nanotechnology_ 8, 839–844 (2013).
  * 4 Tomasello, R. _et al._ A strategy for the design of skyrmion racetrack memories.  _Scientific Reports_ 4, 1–7 (2014).
  * 5 Fert, A., Reyren, N. & Cros, V.  Magnetic skyrmions: Advances in physics and potential applications (2017).
  * 6 Rybakov, F. N., Borisov, A. B. & Bogdanov, A. N.  Three-dimensional skyrmion states in thin films of cubic helimagnets.  _Physical Review B - Condensed Matter and Materials Physics_ 87, 1–4 (2012).
  * 7 Meynell, S. A., Wilson, M. N., Fritzsche, H., Bogdanov, A. N. & Monchesky, T. L.  Surface twist instabilities and skyrmion states in chiral ferromagnets.  _Physical Review B - Condensed Matter and Materials Physics_ 90, 014406 (2014).
  * 8 Leonov, A. O., Loudon, J. C. & Bogdanov, A. N.  Spintronics via non-axisymmetric chiral skyrmions.  _Applied Physics Letters_ 109, 1–5 (2016).
  * 9 Schneider, S. _et al._ Induction Mapping of the 3D-Modulated Spin Texture of Skyrmions in Thin Helimagnets.  _Physical Review Letters_ 120, 217201 (2018).
  * 10 Birch, M. T. _et al._ Real-space imaging of confined magnetic skyrmion tubes.  _Nature Communications_ 11, 1–8 (2020).
  * 11 Yu, X. _et al._ Real-Space Observation of Topological Defects in Extended Skyrmion-Strings.  _Nano Letters_ 20, 7313–7320 (2020).
  * 12 Wolf, D. _et al._ Holographic vector field electron tomography of three-dimensional nanomagnets.  _Communications Physics_ 2, 1–9 (2019).
  * 13 Rößler, U. K., Bogdanov, A. N. & Pfleiderer, C.  Spontaneous skyrmion ground states in magnetic metals.  _Nature_ 442, 797–801 (2006).
  * 14 Rößler, U. K., Leonov, A. A. & Bogdanov, A. N.  Chiral skyrmionic matter in non-centrosymmetric magnets.  In _Journal of Physics: Conference Series_ , vol. 303, 12105 (Institute of Physics Publishing, 2011).
  * 15 Wilhelm, H. _et al._ Confinement of chiral magnetic modulations in the precursor region of FeGe.  _Journal of Physics Condensed Matter_ 24, 294204 (2012).
  * 16 Balkind, E., Isidori, A. & Eschrig, M.  Magnetic skyrmion lattice by the Fourier transform method.  _Physical Review B_ 99, 134446 (2019).
  * 17 Huang, P. _et al._ Melting of a skyrmion lattice to a skyrmion liquid via a hexatic phase.  _Nature Nanotechnology_ 15, 761–767 (2020).
  * 18 Yu, X. Z. _et al._ Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe.  _Nature Materials_ 10, 106–109 (2011).
  * 19 Wilson, M. N. _et al._ Extended elliptic skyrmion gratings in epitaxial MnSi thin films.  _Phys. Rev. B_ 86, 144420 (2012).
  * 20 Rybakov, F. N., Borisov, A. B., Blügel, S. & Kiselev, N. S.  New spiral state and skyrmion lattice in 3D model of chiral magnets.  _New Journal of Physics_ 18, 045002 (2016).
  * 21 Leonov, A. O. _et al._ Chiral Surface Twists and Skyrmion Stability in Nanolayers of Cubic Helimagnets.  _Physical Review Letters_ 117, 087202 (2016).
  * 22 Du, H. _et al._ Interaction of Individual Skyrmions in a Nanostructured Cubic Chiral Magnet.  _Physical Review Letters_ 120, 197203 (2018).
  * 23 Nagaosa, N. & Tokura, Y.  Topological properties and dynamics of magnetic skyrmions.  _Nature Nanotechnology_ 8, 899–911 (2013).
  * 24 Leonov, A. O. & Mostovoy, M.  Edge states and skyrmion dynamics in nanostripes of frustrated magnets.  _Nature Communications_ 8, 1–7 (2017).
  * 25 Lee, M., Kang, W., Onose, Y., Tokura, Y. & Ong, N. P.  Unusual hall effect anomaly in MnSi under pressure.  _Physical Review Letters_ 102, 186601 (2009).
  * 26 Neubauer, A. _et al._ Topological hall effect in the a phase of MnSi.  _Physical Review Letters_ 102, 186602 (2009).
  * 27 Zang, J., Mostovoy, M., Han, J. H. & Nagaosa, N.  Dynamics of Skyrmion crystals in metallic thin films.  _Physical Review Letters_ 107, 136804 (2011).
  * 28 Schulz, T. _et al._ Emergent electrodynamics of skyrmions in a chiral magnet.  _Nature Physics_ 8, 301–304 (2012).
  * 29 Dahir, S. M., Volkov, A. F. & Eremin, I. M.  Interaction of Skyrmions and Pearl Vortices in Superconductor-Chiral Ferromagnet Heterostructures.  _Physical Review Letters_ 122, 097001 (2019).
  * 30 Huang, S. X. & Chien, C. L.  Extended skyrmion phase in epitaxial FeGe(111) thin films.  _Physical Review Letters_ 108, 267201 (2012).
  * 31 Zheng, F. _et al._ Experimental observation of chiral magnetic bobbers in B20-type FeGe.  _Nature Nanotechnology_ 13, 451–455 (2018).
  * 32 Leonov, A. O., Tambovtcev, I. M., Lobanov, I. S. & Uzdin, V. M.  Stability of in-plane and out-of-plane chiral skyrmions in epitaxial MnSi(111)/Si(111) thin films: Surface twists versus easy-plane anisotropy.  _Phys. Rev. B_ 102, 174415 (2020).
  * 33 Jin, C. _et al._ Control of morphology and formation of highly geometrically confined magnetic skyrmions.  _Nature Communications_ 8, 15569 (2017).
  * 34 Kravchuk, V. P. _et al._ Multiplet of Skyrmion States on a Curvilinear Defect: Reconfigurable Skyrmion Lattices.  _Phys. Rev. Lett._ 120, 67201 (2018).
  * 35 Phatak, C., Petford-Long, A. K. & De Graef, M.  Three-dimensional study of the vector potential of magnetic structures.  _Physical Review Letters_ 104, 1–4 (2010).
  * 36 Tanigaki, T. _et al._ Three-dimensional observation of magnetic vortex cores in stacked ferromagnetic discs.  _Nano Letters_ 15, 1309–1314 (2015).
  * 37 Wolf, D. _et al._ 3D Magnetic Induction Maps of Nanoscale Materials Revealed by Electron Holographic Tomography.  _Chemistry of Materials_ 27, 6771–6778 (2015).
  * 38 Streubel, R. _et al._ Retrieving spin textures on curved magnetic thin films with full-field soft X-ray microscopies.  _Nature Communications_ 6, 1–11 (2015).
  * 39 Donnelly, C. _et al._ Three-dimensional magnetization structures revealed with X-ray vector nanotomography.  _Nature_ 547, 328–331 (2017).
  * 40 Hierro-Rodriguez, A. _et al._ Revealing 3D magnetization of thin films with soft X-ray tomography: magnetic singularities and topological charges.  _Nature Communications_ 11, 6382 (2020).
  * 41 Stolt, M. J. _et al._ Electrical Detection and Magnetic Imaging of Stabilized Magnetic Skyrmions in Fe1-xCoxGe (x< 0.1) Microplates.  _Advanced Functional Materials_ 29, 1805418 (2019).
  * 42 Parkin, S. S., Hayashi, M. & Thomas, L.  Magnetic domain-wall racetrack memory (2008).
  * 43 Kovács, A. _et al._ Mapping the magnetization fine structure of a lattice of Bloch-type skyrmions in an FeGe thin film.  _Applied Physics Letters_ 111, 192410 (2017).
  * 44 Leonov, A. O. _et al._ The properties of isolated chiral skyrmions in thin magnetic films.  _New Journal of Physics_ 18, 065003 (2016).
  * 45 Lebech, B., Bernhard, J. & Freltoft, T.  Magnetic structures of cubic FeGe studied by small-angle neutron scattering.  _Journal of Physics: Condensed Matter_ 1, 6105–6122 (1989).
  * 46 Song, D. _et al._ Quantification of Magnetic Surface and Edge States in an FeGe Nanostripe by Off-Axis Electron Holography.  _Physical Review Letters_ 120, 167204 (2018).
  * 47 Butenko, A. B., Leonov, A. A., Rößler, U. K. & Bogdanov, A. N.  Stabilization of skyrmion textures by uniaxial distortions in noncentrosymmetric cubic helimagnets.  _Physical Review B - Condensed Matter and Materials Physics_ 82, 052403 (2010).
  * 48 Beg, M. _et al._ Ground state search, hysteretic behaviour, and reversal mechanism of skyrmionic textures in confined helimagnetic nanostructures.  _Scientific Reports_ 5, 1–14 (2015).
  * 49 Bogdanov, A. & Hubert, A.  Thermodynamically stable magnetic vortex states in magnetic crystals.  _Journal of Magnetism and Magnetic Materials_ 138, 255–269 (1994).
  * 50 Richardson, M., Ingri, N., Salomaa, P., Bloom, G. & Hagen, G.  The Partial Equilibrium Diagram of the Fe-Ge System in the Range 40-72 at. % Ge, and the Crystallisation of some Iron Germanides by Chemical Transport Reactions.  _Acta Chemica Scandinavica_ 21, 2305–2317 (1967).
  * 51 Bosholm, O., Oppermann, H. & Däbritz, S.  Chemischer Transport intermetallischer Phasen IV: Das System Fe-Ge: Chemical Vapour Transport of Intermetallic Phases IV: The System Fe - Ge.  _Zeitschrift fur Naturforschung - Section B Journal of Chemical Sciences_ 56, 329–336 (2001).
  * 52 Boothroyd, C., Kovács, A. & Tillmann, K.  FEI Titan G2 60-300 HOLO.  _Journal of large-scale research facilities JLSRF_ 2, 44 (2016).
  * 53 Wolf, D., Lubk, A., Lichte, H. & Friedrich, H.  Towards automated electron holographic tomography for 3D mapping of electrostatic potentials.  _Ultramicroscopy_ 110, 390–399 (2010).
  * 54 Ghiglia, D. C., Ghiglia, D. C., Pritt, M. D. & Pritt, M. D.  _Two-Dimensional Phase Unwrapping - Theory, Algorithms, and Software_ (Wiley, New York, 1998).
  * 55 Lubk, A. _et al._ Nanometer-scale tomographic reconstruction of three-dimensional electrostatic potentials in GaAs/AlGaAs core-shell nanowires.  _Physical Review B - Condensed Matter and Materials Physics_ 90 (2014).
  * 56 Wolf, D., Lubk, A. & Lichte, H.  Weighted simultaneous iterative reconstruction technique for single-axis tomography.  _Ultramicroscopy_ 136, 15–25 (2014).


# Shear profile in a dense packing of large grains

Alessio Quaresimaa Andrea Platib,c,d Andrea Gnolic,d Alberto Petric,d,e

###### Abstract

We investigate the shear of a dense bed of refracting supermillimetric grains
confined within a transparent horizontal annular cell with a rotating top. The
local time correlation functions of interferometric images allow to
characterize the shear profile close to the wall, with a spatial resolution
well below the grain diameter. For increasing shear, we observe a transition
in the system response and the manifestations of anisotropies in the force
chains. The employed technique is of easy implementation and is especially
suitable for the study of stationary processes.

###### keywords:

Sheared granular bed , Light scattering , Shear band , Interferometry

††journal: Powder Technology

organization=Neurobiology of Language Department, Max Planck Institute for
Psycholinguistics, addressline=Wundtlaan 1, city=Nijmegen, postcode=6525 XD,
country=The Netherlands organization=Université Paris-Saclay, CNRS,
Laboratoire de Physique des Solides, addressline=1 rue Nicolas Appert,
city=Orsay, postcode=91405, country=Italy organization=Institute for Complex
Systems - CNR, addressline=P.le Aldo Moro 2, city=Rome, postcode=00185,
country=Italy organization=Department of Physics, University of Rome Sapienza,
addressline=P.le Aldo Moro 2, city=Rome, postcode=00185, country=Italy

organization=Enrico Fermi Research Center (CREF), addressline=via Panisperna
89, city=Rome, postcode=00184, country=Italy

##  1 Introduction

The response of a granular medium to shear stress is an important and not
completely known property [1], with many practical consequences in natural
phenomena and industry. Different behaviors can be observed, that depend on
many variables [2]. Besides the shear rate itself, these variables include
grain geometry and tribological properties, preparation, volume fraction,
boundary conditions (such as fixed volume or fixed normal stress, etc.) and
possibly others, through dependence that are not always well-understood. Not
only the strain rate profile can change, but highly nonlinear phenomena with
consequent instabilities can take place in different circumstances, like the
formation of shear sub-bands [3], stick-slip [4], etc. Accessing the material
bulk in 3D to measure the strain field is not impossible, but certainly hard
and resources demanding (see e.g. [5, 6, 7] end refs. therein). It can
therefore be useful to have at least a proxy [8] of how the granular medium
deforms and flows under shear in the different regions observing what happens
at the medium boundaries.

A tool that has proven particularly fruitful for the study of granular
dynamics is the scattering of light. Depending on the system and the
phenomenon at hand, different techniques have been developed. One widely
employed is Diffusing Wave Spectroscopy (DWS) [9]. It requires a transparent
medium, in which coherent light can undergo multiple scattering before
reaching the detector, giving rise to a diffusive process generating a grainy
pattern of speckles. DWS can be successfully employed to study the dynamics of
both submicrometer particles [9] and of larger, submillimeter, particles [10],
yielding information on the motion of grains within a certain depth from the
system boundaries. In larger particles, the diffusive process is generated by
multiple reflections and refraction at the grain interfaces.

In the case of supermillimetric particles, absorption prevents light from
being scattered a large number of times and from entering the diffusive
regime. In addition, the grain size can be comparable with the coherence
length of the laser beam. Thus, different methods have been adopted that allow
the study of grains adjacent to a transparent wall or to the system surface.
Among the most used are Particle Tracking Velocimetry (PTV) and Particle
Imaging Velocimetry (PIV) [11]. PTV consists of directly tracking the particle
trajectories of the grains. This technique can in principle attain high
precision, however, it requires that grains remain always visible, which may
be not verified in 3D systems, where they can migrate from the boundary to the
bulk and vice versa, and hence PTV is generally not suitable for long-time
measurements. PIV is based on the comparison of highly contrasted images taken
at different times. Deformation is evaluated by optimizing the overlap of the
images, which can be quite demanding [12]. PIV needs very bright images, often
obtained with the use of very powerful light sources ($\simeq 0.5$ kW) and a
mix of well-contrasting beads, e.g. black and white. Moreover, image
processing requires elaborated software, being performed usually in the
Fourier space. At the same time, for inhomogeneous systems, images must be
split into several windows to get the desired spatial resolution, with
drawbacks on the Fourier Transform. Another technique widely used is the
Digital Image Correlation (DIC) [13] which can give detailed results on the
strain field after preparing the medium to yield highly contrasted images, as
in PTV.

In this paper, we propose an experimental procedure, to the best of our
knowledge original, to investigate deformations and flows at the boundary of a
dense assembly of refracting grains. The investigated system consists of glass
beads contained in a horizontal transparent channel and sheared at the top.
The procedure exploits the interference pattern generated by reflection and
refraction at the air-grain interface and possibly by scattering from rough
grains [14], and is sensible to the motion of particles closer to the
container wall. Information on the vertical profile of the shear band is
obtained from the time correlation of the light intensity. The advantage of
the technique lays in its simple implementation and in its ability in
resolving the dynamics well below the grain size, yielding spatially detailed
information.

The aim of this work is two-fold. On the one hand to show the effectiveness of
the proposed technique. On the other hand to present the main, specific,
features characterizing the shear band in the investigated system. In Sec. 2,
we describe the experimental setup and the data processing. The main results
are reported in Sec. 3 while Sec. 4 contains a discussion and a summary of the
main findings. Several technical and methodological details can be found in
the Supplementary Material (SM), together with the URL where an illustrative
movie is stored, and a Julia code expressly developed to process the images is
made available.

##  2 Material and methods

###  2.1 Experimental apparatus

The granular system under investigation is similar to the one described in
[15] (see SM for details).

Figure 1:  Schematic sketch of the experiment.

Grains are confined in a transparent horizontal annular channel with inner and
outer radii $r=12.5$ cm and $R=19.2$ cm, respectively. To prevent
crystallization, it is composed of a $50\%-50\%$ mixture of hard glass beads,
1.5 and 2.0 mm in diameter. An aluminum plate lays on the top, in direct
contact with the granular material, and is connected to a motor so that it can
shear the system at constant angular speed $\omega$. A layer of grains is
glued to its lower side in order to increase drag. However, it does not extend
to the whole plate surface, but excludes an annulus close to the border, to
prevent free grains from getting stuck between glued grains and container
walls (see the section view in Fig. 1). It will be seen that this mechanical
feature is revealed by optical measurements. As sketched in Fig. 1,
interferometry is performed in backscattering through a single mode Nd:Yag
laser yielding elliptically polarized light at a wavelength of 532 nm, with a
maximum nominal power of 500 mW. The horizontal laser beam is filtered by a
polarizer and is expanded by a diverging lens. The back-scattered light is
then focused towards a perpendicular white screen by another lens and
eventually passes through a crossed polarizer, to minimize the amount of light
directly reflected by the surfaces of the container wall. The resulting
interference image is framed by a digital camera equipped with a CMOS array,
at a rate of 15 frames per second. The size of the probed region is limited by
at least three factors: a) the laser power, which restricts the scattering
area and the thickness of the granular medium involved. With the employed
equipment, a sufficient intensity contrast is still obtained with a spot of
about 20 mm in linear size, i.e. 10-13 bead diameters, independently of the
chosen optical magnification for the image; b) the coherence length of the
beam, which determines the optical depth of the region where interference can
be generated; c) the curvature of the container, which gives rise to optical
aberration and consequently requires a sufficient narrow imaged area. In order
to get a good focusing, the image has been cropped horizontally.

###  2.2 Data processing and analysis: coarse-graining and noise reduction

The size of the interference pattern depends on the width of the laser spot
and on the optical magnification $M$ of the scattered image, that can be
chosen opportunely. An example of the emerging pattern is shown in Fig. 2. The
image is made of 240 (horizontal) x 420 (vertical) pixels whose intensities
have been normalized as explained in the SM. One pixel corresponds to a region
of linear size 0.04 mm on the container wall. The surface of the material in
contact with the rotating plate coincides with the top of the image ($z=0$).

Figure 2: An interference image showing the spatial reference system of our
setup. The position of the rotating plate coincides with $z=0$.

To detect how the granular system deforms under different shear rates at
different depths, video records of the interference pattern have been shot at
a rate of 15 frames per second for about 11 minutes (see SM for more details).
So, each video supplies about 105 time series of 9.9k points each.

$\omega$ (s-1) | 0.05 | 0.07 | 0.1 | 0.5 | 0.7 | 1.0 | 1.2 | 1.5
---|---|---|---|---|---|---|---|---
$v$ (mm s-1) | 9.6 | 13.4 | 19.2 | 96 | 134 | 192 | 230 | 288
Table 1: Values of plate speed considered in the experiments in radiant per
second. The lower line reports the plate speed in mm s-1 at the outer wall of
the container.

However, the signal obtained from single pixels are generally highly noisy.
For this reason, it is often convenient to perform an average of the intensity
over larger regions, covering $S^{2}$ pixels, also called metapixels [16],
chosen on the base of the experimental conditions. In the present case the
size $S$ of the metapixel has been chosen through a procedure (detailed in the
SM) aimed at getting a trade-off between spatial resolution and intrinsic
noise, yielding $S=10$ pixels, corresponding to $\approx 1/4\div 1/2$ bead
diameters. Thus, the final interference images are composed of 24 (horizontal)
x 42 (vertical) metapixels. Each metapixel is identified by an integer
horizontal coordinate $x\in[1,24]$, from left to right, and a vertical
coordinate $z\in[1,42]$, from top to bottom.

For every value of the plate speed reported in Tab. 1 the instantaneous
intensity on each metapixel, $I(x,z,t)$ has been obtained averaging the
intensities on the corresponding metapixels; Then, the connected
autocorrelation functions have been computed on each metapixel as a function
of position and time; in the following equation, the variables x,z are omitted
for clarity, i.e., $I(x,z,t)=I(t)$.

$\displaystyle g_{\omega}(x,z,t)=\cfrac{\langle
I(t_{0})I(t_{0}+t)\rangle_{t_{0}}-\langle I(t_{0})\rangle^{2}_{t_{0}}}{\langle
I^{2}(t_{0})\rangle_{t_{0}}-\langle I(t_{0})\rangle^{2}_{t_{0}}}$ (1)

where, as indicated by the subscript, averages are made over different initial
times $t_{0}$. Their values at discrete times $t=n~{}\cdot~{}\Delta t$ with
$\Delta t=1/15$ $s$, have been obtained through the corresponding empirical
estimates (see SM). A purposely devised code with a parallelized structure has
been developed to process the data111The code is made available for public
use, see on-line SM.

Figure 3: Intensity matrices averaged over time for the entire duration of the
experiments. Images become more blurred for increasing plate speed (from top
left to bottom right).

##  3 Results

###  3.1 Phenomenology

Before discussing autocorrelation functions, we first qualitatively
characterize the phenomenology of the system by looking at the time-averaged
intensity matrix. We recall that different mechanisms may contribute to the
time evolution of scattered light intensity. On the one hand, grains slowly
diffuse, changing their relative position and thence the interference pattern
generated by refractions and reflections at the glass-air interfaces. On the
other hand, each single grain moves as a rigid body and can originate changes
in the pattern due to the light scattering occurring at its rough surface
[14].

Figure 3 shows the time average of the intensity matrix for some different
velocities of the driving plate. Such an average has been performed over the
entire experiment duration producing the same effect of a high exposure time
in a photo. Indeed, the image appears blurred in the regions where intensity
changes in time 222This effect is exploited in the so-called Speckle
Visibility Spectroscopy or Single Exposure Speckle Photography. In the
uppermost part (closer to the rotating plate) the image is blurred rather
homogeneously, while in the bottom area, an intensity pattern can be neatly
seen as for the system at rest. Homogeneously blurred regions reflect a fast
and disordered motion of the particles while spotted areas contain grains that
are almost still during the experiment. As expected, the region in which
flowing motion is visible increases in size for increasing shear velocity,
extending deeper. This allows to discriminate among different behaviours of
the interference figure at different depths and to identify a shear band. An
illustrative movie showing intensity over the entire time duration of the
experiment is linked to the on-line SM.

Figure 4: Autocorrelation value after 15 frames (1 second) for each metapixel,
for plate speed $\omega=1$ $s^{-1}$.

###  3.2 Shear band characterization

####  3.2.1 Autocorrelation decay

The above considerations can be made more quantitative by analyzing the time
behavior of the intensity autocorrelation Eq. (1) for different depth $z$ and
plate speed $\omega$. Figure 4 shows an example ($\omega=1$ $s^{-1}$) of how
fast intensity decorrelates at each meta-pixel. The color map encodes the
autocorrelation values $g_{\omega}(x,z,T)$ after $T=1$ second (15 frames). The
map shows that the metapixels decorrelation speed depends on the depth but not
on the horizontal position, as expected from the symmetry of the system (see
also the illustrative movie linked to the SM).

In the following analysis, being the system rotationally invariant, we have
considered the curves $g_{\omega}(z,t)$ obtained by averaging
$g_{\omega}(x,z,t)$ over $x$. Examples of the curves obtained in this way for
different speeds are shown in Fig. 5A at a fixed depth $z=10$ mm. As expected,
the decay is faster at larger speeds. In a similar fashion, the
autocorrelation at different depths slows down for increasing distance from
the plate. Examples are shown in Fig. 5B, for $\omega=0.1$ s-1.

Figure 5: A) Autocorrelation function for different driving velocities
$\omega$ at a depth $z=10$ mm. Curves are shown without statistical errors for
clarity. B) Autocorrelation functions at different depths $z$ and fixed
$\omega=0.1$ s-1. C): Comparison between autocorrelations of raw intensities
(C1) and the ones obtained from the same signals after a running average over
5 frames $\simeq$ 0.3 seconds (C2).

####  3.2.2 Shear profile: Best fit estimates

In order to characterize the shear profile we performed a more detailed
analysis of the autocorrelation decay. To this aim, we estimated a
characteristic decay time $\tau_{C}(z)$ by fitting the initial part of
$g_{\omega}(z,t)$. Considering that the decorrelation of the interference
pattern is caused by grains motion, we can use $\tau_{C}(z)$ as a proxy of the
typical time needed for the system to perform a local spatial rearrangement at
depth $z$. Then, the shear profile can be defined through the decorrelation
rate $\nu_{C}(z)=\tau^{-1}_{C}(z)$. Since the finite acquisition rate (i.e. 15
fps) can cause abrupt decorrelations at short times, the fitting procedure has
been performed on autocorrelation functions of smoothed intensity signals.
Such smoothing has been done with a running average over a time window $\tau$.
Examples of autocorrelation functions obtained from smoothed signals are
reported in the bottom panel of figure 5C. Their decay at short times can be
fitted by a Gaussian function $\propto\exp(-(t/\tau_{G}(z))^{2})$. However,
the fitted time $\tau_{G}(z)$ is not equivalent to the physical decorrelation
time $\tau_{C}(z)$ since the short time behavior of the autocorrelation
functions is affected by the smoothing procedure (which in turn introduces an
additional timescale $\tau$). The analytical expression of the autocorrelation
function obtained from a smoothed signal reveals that the physical decay time
depends on $\tau_{G}(z)$ and $\tau$ through the relation
$\tau_{C}(z)=2\tau^{2}_{G}(z)/\tau$ (details of the calculations are reported
in the SM). Then, fitting $\tau_{G}(z)$ and knowing $\tau$ allows to estimate
$\tau_{C}(z)$. The shear profiles $\nu_{C}(z)=\tau_{C}^{-1}(z)$ obtained from
autocorrelation functions of signals smoothed over $\tau=0.33$ s (5 frames)
are reported in Fig. 6. As expected, $\nu_{C}$ at a fixed depth displays
larger values at higher shear rates.

A first remarkable result is that the technique is able to resolve the
dynamics well below the grain size. In fact, while the image of one grain
spans 2-4 metapixels (see SM), $\nu_{C}$ is seen to significantly change on
the scale of the single metapixel. It is also interesting to note that, for
any imposed shear velocity $\omega$, the decorrelation rate $\nu_{C}$ stays
finite at all the observed distances from the shearing plate, being generally
larger for larger $\omega$. Therefore, although very slowly, the grains flow
at any depth. This makes it difficult to define the shear band in a strict
manner, however, it is seen that below the surface $\nu_{C}$ displays a local
maximum, after which it decays until an inflection point. One can
conventionally identify this point with the shear band limit $Z_{B}$ (see SM).

Figure 6: Characteristic inverse times $\nu_{C}=\tau_{C}^{-1}$, representative
of the granular strain rate from Gaussian fits of smoothed autocorrelations at
short times for the different $\omega$ after averaging on the horizontal
coordinate. Figure 7: A) Decorrelation rates from Gaussian fits $\nu_{C}$ at
different depths and plate speeds. B) Half-height decorrelation rates
$\nu_{H}$ at different depths and plate speeds.

Another feature displayed by the different curves in Fig. 6 is that $\nu_{C}$
does not attain its maximum $\nu_{C}^{max}$ just below the plate but at a
certain depth $Z_{max}$. This can be understood considering that, as described
in §2.1, the lower side of the plate, in contact with the grains, is smoother
close to the container walls, and therefore its dragging is less effective at
the boundary than in the bulk. The appearance of a maximum shear at a certain
distance from the plate is consistent with the transmission of the shear
stress through a network of force chains which results to follow an average
inclined path from the central region of the plate to the wall. At the same
time it may seem weird that, increasing the plate velocity, the region of
maximum shear rate deepens, as seen in Fig. 6. It is known [17, 18] that
systems develop anisotropy in the contact force network along the direction of
the shear so that stress propagates at an average angle with respect to the
vertical that increases as a function of the shear rate. Likely, anisotropy
affects force propagation in all directions and the observed increase of
$Z_{max}$ with $\omega$ suggests that force chains in the direction
transversal to the shear undergo an opposite fate, decreasing their angle with
respect to the vertical for increasing shear rate. We have not been able to
find related results in the literature, so this would be the first evidence of
such a mechanism, to be corroborated by additional experiments and
simulations.

A pictorial view of the shear band is shown in Fig. 7A. It consists of a color
map of the same $\nu_{C}$ reported in Fig. 6 after normalizing all values to
their maximum for each plate velocity. The vertical scale reports the deepness
$z$. It can then be noticed that for each speed there is an upper region where
$\nu_{C}$ is larger and changes more rapidly. This region extends deeper for
increasing plate speed evidencing a main shear band.

The behavior of the main quantities characterizing the shear band as a
function of the plate velocity is resumed in Fig. 8. Noticeably, both
quantities display a similar dependence. Two different linear regimes can be
envisaged, separated by a crossover velocity $\omega_{c}\approx 0.1$ s-1. In
the case of a compliant coupling between the motor and the plate, $\omega_{c}$
corresponds to a transition from plate stick-slip to continuous sliding [15].
This transition does not depend on the stiffness of the coupling [4] but is
intrinsic to the medium response to shear in this system. Therefore it gives
rise to the appearance of different regimes even in the present case of a
rigid coupling where the plate cannot perform stick-slip.

Figure 8: Characterizing quantities for the shear band as functions of
$\omega$. The automatized procedure to detect such regions is explained in the
SM together with the definition of the error bars for $Z_{B}$.

####  3.2.3 Shear profile: Fast characterization

For many practical applications, it can be desirable to have a fast way for
characterizing the autocorrelation decays, possibly avoiding performing many
best fits to all the curves. To this aim, we have evaluated the half-height
decorrelation rate $\nu_{H}=\tau_{H}^{-1}$, where $\tau_{H}$ is the time at
which the value of $g_{\omega}(z,t)$ is halved (we point out that $\nu_{C}$ is
obtained from autocorrelations of smoothed signals while $\nu_{H}$ is not). As
a first qualitative comparison with the best fits, we report in Fig. 7B the
shear-band colormap obtained with $\nu_{H}$, which exhibits a phenomenology in
good agreement with Fig. 7A. Indeed, for each omega, $\nu_{H}$ is able to
distinguish between the slow-flowing region at large depth and the fast-
flowing one at lower $z$. However, within this latter region, $\nu_{H}$
reveals to be less accurate than $\nu_{C}$ because of the finite acquisition
rate (see also Fig. S4 in the SM). A more quantitative test can be done by
comparing the depth of the maximum decorrelation rate obtained from $\nu_{H}$
with the one obtained from the best-fit parameter $\nu_{C}$ for each value of
$\omega$. Figure 9 reports a regression between the two quantities, showing a
good degree of correlation. The half-life times can therefore be employed to
characterize the shear band in a fast and simple way.

Figure 9: Correlation between the height of the characteristic maxima of
decorrelation rates $\nu_{H}$ and $\nu_{C}$.

##  4 Conclusions

We have experimented an interferometry-based imaging technique to investigate
the shear in a dense packing of supermillimetric transparent grains confined
in a transparent channel and sheared by a rotating top lid. The technique is
of easy implementation and capable of resolving the shear profile well below
the grain size. Different characterizing parameters can be extracted from the
analysis of local intensity time autocorrelation functions. Their similar
dependence on the plate speed shows that for the investigated system they are
equivalent. The technique has also allowed to detect the existence of two
shear regimes, separated by a crossover which corresponds to the transition
from granular stick-slip to continuous sliding. Shifting of the maximum shear
to larger depths for increasing plate speed reveals the development of
anisotropies in the force chains. Our comparative analysis of the parameters
shows that a quick and robust characterization of the shear profile can also
be obtained straightly from the half-life time of the correlation functions.
The technique can potentially be employed in different situations, especially
for the study of stationary processes, and more internal layers could be
probed increasing the coherence length and power of the light source.

##  5 Acknowledgments

The authors acknowledge R. Cerbino, M. Leonetti and G. C. Ruocco for useful
discussions.

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## References

  * [1] Yoël Forterre and Olivier Pouliquen.  Flows of dense granular media.  Annual Review of Fluid Mechanics, 40(1):1–24, 2008.
  * [2] GDR MiDi.  On dense granular flows.  The European Physical Journal E, 14, 2004.
  * [3] Priyanka Shukla, Lima Biswas, and Vinay Kumar Gupta.  Shear-banding instability in arbitrarily inelastic granular shear flows.  Phys. Rev. E, 100:032903, Sep 2019.
  * [4] F. Dalton, F. Farrelly, A. Petri, L. Pietronero, L. Pitolli, and G. Pontuale.  Shear stress fluctuations in the granular liquid and solid phases.  Phys. Rev. Lett., 95:138001, Sep 2005.
  * [5] Xiaoyu Yang, Chao Huan, D. Candela, R. W. Mair, and R. L. Walsworth.  Measurements of Grain Motion in a Dense, Three-Dimensional Granular Fluid.  Physical Review Letters, 88(4):4, 2002.
  * [6] Axelle Amon, Philip Born, Karen E. Daniels, Joshua A. Dijksman, Kai Huang, David Parker, Matthias Schröter, Ralf Stannarius, and Andreas Wierschem.  Preface: Focus on imaging methods in granular physics.  Review of Scientific Instruments, 88(5), 2017.
  * [7] H. T. Fabich, T. I. Brox, D. Clarke, J. D. Seymour, S. L. Codd, P. Galvosas, J. Brown, A. J. Sederman, and D. J. Holland.  Measurements of the velocity distribution for granular flow in a couette cell.  Phys. Rev. E, 98:062901, Dec 2018.
  * [8] L. Bocquet, W. Losert, D. Schalk, T. C. Lubensky, and J. P. Gollub.  Granular shear flow dynamics and forces: Experiment and continuum theory.  Phys. Rev. E, 65:011307, Dec 2001.
  * [9] D. A. Weitz and D. J. Pine.  Diffusing-wave spectroscopy, chapter Dynamic Light Scattering: The Method and Some Applications.  Oxford University Press, 1993.
  * [10] A. Amon, A. Mikhailovskaya, and J. Crassous.  Spatially-resolved measurements of micro-deformations in granular materials using dws.  Review of Scientific Instruments, 88, 2017.
  * [11] T. Hagemeier, M. Börner, A. Bück, and T. Evangelos.  A comparative study on optical techniques for the estimation of granular flow velocities.  Chemical Engineering Science, pages 63–75, 2015.
  * [12] L. Sarno, Y-C Tai, A. Carravetta, R. Martino, M. N. Papa, and C-Y Kuo.  Challenges and improvements in applying a particle image velocimetry (piv) approach to granular flows.  J. Phys.: Conf. Ser., 1249:012011, 2019.
  * [13] Dong Lei, Jinfeng Huang, Wenxiang Xu, Wenchao Wang, and Pei Zhang.  Deformation analysis of shear band in granular materials via a robust plane shear test and numerical simulation.  Powder Technology, 383:385–392, 2018.
  * [14] Lisa Dossow, Raphael Kessler, Matthias Sperl, and Philip Born.  Dynamic light scattering from single macroscopic particles.  Applied Optics, 60:10160–10167, 2021.
  * [15] A. Baldassarri, M. A. Annunziata, A. Gnoli, G. Pontuale, and A. Petri.  Breakdown of scaling and friction weakening in intermittent granular flow.  Scientific Reports, 9:16962, 2019.
  * [16] Virgile Viasnoff, François Lequeux, and D. J. Pine.  Multispeckle diffusing-wave spectroscopy: A tool to study slow relaxation and time-dependent dynamics.  Review of Scientific Instruments, 73(6):2336, 2002.
  * [17] R. Behringer and B. Chakraborty.  The physics of jamming for granular materials: a review.  Reports on Progress in Physics, 82:012601, 2019.
  * [18] Fanjing Meng, Huabo Liu, Shaozhen Hua, and Minghua Pang.  Force chain characteristics of dense particles sheared between parallel-plate friction system.  Results in Physics, 25:104328, 2021.

## Supplemental Material

## Experimental setup

### Mechanics

The mechanical experimental setup used for this research (Fig. S1) is based on
a circular horizontal PPMI channel of outer and inner radii $R$ = 19.2 cm and
$r$ = 12.5 cm respectively, and 12 cm height. A bidisperse mixture of
$50\text{\,}\%$ \- $50\text{\,}\%$ glass beads, of radii $r_{1}=$
$1.5\text{\,}\mathrm{mm}$ $\pm$ $10\text{\,}\%$ and
$r_{2}=$$2\text{\,}\mathrm{mm}$ $\pm$ $10\text{\,}\%$, fills the channel to
about $9\text{\,}\mathrm{cm}$ in height. A horizontal plate fitting the
channel can shear the granular medium from the top. It is free to move
vertically and can be rotated. A layer of grains is glued to its lower face in
order to better drag the underlying granular medium. It does not extend to the
whole plate surface, but excludes an annulus close to the border, to prevent
free grains from getting stuck between glued grains and container walls (Fig.
S10). The plate has mass $M=$$1200\text{\,}\mathrm{g}$ and moment of inertia
$I=$$0.026\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{2}$, implying that in our
experiments the medium is under a nominal pressure
$p=Mg/[\pi(R^{2}-r^{2})]\simeq$ $176\text{\,}\mathrm{Pa}$. The plate can be
rotated through a stepping motor, at constant angular speed $\omega_{p}$ in
the interval of
$0.05\text{\,}\mathrm{rad}\mathrm{/}\mathrm{s}1.5\text{\,}\mathrm{rad}\mathrm{/}\mathrm{s}$.

### Optics

The optical setup is also sketched in Fig. S10. A Nd:YAG laser supplies an
elliptically polarized light beam at a wavelength of
$532\text{\,}\mathrm{nm}$, with a maximum nominal power of
$500\text{\,}\mathrm{mW}$. The beam is filtered by a polarizing filter and is
expanded by a lens with $f=50$ mm, hitting the channel wall horizontally. The
back-scattered light is then focused towards a perpendicular white screen by
another $f=50$ mm lens. To minimize the amount of light reflected by the
perspex walls and by the first vertical layer of grains, the beam eventually
passes through another polarizer, crossed with the first.

Figure 10: Figure S10: a) Sketch of the experimental set-up with a detail of
the top plate.

The resulting interference figures have been recorded with a D90 Nikon camera
equipped with a CMOS, of 640$\times$424 pixels in the horizontal and vertical
direction respectively, and a dynamics of 8 bits. A lens with $f=50$ mm has
been employed, with an acquisition rate of
$15\text{\,}\text{frames}\mathrm{/}\mathrm{s}$. Due to the channel curvature
and the laser spot dimension, the useful area of the interference pattern on
the detector had a size of 250$\times$420 pixels. Calibration measures have
shown that one pixel corresponds to a region of the size of
$0.04\text{\,}\mathrm{mm}$ on the container wall, so the analyzed region of
the granular system has an external surface of 10 x 16.8 . The grains diameter
is, on average, of about $1.75\text{\,}\mathrm{mm}$, corresponding to
$44\text{\,}\text{pixels}$. The main typical linear sizes characterizing our
system are resumed in Tab. S1.

element | mm | pixel
---|---|---
pixel | 0.04 | 1
grain | 1.5 - 2.0 | 38 - 50
interference figure (width $\times$ height) |  10.0 x 16.8 |  250x420
metapixel | 0.4 | 10
Table 2: Table S1: Some sizes characterizing the system

## Preprocessing of video

##### Scaling of pixel intensities.

A video of duration $11\text{\,}\mathrm{min}$ has been recorded for different
speeds of the rotating plate. Each video is encoded as a three-dimensional
matrix whose indexes span the two spatial dimensions and the time. This way,
the time evolution of the interference pattern is contained in a three-
dimensional array of dimensions (250, 420, 9900) whose entries represent the
back-scattered light intensity recorded at a pixel $(x,z)$ at a time $t$.
Because the beam intensity is not uniform over the channel wall, one could
normalize data with the average intensity of each pixel. However, the presence
of very bright spots in the interference pattern makes this kind of
normalization ineffective. Thus we proceeded by taking a reference value at
each pixel which has been computed in the following way. First, we averaged
all the intensities over time obtaining a two-dimensional array (250, 420).
Then, we performed averages of this array over rows and columns obtaining two
one-dimensional arrays. Finally, the normalization matrix was obtained through
the external product of such vectors and taking the square root of each entry
(see Fig. S11A).

##### Spatial coarse-graining and choice of the metapixel.

The intensity data collected at every pixel have been averaged in space over
squares of side $S=$ $10$ pixels, the metapixel (MP). To determine the optimal
size, a set of $5\text{\,}\mathrm{min}$ long data has been collected with the
plate at rest and the following steps followed:

  * i)

Coarse-graining was performed on each frame by averaging intensities over the
pixels of a square block of size $S$. At the beginning $S=$$1$;

  * ii)

for each block, the variance $\sigma^{2}$ over all the frames of the video was
computed;

  * iii)

The computed variances were averaged over the blocks of the entire
interference image, yielding a value $\langle\sigma^{2}(S)\rangle$ for each
choice of $S$

  * iv)

the block size $S$ was increased by one and the procedure repeated.

Notice that in the absence of noise on the detector, each variance should be
strictly zero since we are considering the case with the plate at rest. Thus,
the average variance $\langle\sigma^{2}(S)\rangle$ indicates the amount of
background noise; being uncorrelated, this noise is expected to decrease for
increasing sizes of the metapixel $S$. The results obtained in our setup are
shown in Fig. S11B). Panel B shows that the standard deviation decreases as
$S^{-1}$ as expected. However, the noise reduction obtained in this way is
paid for by a loss of spatial resolution. Indeed, once we construct an
interference image with metapixels of size $S$ we cannot resolve anymore any
spatial modulation of the intensity over scales $<S$. This is a crucial point
since we are interested in the characterization of the shear profile. The
criterion we used to stop the coarse-graining procedure is to choose an
optimal $S$ which is smaller than the typical spatial correlation length of
the intensity pattern in the video with the plate at rest. In this way, we
minimize the resolution loss since each metapixel contains on average only
correlated pixels which bring the same local information. This procedure
allows us to only filter the noise around local intensity values without
losing information on spatial intensity modulations. In Fig. S11C we show the
spatial correlation functions of the intensity pattern for $\omega=0$
$s^{-1}$. We compared the ones obtained considering distances over both $x$
and $y$, computed on single frames and then averaged over 100 frames. From
this plot, it is seen that a reasonable choice for the metapixel size is
$S=10$ since intensities within this distance are correlated in both
directions.

Figure 11: Figure S11: A) Normalization matrix obtained from the external
product of intensities averaged over rows and columns (see the text for the
specific procedure). B) Total time variance for $\omega=0$, as a function of
the metapixel size $S$; C) Spatial correlation functions of the intensity
pattern for $\omega=0$ on the $x$ and $y$ directions.

## Autocorrelation function

For each value of the plate rotational speed, we computed the connected time
autocorrelation function for the metapixel centred in $(x,z)$ (horizontal and
vertical coordinates):

$g(x,z,t)=\frac{\langle I(x,z,t_{0})I(x,z,t_{0}+t)\rangle_{t_{0}}-\langle
I(x,z,t_{0})\rangle^{2}_{t_{0}}}{\langle
I^{2}(x,z,t_{0})\rangle_{t_{0}}-\langle I(x,z,t_{0})\rangle^{2}_{t_{0}}}$ (2)

averaged over discrete initial times $t_{m}=m\Delta t$ as

$\displaystyle\langle
I(x,z,t_{0})I(x,z,t_{0}+t)\rangle_{t_{0}}=\frac{\sum_{m=1,N_{t}}I(x,z,t_{m}+t)I(x,z,t_{m})}{N}$
(3) $\displaystyle\langle
I^{2}(x,z,t_{0})\rangle_{t_{0}}=\sum_{m=1,N}I^{2}(x,z,t_{m})/N_{t}$
$\displaystyle\langle
I(x,z,t_{0})\rangle^{2}_{t_{0}}=\sum_{m=1,N}I(x,z,t_{m})/N,$ (4)

with $\Delta t=1/15$, $N=9.9\cdot 10^{3}$, corresponding to a time window of
$11\text{\,}\mathrm{min}$, and $N_{t}=N-\lceil t/\Delta t\rceil$.

The number of operations required for computing the autocorrelation function
grows at least as $N\log N$ if it is done via a fast Fourier transform. In the
present case, it was preferred a straight computation, even if it grows as
$N^{2}$, in order to avoid other drawbacks related to the calculation of the
Fourier transform of a finite series. To this aim, we developed a Julia
package named interference.jl described below. The autocorrelation function is
then obtained by first averaging the intensities over each metapixel, thus
reducing to a matrix of dimension (25,42,9900), and then as direct products of
3-d matrices (see Sec. 6). Being the system rotationally invariant, the
$g(x,z,t)$ have been averaged over $x$ ($x\in[1,25]$), resulting in the
$g(z,t)$, reported in the main text.

## Smoothing and fitting procedure

Using raw intensity matrices over time, we often observed a sharp jump between
the first and the second point of the autocorrelation function for some
metapixels, suggesting the presence of some residual noise faster than our
acquisition rate (15 fps). In order to address this issue, the inverse times
$\nu_{C}=\tau_{C}^{-1}$, shown in Fig. 5, have been estimated by fitting with
the positive part of a Gaussian $\propto\exp{(-2t^{2}/\tau\tau_{C})}$ the
short-time regime of autocorrelation functions obtained from smoothed
intensity signals. Here $\tau_{C}$ is the characteristic physical time related
to the autocorrelation decay while $\tau$ is the time window duration of the
running average used to smooth the signals. In the following, we provide a
simple mathematical argument to justify this procedure.

The smoothed intensities can be expressed as

$I_{s}(t)=\frac{1}{\tau}\int_{t}^{t+\tau}dt^{\prime}I(t^{\prime})$ (5)

and the non-normalized smoothed autocorrelations read

$\tilde{g}_{s}(t)=\frac{1}{T}\int_{0}^{T}dt^{\prime}I_{s}(t^{\prime})I_{s}(t^{\prime}+t).$
(6)

Here we assume $I(t)$ to be the signals with the time averages subtracted. To
simplify it, we avoid including the spatial and the $\omega$ dependence in the
notation. Let us now show analytically how this smoothing suppresses the fast-
decaying spurious contribution and why using a Gaussian fit for
$\tilde{g}_{s}(t)$ is reasonable. Using Eq. (5) we can write:

$\begin{split}\tilde{g}_{s}(t)=\frac{1}{T}\frac{1}{\tau^{2}}\int_{0}^{T}dt^{\prime}\int_{t^{\prime}}^{t^{\prime}+\tau}dt^{\prime\prime}I(t^{\prime\prime})\int_{t^{\prime}+t}^{t^{\prime}+t+\tau}dt^{\prime\prime\prime}I(t^{\prime\prime\prime})\\\
=\frac{1}{\tau^{2}}\int_{0}^{\tau}dt_{1}\int_{0}^{\tau}dt_{2}\frac{1}{T}\int_{0}^{T}dt^{\prime}I(t_{1}+t^{\prime})I(t_{2}+t^{\prime}+t)\end{split}$
(7)

where we changed variables as $t_{1}=t^{\prime\prime}-t^{\prime}$ and
$t_{2}=t^{\prime\prime\prime}-t^{\prime}-t$. The last integral of Eq. (7) can
be expressed in terms of the original (i.e. non smoothed) autocorrelation
function $\tilde{g}(t)$ only if the argument of one of the two integrated
intensities is larger than the other. Indeed we have that

$\frac{1}{T}\int_{0}^{T}dt^{\prime}I(t_{1}+t^{\prime})I(t_{2}+t^{\prime}+t)=\begin{cases}\tilde{g}(t_{1}-t_{2}-t)\quad\text{if}\quad
t_{1}\geq t_{2}+t\\\ \tilde{g}(t_{2}+t-t_{1})\quad\text{if}\quad t_{1}\leq
t_{2}+t\end{cases}\quad\text{with}\quad t>0$ (8)

From the above equation we see that the smoothed autocorrelation can be
written as a sum of two contributions:

$\tilde{g}_{s}(t)=\frac{1}{\tau^{2}}\int_{0}^{\tau}dt_{2}\int_{t_{2}+t}^{\tau}dt_{1}\tilde{g}(t_{1}-t_{2}-t)+\frac{1}{\tau^{2}}\int_{0}^{\tau}dt_{1}\int_{t_{1}-t}^{\tau}dt_{2}\tilde{g}(t_{2}+t-t_{1})$
(9)

For $t\leq\tau$, which is representative of the short-time decay of
$\tilde{g}_{s}$, we have to consider both contributions while for $t\geq\tau$
only the second one holds. This is because $t_{1}\in[0,\tau]$ so if
$t\geq\tau$ the inequality $t_{1}\geq t_{2}+t$ is never true.

In order to go on with the calculations, we need to assume a functional form
for $\tilde{g}(t)$. The simplest choice is taking exponentially decaying
correlations, with a characteristic time $\tau_{C}$:
$\tilde{g}(t)=A\exp(-t/\tau_{C})$. It coincides with modelling the intensity
signals $I(t)$ as Ornstein-Uhlenbeck processes. We make this assumption and
later verify _a posteriori_ his compatibility with empirical observations. By
substituting exponential autocorrelations into Eq. (9) we get

$\tilde{g}_{s}(t)=2A\frac{\tau_{C}}{\tau}\left[1-\frac{\tau_{C}}{\tau}\left(1-e^{-\tau/\tau_{C}}\right)\cosh(t/\tau_{C})\right]\quad\text{for}\quad
t<\tau.$ (10)

The first thing to note is that the above expression is proportional to
$\tau_{C}/\tau$ when this ratio is small. Thus, if the total signal contains a
statistically independent fast spurious component with a characteristic time
$\tau_{C}^{\prime}\ll\tau\ll\tau_{C}$, then the autocorrelation of the
smoothed signal will be dominated by the component proportional to
$\tau_{C}/\tau$. Another important feature of Eq. (10) is that it presents a
zero first derivative and a negative second derivative for $t\sim 0$,
consistent with the behavior observed in Fig. 4 of the main text (compare
panel C1 and C2). Considering now the normalized autocorrelation function
$g_{s}(t)=\tilde{g}_{s}(t)/\tilde{g}_{s}(0)$ and expanding it at the leading
order we find:

$\tilde{g}_{s}(t)=1-\frac{2t^{2}}{\tau\tau_{C}}.$ (11)

which is equivalent to the second-order expansion of the Gaussian
$\exp(-(t/\tau_{G})^{2})$, where $\tau_{G}=\sqrt{\tau\tau_{C}/2}$. Such result
justifies the use a Gaussian to fit the short-time part of the experimental
autocorrelation function and explains how $\tau$ and $\tau_{C}$ appear in it.
In principle one could perform the fit directly with Eq. (10), but we have
found more efficient the use of a Gaussian since it properly interpolates the
$t\sim\tau$ regime of Eq. (9) where Eq. (10) ceases to be valid.

## Characteristic time and shear band limit detection

For each considered plate velocity, $\omega_{p}$, the characteristic
decorrelation time at each depth $z$, $\tau_{C}(z,\omega_{p})$, has been
determined with two methods, obtaining compatible values.

### Gaussian fit

The method that yields the most accurate results is analyzing the
autocorrelation function of the smoothed intensity signal. The short-time
decay in the autocorrelation is fit with the positive part of a Gaussian
function (coherent with Eq. 11), which allows extracting the characteristic
time as a function of the metapixel height, $\tau_{C}(z)$, for each of the
eight plate speeds considered. In the analysis, the inverse of the
characteristic time $\nu_{C}(z)$ has been used because less noisy.

We have used a least-square method from the module LsqFit of Julia programming
language to fit the Gaussian. The code implements the Levenber-Marquardt
algorithm for non-linear fitting. The fitting function is a Gaussian with only
the characteristic time $\tau_{G}=\sqrt{\tau\tau_{C}/2}$ as a parameter:

$G(x,\tau_{C},\tau)=\exp\left[-2\left(\frac{x}{\sqrt{\tau_{C}\tau}}\right)^{2}\right]$
(12)

where the $\tau$ is the smoothing time, as discussed in the previous section,
set equal to 5 (equivalent to $1/3$ s). This allows extracting the
characteristic $\tau_{C}(z)$, for each of the eight plate speeds considered.
Its inverse, $\nu_{C}(z)$, expresses the decorrelation rate and has been
employed in the analysis

##### Shear band detection.

The decorrelation rates $\nu_{C}(z)$ allow to individuate the shear band
characterizing the quasi solid and liquid phases. As seen in Fig. 5 of the
main text, following the first maximum, $\nu_{C}(z)$ first decays and then
enters a plateau associated with the solid phase. We assume the shear band as
the region between the maximum of $\nu_{C}(z)$ and the beginning of the
plateau. The starting of the shear band is easy to localize. Conversely, to
individuate the plateau is less straightforward because of the noise in the
data. We determine, within a certain approximation, the beginning of the
plateau by associating it with the first relative minimum of $\nu_{C}(z)$
after the maximum. Assuming the curve to be smooth, to find the first minimum,
one could compute the finite differences of the decorrelation rates
$\Delta\nu_{C}(z)=\nu_{C}(z+1)-\nu_{C}(z)$, and look for the first point where
the difference is non-negative. However, the data is noisy, and the finite
differences can be non-negative also before the actual minimum in the data.
Hence, we averaged the finite differences over a window of size $w$, and
looked for the first point where the average is positive. A good estimate is
obtained by averaging over all the window lengths in the range
$3\text{\,}\text{metapixels}13\text{\,}\text{metapixels}$. The maxima and the
minima individuated with this method are shown in Fig. S12, together with the
corresponding $\nu_{C}(z)$ profiles. The error bars in Fig 8 of the main text
represent the standard deviation across the seven window lengths considered.

Figure 12: Figure S12: Profiles of the characteristic time for all the plate
velocities considered. The minima and maxima of the inverse characteristic
time are shown as stars and diamonds, respectively. The shear band is defined
as the region between the minimum and the first maximum. The data shown here
are the same used for Fig. 6 of the main text.

### Half-height method

A second method to analyze the data was also tested; this method requires less
analytical insight into the smoothing procedure and can be performed without
fitting the auto-correlation decays. To estimate the characteristic time, we
used the half-height method. The characteristic time $\tau_{H}(z)$ is when the
(raw) auto-correlation function decays to half its initial value. The half-
height method is less accurate than the Gaussian fit and more sensitive to
noise; however, it can be performed on non-smoothed data and is therefore
simpler and virtually unbiased by any timescale introduced for the analysis.
The profiles of $\tau_{H}(z)$ and their inverse $\nu_{H}(z)$ are shown in Fig.
S13. Because of the relatively coarse temporal resolution (determined by the
camera frame rate), the resolution of the half-height method saturates on the
short timescales associated with the beginning of the shear band.
Consequently, the profile of $\nu_{H}(z)$ has a plateau near the maximum, and
it is impossible to precisely detect the upward shear band limits. The right
panel in Fig. S13 illustrates this limitation.

Figure 13: Figure S13: Characteristic time (and its reciprocal in the right
panel) computed with the half-height method. The right panel better shows the
saturation of the method in the short timescales corresponding to the start of
the shear-band.

However, it is possible to compare the two methods and verify whether they
yield compatible results. The scatter plots in Fig. 14 illustrate the
characteristic times and their inverse for each metapixel and each plate
velocity, showing that there is good agreement between the two measures.

Figure 14: Figure S14: Comparison of the decorrelation rate computed with the
Gaussian fit $\nu_{C}$, and the one computed with the half-height method
$\nu_{H}$. The two methods yield compatible results. The diamond indicates the
maximum value of the decorrelation rates, which is used to define the starting
of the shear band. For large shear rates, $\nu_{H}$ saturates (see Fig. S13B)

## Code for interference image analysis

The code interference.jl has been developed for the experiment itself and is
available on GitHub333https://github.com/aquaresima/interference.jl. The
software uses bindings to Python Numpy to transform the `.avi` recordings into
`.npy` matrices, which are then imported in Julia. We preferred to use Julia
because the computation of the coarse-graining and the correlation is faster.
The original matrices are sized 250x420 $\times$ T, where T is not fixed. The
first two dimensions are the corresponding width and height pixel number of
the single frame, the third dimension is the number of frames. Each element of
the matrix represents the intensity at the corresponding pixel encoded in 8
bits - thus spanning from 0 to 255.

The code interference.jl offers a set of tools for data processing.

##### Signal coarse-graining and autocorrelation computation.

The matrices are coarsened in time and space with a non-overlapping average-
pooling algorithm - i.e. the value of coarsen matrix is the average of pixels’
intensity in the corresponding region of the fine grain matrix. The spatial
coarsen reduces raw signal noise, the coarsen in time also removes
uninformative correlations on time scales that are irrelevant for the present
study. The choice of the coarsening size is left to the user. The coarsen
matrices are then stored in `Float32` format to limit the memory occupation.

From the obtained three-dimensional matrix, the autocorrelation functions are
then computed as direct products of 3D matrices realized in a parallelized
way.

##### Estimation of the correlation time

The time correlation function is computed at each metapixel within the
assumption of stationary dynamics. Being $v_{ij}(t)=v_{ijt}$ the intensity at
each metapixel, the correlation algorithm returns a tridimensional matrix:

$\hat{C}_{i,j,s}=T^{-1}\sum_{t=0}^{T}\frac{v_{ij}(t)v_{ij}(t+s)-T^{-1}\sum_{t=0}^{T}v_{ij}^{2}(t)}{T^{-1}\sum_{t=0}^{T}v_{ij}^{2}(t)-(T^{-1}\sum_{t=0}^{T}v_{ij}(t))^{2}}.$
(13)

The time dimension T -the maximum correlation time- is determined by the
available number of frames $N_{f}$ and the maximum lag $s$: $T=N_{f}-s$. Thus,
each matrix element corresponds to the intensity autocorrelation at the
metapixel $i,j$ at the discrete time $s$. Assuming translational symmetry on
the $x$-axis, the correlation function is averaged over the horizontal
coordinate $i$ to obtain a 2-d matrix of scalar functions of $s$,
$\hat{C}_{j,s}$.

The characteristic decorrelation time $\tau_{j}$ is then determined with two
methods:

  * 1.

Half-height characteristic time: A decaying characteristic time scale of each
height is determined as the value of $s$ at which the function assumes half of
its initial value (i.e. 0.5). One is then left with a vector of characteristic
times: $\hat{\tau}_{j}$.

  * 2.

Gaussian characteristic time: The first 10 points of the correlation function
are fitted with a Gaussian function. The characteristic time is then computed
as described in Eqs. (12).

##### Shear band detection.

The detection of the jump between solid and fluid state, respectively left and
right in Fig. S12, corresponds to detect the maximum and first minimum of the
function $\tau_{i}$. The algorithm described in Sec. Characteristic time and
shear band limit detection is implemented in the code.