# Temperature Dependence of Magnetic Anisotropy and Domain Wall Tuning in
BaTiO3(111)/CoFeB Multiferroics

R. G. Hunt  School of Physics and Astronomy, University of Leeds, LS2 9JT,
United Kingdom  Bragg Centre for Materials Research, University of Leeds, LS2
9JT, United Kingdom  K. J. A. Franke  School of Physics and Astronomy,
University of Leeds, LS2 9JT, United Kingdom  P. S. Keatley  Department of
Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL,
United Kingdom  P. M. Shepley  School of Physics and Astronomy, University of
Leeds, LS2 9JT, United Kingdom  Bragg Centre for Materials Research,
University of Leeds, LS2 9JT, United Kingdom  M. Rogers  School of Physics and
Astronomy, University of Leeds, LS2 9JT, United Kingdom  Bragg Centre for
Materials Research, University of Leeds, LS2 9JT, United Kingdom  T. A. Moore
School of Physics and Astronomy, University of Leeds, LS2 9JT, United Kingdom
Bragg Centre for Materials Research, University of Leeds, LS2 9JT, United
Kingdom

###### Abstract

Artificial multiferroics consist of two types of ferroic materials, typically
a ferroelectric and ferromagnet, often coupled interfacially by
magnetostriction induced by the lattice elongations in the ferroelectric. In
BaTiO3 the magnitude of strain induced by these elongations is heavily
temperature dependent, varying greatly between each of the polar crystal
phases and exerting a huge influence over the properties of a coupled magnetic
film. Here we demonstrate that temperature, and thus strain, is an effective
means of controlling the magnetic anisotropy in BaTiO3(111)/CoFeB
heterostructures. We investigate the three polar phases of BaTiO3: tetragonal
(T) at room temperature, orthorhombic (O) below 280 K and rhombohedral (R)
below 190 K, across a total range of 77 K to 420 K. We find two distinct
responses; a step-like change in the anisotropy across the low-temperature
phase transitions, and a sharp high-temperature reduction around the
ferroelectric Curie temperature, measured from hard axis hysteresis loops.
Using our measurements of this anisotropy strength we are then able to show by
micromagnetic simulation the behaviour of all possible magnetic domain wall
states and determine their scaling as a function of temperature. The most
significant changes occur in the head-to-head domain wall states, with a
maximum change of 210 nm predicted across the entire range effectively
doubling the size of the domain wall as compared to room temperature. Notably,
similar changes are seen for both high and low temperatures which suggest
different routes for potential control of magnetic anisotropy and elastically
pinned magnetic domain walls.

##  I Introduction

Ferroelectric and multiferroic materials are becoming an increasingly
important set of materials for next-generation devices due to the ability to
manipulate non-volatile states by application of electric fields[1, 2]. The
integration of these materials into spintronic devices offers the potential to
improve existing device designs by replacing current-controlled operations
with electric field-controlled ones which could reduce the energy required to
operate devices by orders of magnitude[3, 4].

Associated with the ferroelectric order is a substantial lattice elongation,
on the order of 1% of the in-plane lattice parameter, that makes these
materials suitable for imparting a large strain. Conventional magnetic
multilayer devices, coupled to a ferroelectric material, can be modified and
controlled by applying an electric field to modify the strain, lowering the
power needed to operate a device[5]. This makes ferroelectrics a core material
in the field of straintronics[6, 7] in which devices can be modified or
manipulated by strain. The scope of these modifications is very broad, but the
applications in which ferroelectrics play a particularly large role include
thin film systems where ferroelectric materials can be used to control the
interfacial strain. Indeed, research involving piezoelectric and ferroelectric
materials has already shown that it is possible to control magnetic domain
wall velocities[8], domain structure[9], and the magnitude of micromagnetic
parameters such as effective anisotropy[10].

In our previous work[11], we have shown that growing a ferromagnetic film of
Co40Fe40B20 on a (111)-oriented BaTiO3 substrate couples the magnetic domain
structure to the ferroelectric domain structure and leads to two in-plane
magnetoelastic anisotropy configurations (illustrated in Fig. 1) in which the
imprinted magnetoelastic anisotropy rotates by either 60∘ or 120∘ between
adjacent ferroelectric domains, and that this impacts the magnetic domain
structure and the magnetic field response of magnetic domain walls, giving
rise to the 60U, 60C, 120U and 120C domain wall states. Previous
investigations of this system were limited to static effects at room
temperature and without applied electric field. In this paper we will detail
the effects of temperature.

BaTiO3 (BTO) exhibits three crystal phase transitions from rhombohedral to
orthorhombic at 190 K, orthorhombic to tetragonal at 280 K, and tetragonal to
cubic at 400 K[12] (dependent on crystal quality and if the sample is being
heated or cooled), with the rhombohedral (R), orthorhombic (O), and tetragonal
(T) crystal phases all displaying ferroelectric order. Each of these phases
has different preferred orientations of ferroelectric polarization and so
these phase transitions are accompanied by large changes in the lattice
parameters of BTO, particularly along the directions of ferroelectric
polarization where the lattice is elongated. Interestingly, the projection of
polarization onto the (111) plane in all phases results in the same two
possible rotations of in-plane projection of the polarization, where the
polarization rotates by either 60∘ or 120∘ through the domain wall. This angle
directly affects the magnetic domain wall angle, reducing the rotation of
magnetization between adjacent domains and creating the possibility for the
domain walls to be either head-to-head (charged) or head-to-tail (uncharged).

Here, we investigate the temperature response of coupled BTO(111)/CoFeB films
both in bulk and at the individual domain level. By sweeping the temperature
we can access different polar phases with large changes in lattice elongation
and assess the viability of strain-dependent devices based on BTO(111)
substrates, measured by the response of the properties of the ferromagnetic
film deposited on top. Similar to the work performed in the literature[13], we
find that the magnetic thin film remains coupled in all crystal phases
regardless of the change in strain.We also determine values for the
magnetoelastic anisotropy in CoFeB within individual coupled ferroelectric
domains across the entire temperature-dependent structural range for BTO, from
77 K (below which we expect no significant changes in the lattice constants)
to 417 K (above which there is no ferroelectric order). Finally, we perform
micromagnetic simulations of the domain wall width to show the impact these
changes represent and find that at the very extremes of the temperature range
the difference between domain wall configurations becomes significant for
charged domain wall structures. This will have important implications for
devices making use of them such as domain wall resonators[14, 15].

Figure 1: Schematic illustration of the magnetization and anisotropy
configurations in the strain-coupled ferromagnet. The domain wall
configurations are labelled by the underlying angle between the magnetoelastic
anisotropy axes (60∘ or 120∘) and the charged or uncharged nature of the
domain wall structure leading to four configurations: 60U, 60C, 120U and 120C.
Black double-headed arrows indicate the direction of magnetoelastic
anisotropy, blue arrows represent the direction of local magnetization.

##  II Methods

Samples were grown by DC sputtering onto commercially available BaTiO3(111)
substrates purchased from SurfaceNet. These substrates display ferroelectric
domains as received and preparation involved cleaning by sonication in acetone
and then isopropanol. Using the Royce deposition system[16] at Leeds, a
ferromagnetic Co40Fe40B20 layer is deposited at 300∘C above the Curie
temperature of the substrate to promote more efficient strain transfer. The
full sample structure is BaTiO3(111)/Co40Fe40B20(20nm) @ 300∘C/Pt(5nm) @ 27∘C.
As in our previous work, the as-received domain structure of the BTO is
diverse and contains regions in which the polarization rotates through either
60∘ or 120∘ between adjacent ferroelectric domains. This results in a
diversity of the accompanying magnetic configurations once the ferromagnetic
layer is deposited. More details can be found in our previous work[11].

Magnetometry was performed using a Quantum Design MPMS3 superconducting
quantum interference device-vibrating sample magnetometry (SQUID-VSM) to
investigate the bulk changes in magnetic moment and magnetization over a
temperature range of 77-300 K in magnetization v temperature (MvT)
measurements. Values of saturation magnetization were extracted from $\pm$ 1 T
magnetic hysteresis loops performed at the respective temperatures with
background contributions subtracted. A temperature-dependent Bloch
$T^{\frac{3}{2}}$ fit[17] to the data is performed from the values of
$M_{sat}$ obtained and used to extrapolate values of $M_{sat}$ between these
data points and above room temperature up to 417K. MvT measurements are
presented with no background contributions removed.

Local hysteresis measurements below room temperature were carried out using a
wide-field Kerr microscope[18] with an optical cryostat at the EXTREMAG[19]
facility in Exeter with a 60$\times$ objective lens. The sample is cooled to a
base temperature of 77 K and then warmed to achieve the desired temperature.
Magnetic hysteresis loops at each temperature are taken along the magnetic
hard axis of an individual region coupled to a ferroelectric domain and the
anisotropy field is calculated from the resulting loop. Measurements above
room temperature were performed at Leeds with a heater stage in the same
manner. All measurements are taken in the positive (warming) temperature
direction.

##  III Results

We first study the change in magnetic moment with temperature using SQUID-VSM
magnetometry. The temperature response of the sample is measured in an applied
field of 200 mT between room temperature and 77 K. This temperature is well
below the rhombohedral phase transition at 190 K. A magnetic field of this
strength, greater than the local anisotropy field measured previously at room
temperature, is chosen to bias the moment towards the magnetic field direction
across phase transitions. The results in Fig. 2 show two hysteretic changes in
magnetic moment around 190 K and 280 K, corresponding to the rhombohedral-
orthorhombic and orthorhombic-tetragonal phase transitions of the BTO
respectively. This hysteretic behaviour is an expected property of the BTO
substrates, which have been widely reported[12, 20, 21] to have hysteretic
phase transitions in a variety of properties including the lattice parameters,
which are closely coupled to the lattice elongations of the ferroelectric
order. Similar behaviour has been observed in epitaxial BTO-ferromagnet
systems in the literature and attributed to either the change in strain[20] or
the change in ferroelectric domain structure (and so, the resulting magnetic
easy axes) [21] and demonstrates strong strain coupling between the substrate
and film in all crystal phases.

Next we investigate the change in saturation magnetization with temperature.
In Fig. 3 we show the saturation magnetization as a function of temperature.
The result of this is values that fit well to a standard $T^{3/2}$ fit with no
discontinuities in $M_{s}(T)$ indicating that the strain is insufficient to
change the saturation magnetization. This means that the temperature
hysteresis observed in Fig. 2 most likely corresponds to the abrupt change in
ferroelectric domain structure at the phase transitions of the BTO substrate.
While 200mT is in excess of the coercive field at all temperatures, it is not
sufficient to pin the magnetization to the field direction as the
ferroelectric domain structure, and consequently the local magnetoelastic
anisotropy, abruptly changes. Performing the same MvT with a larger 1T field
suppresses the change in moment across these transitions (see Supplementary
Fig. 1). In the case of these amorphous thin-film ferromagnets, it would seem
that the change in domain structure is more important than the change in
strain for these whole-sample measurements.

Figure 2: MvT measurement in an applied field of 200mT. The blue curve
corresponds to data taken during cooling, with the red curve taken while
heating to room temperature. The regions corresponding to the rhombohedral
(R), orthorhombic (O), and tetragonal (T) phases of BaTiO3 are indicated.
Figure 3: Saturation magnetization (points) for 20nm CoFeB thin film coupled
to a BTO(111) substrate, extracted from hysteresis loops at each temperature
and the corresponding Bloch $T^{3/2}$ law fit (dashed line). Indicated are the
temperatures at which the rhombohedral (R), orthorhombic (O) and tetragonal
(T) phase transitions occur. No significant change around these points is
observed. Figure 4: Ferromagnetic domain images taken in the same vicinity for
the a) rhombohedral (R), b) orthorhombic (O), and c) tetragonal phases (T).
The magnetic stripe orientation changes to match the ferroelectric domain
structure demonstrating good coupling between films across phase transitions.

The results from the SQUID-VSM measurements in Figs. 2 and 3 represent the
volume-averaged response of the sample. As mentioned previously there are many
ferroelectric domains within the sample with different orientations of lattice
elongation that leads to a variety imprinted magnetoelastic easy axes in the
CoFeB film. While we are able to infer strain coupling between the
ferroelectric substrate and ferromagnetic film, it is not possible to probe in
detail the effect of the substrate on the ferromagnetic film as we do not have
a monodomain substrate, making it more challenging to interpret the results.
To understand in more detail the effect of temperature on the strength of the
magnetoelastic anisotropy within one ferroelectric domain we performed local
measurements using a Kerr microscope to focus our measurements down to one
region of ferroelectric domain structure and coupled regions corresponding to
individual ferroelectric domains.

Using optical cryostat attachments we perform wide-field Kerr microscopy on
the sample in two setups, a cryo-stage with liquid nitrogen that allows us to
cool to 77 K (well below the orthorhombic and rhombohedral phase transitions)
and a heater-stage in which we heat to above 420K (the tetragonal to cubic
phase transitions). Fig. 4 shows the result of domain imaging in the same
region across three different phases, imaged here in zero applied magnetic
field. The presence of stripe domains in all three polar phases are in
excellent agreement with the SQUID data presented previously and work carried
out on (100)-oriented substrates in the literature[13], and shows that the
strain is sufficient to couple the domains at all temperatures.

In this experimental setup, it is not possible to precisely ascertain the
ferroelectric (and anisotropic) configuration of the region before and after
the phase transition as the sample cannot be freely rotated so a rigorous
extrapolation of the easy axes cannot be performed. However, we can gain some
understanding of the configuration from the orientation of the domain walls
and the angle at which 180∘ domain walls are canted relative to the stripe
axis. In the tetragonal and rhombohedral phases, the projection of lattice
elongation onto the (111) surface lies in [11$\bar{2}$]-like directions and so
we expect that stripes oriented in the same direction in both crystal phases
will maintain the same ferroelectric configuration with the rotation between
the easy axes remaining the same[11]. In the orthorhombic phase, the lattice
elongation lies along [01$\bar{1}$]-like directions, which is the orientation
of domain walls in the tetragonal and rhomobohedral phases. As a result, a
rotation of the domain wall by 30∘ or 90∘ indicates that the ferroelectric
configuration remains the same, and a rotation of 0∘ or 60∘ is indicative of a
change in the ferroelectric configuration from either a 60∘ to 120∘ or vice
versa.

In the domain images presented here, the 180∘ reversal domains indicate that
the T-phase region has a 60∘ configuration, and that the O-phase image has a
120∘ configuration, and indeed the stripe orientation rotates by 60∘ as we
expect it should for a change in ferroelectric domain type. The stripe
orientation is in the same crystal direction in the R and T phases from which
we infer that the ferroelectric configuration is the same in both phases. We
have shown the case where the ferroelectric configuration remains the same
between the R and T phases - we have however also observed a rotation of 90∘
suggesting that both outcomes are equally likely. Separately, in the heater
stage we observe that the magnetic stripe domain structure vanishes above the
Curie temperature of BTO ($\approx$420K).

Starting from a base temperature of 77 K the value of the anisotropy field is
extracted from hard-axis hysteresis loops measured from individual stripe
domains and the value of anisotropy calculated from
$H_{k}=\frac{2K_{eff}}{\mu_{0}M_{s}}$. Values for $M_{s}$ are used from our
previously presented data where possible and the obtained fit is used to
extrapolate values for high temperatures.

Figure 5: Local magnetoelastic anisotropy imprinted in the CoFeB film by the
BaTiO3(111) substrate. Dashed lines indicate the phase transitions between the
rhombohedral (R), orthorhombic (O) and tetragonal (T) phases of the substrate.
Measurements in blue ($\filledmedtriangleup$) are taken at low temperature
with the use of an optical cryostat, and measurements in red
($\filledmedtriangledown$) were taken separately in a heater stage.

The results are summarized in Fig. 5. The change in anisotropy far from the
phase transitions shows a step-like change for the measurements performed in
the cryo-stage, with different transition behaviour. From the rhombohedral to
orthorhombic phase there is a change in the magnitude of the magnetoelastic
anisotropy from an average value of (14$\pm$1) kJ/m3 in the rhombohedral phase
to a value of (20$\pm$1) kJ/m3 in the orthorhombic phase. In the region around
the O-R phase transition the magnitude of the magnetoelastic anisotropy
smoothly transitions between these two average values. Whereas, when going
from the orthorhombic phase to the tetragonal the increase is sharp with
presumably some saturation to a maximum value around room temperature. To
examine this in more detail we calculate in Fig. 6 the expected magnitude of
the lattice elongations (relative to the in-plane lattice constant) in the
(111) plane based on the reported results from Kwei et al[22].

In the rhombohedral phase there is minimal change in the lattice parameters
and indeed we observe no significant changes to the measured magnetoelastic
anisotropy far from the orthorhombic phase transition. In the orthorhombic
phase, we expect two distinct branches corresponding to directions of
ferroelectric polarization with either the major or minor component of the
biaxial lattice elongation lying in the (111) plane. We expect that the
shorter lattice elongation should show a continuous change across the phase
transition, however in our results the same behaviour is instead reflected in
what must correspond to the larger lattice elongation as the magnetoelastic
anisotropy increases from the rhombohedral to orthorhombic phase. The
discontinuous ‘jump’ in the lattice elongation at the orthorhombic-tetragonal
transition is well reflected in our measurements of $K_{me}$.

Figure 6: Calculated magnitude of the lattice elongations on the (111) plane
using structural data obtained by Kwei et al[22], expressed as a percentage of
the in-plane lattice parameter.

For measurements performed on the heater stage, there is a large discrepancy
that leads to a discontinuous transition between the two datasets. We
attribute this to effects observed in previous work[11], where there can be a
large difference in the magnitude of Kme between the 60∘ or 120∘ domain
states, and measurements here were performed on different regions of the
sample. Regardless, it shows the changes from the initial state that we expect
with an approximately linear reduction from 315-385 K corresponding to the
decrease in lattice elongations and then a sharp drop off close to the
ferroelectric Curie temperature at approximately 420 K (in agreement with the
Curie temperature seen in the work by Lahtinen et al[13]) where the polar
order reduces. Above this Curie temperature, we do not observe any magnetic
stripe domains as the ferroelectric order is no longer present and there is no
meaningful in-plane anisotropy.

Using our measured values of $K_{me}$ we explore the resulting implications
that the change in magnetoelastic anisotropy has for the magnetic domain wall
width (DWW). In a micromagnetic[23] framework we have studied the effect of
temperature on the DWW of charged and uncharged domain wall configurations.
The simulation is divided into three distinct regions of uniaxial anisotropy
corresponding to different directions of lattice elongation in adjacent
ferroelectric domains. This is illustrated in Supplementary Fig. 2. The
central stripe width is set to be 2 $\mu$m and the micromagnetic cell size is
2 nm x 2 nm x 20 nm, and periodic boundary conditions are used in the $x$ and
$y$ directions. Values of $K_{u1}$, representing the parameter previously
defined as $K_{me}$, and $M_{sat}$ are informed from the above experiment with
the assumption made that similar magnitudes of $K_{me}$ can be expected for
all ferroelectric domains. Values of $A_{ex}$ are taken from previous work
examining the temperature dependence of exchange stiffness[24] of Co which is
likely to be an overestimate for a CoFeB system. This is intended not to
predict the exact values of the DWW but instead to examine how it scales with
temperature taking into account the relative scaling of all relevant
micromagnetic parameters.

For all phases of BTO the projection of the polarization onto the (111) plane
leads to an angle between adjacent lattice elongations of either 60∘ or 120∘,
depending on which axes the polarization switches between. In each
ferroelectric domain, the magnetization is strongly pinned to the imprinted
magnetoelastic anisotropy axis which leads to magnetic domain walls that are
pinned to the ferroelectric domain walls with a reduced wall angle dependent
upon how the bulk magnetization rotates between the ferroelectric domains. The
magnetization can rotate in a head-to-tail or head-to-head fashion leading to
charged or uncharged domain wall structures. In total, this means that there
are four domain wall structures to consider by combining the possible
ferroelectric rotations and charged or uncharged states: a rotation of 60∘
with charged (60C) or uncharged (60U) character and a rotation of 120∘ with
charged (C) or uncharged (U) character. These were previously illustrated in
Fig. 1.

Figure 7: Domain wall width in the charged (C) and uncharged (U)
configurations as a function of temperature. Numbers in the legend represent
the total angle through the domain wall between the magnetoelastic anisotropy
axes in adjacent stripe domains corresponding to the rotation of polarization.

Within each of the pre-defined stripes, the magnitude of the in-plane
anisotropy remains the same but the orientation of the magnetoelastic
anisotropy is varied (rotating by either 60∘ or 120∘). Periodic boundary
conditions are used in both the $x$ and $y$ dimensions and the domain wall is
homogenous in the $y$-axis. To calculate the domain wall width, $\delta$, we
use the method outlined in previous work[11, 25] using the integral definition
of the domain wall width:

$\delta=\int_{-\infty}^{\infty}\cos^{2}(\phi^{\prime})\,dx,$ (1)

with $\phi^{\prime}$ being the reduced magnetization angle,

$\phi^{\prime}=\bigg{(}\phi-\frac{\mid\phi_{\frac{d}{2}}-\phi_{\frac{-d}{2}}\mid}{2}\bigg{)}\frac{180}{\mid{\phi_{\frac{d}{2}}-\phi_{-\frac{d}{2}}\mid}},$
(2)

with $\phi$ being the magnetization angle measured relative to the direction
of the easy axis and $\mid\phi_{\frac{d}{2}}-\phi_{\frac{-d}{2}}\mid$ being
the change in magnetization angle between magnetic domains measured far from
the domain wall.

The variation of each of these domain wall widths with temperature is shown in
Fig. 7. In all cases, the scaling is most strongly dominated by the change in
$K_{me}$ which has an inversely proportional change in the magnetic domain
wall width. The result is a profile that shows the opposite temperature
scaling to that shown in Fig. 5, with the very extremes of the measured
temperatures (77 K, 417 K) displaying the largest values of domain wall width.
The charged domain wall structures show the largest absolute change in DWW,
representing something that would in principle be easier to measure using a
high-resolution technique such as photo-emission electron microscopy to probe
the domain wall profile.

Differences resulting from the different magnetoelastic configurations become
apparent only at the temperature extremes with changes in all uncharged
configurations being extremely small, on the order of a few nanometers of
difference, while the difference between the 60C and 120C is on the order of
50nm at the largest. The differences at these extremes stem only from the
change in angle between the easy axes, indicating that the charged domain wall
width could be a useful way of measuring the two distinct states.

##  IV Conclusion

In summary, we have studied the temperature dependence of BaTiO3(111)/CoFeB
heterostructures using SQUID-VSM and Kerr microscopy and determined the change
in magnetoelastic anisotropy as a result of temperature-dependent strain
transfer. We find strong coupling between the substrate and thin film at all
temperatures, with the local anisotropy varying as the lattice parameters of
the BaTiO3 change with temperature. The absolute changes in the magnetoelastic
anisotropy either side of room temperature are approximately equal, and so
this presents two routes by which its magnitude could be changed in a device.
Following this, we performed simulations informed by our experiments to
determine how this dependence dominates the domain wall width of magnetic
domains. These domain wall widths were found to be most strongly dependent on
the magnetoelastic anisotropy, with the low-anisotropy charged domain walls
showing a large difference in domain wall widths between the 60∘ and 120∘ \- a
result purely of the difference in anisotropy configuration and resulting wall
angles of the magnetic domain wall structure.

These measurements demonstrate two distinct regions of interest around the
rhombohedral-orthorhombic and orthorhombic-tetragonal transitions with the
low-temperature results showing step-wise changes in the magnitude of the
anisotropy, and high-temperature results demonstrating a drop off as the polar
ordering weakens. Both routes to manipulating the magnetic anisotropy lead to
similar absolute changes in the magnitude of the anisotropy. Control of
magnetoelastic anisotropy is fundamental to devices based on multiferroic
heterostructures. The anisotropy tuning that we show here will be useful for
future devices based on BaTiO3(111).

###### Acknowledgements.

We acknowledge support from the Henry Royce Institute and funding from the
Engineering and Physical Sciences Research Council (EPSRC) Grant No.
EP/M000923/1. Variable-temperature, wide-field Kerr microscopy was performed
at the Exeter Time-Resolved Magnetism Facility (EXTREMAG) (EPSRC Grant
Reference EP/R008809/1 and EP/V054112/1). This project received funding in
part from the European Union’s Horizon 2020 research and innovation programme
under the Marie Sklodowska-Curie grant agreement No 750147. R. G. H.
acknowledges the support of an EPSRC DTA studentship.

## Author Declarations

## Conflict of interest

The authors have no conflicts to disclose.

## Author Contributions

R. G. Hunt: Investigation (lead); methodology (lead); conceptualization
(equal); writing - original draft (lead); visualization (lead); formal
analysis (lead); writing - review and editing (equal). K. J. A. Franke:
Conceptualization (equal); formal analysis (supporting); writing - review and
editing (equal). P. S. Keatley: Investigation (supporting), Methodology
(supporting), writing - review and editing (equal). P. M. Shepley:
Investigation (supporting); methodology (supporting), writing - review and
editing (equal). M. Rogers: Investigation (supporting); methodology
(supporting); writing - review and editing (equal). T. A. Moore: Project
administration (lead); funding acquisition (lead); conceptualization (equal);
methodology (supporting); investigation (supporting); writing - review and
editing (equal).

## Data Availability Statement

The data associated with this paper is available from the Leeds research
depository at https://doi.org/10.5518/1349.

## References

  * Béa  _et al._ [2008] H. Béa, M. Gajek, M. Bibes, and A. Barthélémy, “Spintronics with multiferroics,” Journal of Physics: Condensed Matter 20, 434221 (2008).
  * Ajayan  _et al._ [2023] J. Ajayan, P. Mohankumar, D. Nirmal, L. L. Joseph, S. Bhattacharya, S. Sreejith, S. Kollem, S. Rebelli, S. Tayal, and B. Mounika, “Ferroelectric field effect transistors (FeFETs): Advancements, challenges and exciting prospects for next generation non-volatile memory (NVM) applications,” Materials Today Communications 35, 105591 (2023).
  * Franke  _et al._ [2015] K. J. Franke, B. Van de Wiele, Y. Shirahata, S. J. Hämälainen, T. Taniyama, and S. van Dijken, “Reversible electric-field-driven magnetic domain-wall motion,” Physical Review X 5, 1–9 (2015), arXiv:1411.6798 .
  * Yin and Li [2017] Y. Yin and Q. Li, “A review on all-perovskite multiferroic tunnel junctions,” Journal of Materiomics 3, 245–254 (2017).
  * Damjanovic [1998] D. Damjanovic, “Ferroelectric, dielectric and piezoelectric properties of ferroelectric thin films and ceramics,” Reports on Progress in Physics 61, 1267–1324 (1998).
  * Roy, Bandyopadhyay, and Atulasimha [2011] K. Roy, S. Bandyopadhyay, and J. Atulasimha, “Hybrid spintronics and straintronics: A magnetic technology for ultra low energy computing and signal processing,” Applied Physics Letters 99, 063108 (2011).
  * Bukharaev  _et al._ [2018] A. A. Bukharaev, A. K. Zvezdin, A. P. Pyatakov, and Y. K. Fetisov, “Straintronics: a new trend in micro- and nanoelectronics and materials science,” Physics-Uspekhi 61, 1175–1212 (2018).
  * Shuai  _et al._ [2022] J. Shuai, M. Ali, L. Lopez-Diaz, J. E. Cunningham, and T. A. Moore, “Local anisotropy control of Pt/Co/Ir thin film with perpendicular magnetic anisotropy by surface acoustic waves,” Applied Physics Letters 120, 252402 (2022).
  * Lahtinen, Franke, and Van Dijken [2012] T. H. Lahtinen, K. J. Franke, and S. Van Dijken, “Electric-field control of magnetic domain wall motion and local magnetization reversal,” Scientific Reports 2, 1–6 (2012), arXiv:1109.5514 .
  * Shepley  _et al._ [2015] P. M. Shepley, A. W. Rushforth, M. Wang, G. Burnell, and T. A. Moore, “Modification of perpendicular magnetic anisotropy and domain wall velocity in Pt/Co/Pt by voltage-induced strain,” Scientific Reports 5 (2015), 10.1038/srep07921.
  * Hunt  _et al._ [2023] R. G. Hunt, K. J. A. Franke, P. M. Shepley, and T. A. Moore, “Strain-coupled domains in BaTiO3(111)-CoFeB heterostructures,” Phys. Rev. B 107, 014409 (2023).
  * Acosta  _et al._ [2017] M. Acosta, N. Novak, V. Rojas, S. Patel, R. Vaish, J. Koruza, G. A. Rossetti, and J. Rödel, “BaTiO3-based piezoelectrics: Fundamentals, current status, and perspectives,” Applied Physics Reviews 4 (2017), 10.1063/1.4990046.
  * Lahtinen and van Dijken [2013] T. H. E. Lahtinen and S. van Dijken, “Temperature control of local magnetic anisotropy in multiferroic CoFe/BaTiO3,” Applied Physics Letters 102, 112406 (2013).
  * Hämäläinen  _et al._ [2018] S. J. Hämäläinen, M. Madami, H. Qin, G. Gubbiotti, and S. van Dijken, “Control of spin-wave transmission by a programmable domain wall,” Nature Communications 9, 1–8 (2018).
  * de Wiele  _et al._ [2016] B. V. de Wiele, S. J. Hämäläinen, P. Baláž, F. Montoncello, and S. van Dijken, “Tunable short-wavelength spin wave excitation from pinned magnetic domain walls,” Scientific Reports 6 (2016), 10.1038/srep21330.
  * [16] The Royce Deposition System is a multi-chamber, multi-technique thin film deposition tool based at the University of Leeds as part of the Henry Royce Institute.
  * Argyle, Charap, and Pugh [1963] B. E. Argyle, S. H. Charap, and E. W. Pugh, “Deviations from ${T}^{\frac{3}{2}}$ law for magnetization of ferrometals: Ni, Fe, and Fe+3% Si,” Phys. Rev. 132, 2051–2062 (1963).
  * Soldatov and Schäfer [2017] I. V. Soldatov and R. Schäfer, “Selective sensitivity in Kerr microscopy,” Review of Scientific Instruments 88 (2017), 10.1063/1.4991820, arXiv:1612.02027 .
  * [19] The Exeter time-resolved magnetism (EXTREMAG) facility is an EPSRC funded ultrafast laser facility that offers its users low temperature and high magnetic field environments in a variety of time-domain optical pump-probe, microscopy, and terahertz experiments.
  * Vaz  _et al._ [2009] C. A. Vaz, J. Hoffman, A. B. Posadas, and C. H. Ahn, “Magnetic anisotropy modulation of magnetite in Fe3O4/ BaTiO3 (100) epitaxial structures,” Applied Physics Letters 94, 92–95 (2009).
  * Venkataiah  _et al._ [2012] G. Venkataiah, Y. Shirahata, I. Suzuki, M. Itoh, and T. Taniyama, “Strain-induced reversible and irreversible magnetization switching in Fe/BaTiO3 heterostructures,” Journal of Applied Physics 111, 033921 (2012).
  * Kwei  _et al._ [1993] G. H. Kwei, A. C. Lawson, S. J. Billinge, and S. W. Cheong, “Structures of the ferroelectric phases of barium titanate,” Journal of Physical Chemistry 97, 2368–2377 (1993).
  * Vansteenkiste  _et al._ [2014] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, “The design and verification of MuMax3,” AIP Advances 4 (2014), 10.1063/1.4899186, arXiv:1406.7635 .
  * Moreno  _et al._ [2016] R. Moreno, R. F. L. Evans, S. Khmelevskyi, M. C. Muñoz, R. W. Chantrell, and O. Chubykalo-Fesenko, “Temperature-dependent exchange stiffness and domain wall width in Co,” Phys. Rev. B 94, 104433 (2016).
  * Franke  _et al._ [2014] K. J. Franke, D. López González, S. J. Hämäläinen, and S. Van Dijken, “Size dependence of domain pattern transfer in multiferroic heterostructures,” Physical Review Letters 112, 1–5 (2014), arXiv:1307.0605 .


# Stability of a tilted granular monolayer: How many spheres can we pick
before the collapse?

Eduardo Rojas<EMAIL_ADDRESS>Department of Mechanical Engineering,
Universidad de Antofagasta. Antofagasta, Chile.  Héctor Alarcón  Instituto de
Ciencias de la Ingeniería, Universidad de O’Higgins. Rancagua, Chile.  Vicente
Salinas  Instituto de Ciencias Químicas Aplicadas, Facultad de Ingeniería,
Universidad Autónoma de Chile. Santiago, Chile.  Gustavo Castillo  Instituto
de Ciencias de la Ingeniería, Universidad de O’Higgins. Rancagua, Chile.
Pablo Gutiérrez<EMAIL_ADDRESS>Instituto de Ciencias de la
Ingeniería, Universidad de O’Higgins. Rancagua, Chile.

###### Abstract

A new way of triggering avalanches is investigated using discrete element
simulations for a process of random extraction of spheres. A monolayer, formed
by identical spheres in a hexagonal configuration, is placed on a tilted plane
surrounded by a small fence that sustains the spheres, mimicking the disposal
of fruits in the market. Then, a random continuous extraction process of
spheres is imposed until the collapse. For this simple numerical experiment, a
phase diagram was obtained to visualize the occurrence of avalanches triggered
by vacancies as a function of the tilting angle and the size of the system.
More importantly, a sub-zone was found where we can predict the critical
number of extractions until the avalanche takes place. The prediction is made
from an evolution model of the average coordination number and elastic energy,
adding an empirical relationship between initial and critical average elastic
energy. The theoretical prediction is valid for a limited range of angles,
where avalanches occur at a critical packing fraction.

##  I Introduction

Granular avalanches are large-scale phenomena that can be triggered at small
scales [1, 2, 3]. However, the effectiveness of a small-scale precursor (for
example, a particle hitting another after falling a small distance) depends on
its ability to produce a chain effect that collectively breaks the strength of
the entire system. This duality between small-scale precursors and their
large-scale effect is still an open question. Some progress has been made in
revealing the small-scale phenomenology by imposing shear through the system’s
walls [4, 5, 6, 7]. This forcing corresponds to a large-scale perturbation
imposed instantaneously to the whole system, as it is also when inclining,
rotating, or shaking the entire granular assembly. Motivated by everyday
phenomena, we want to see how random-small-scale perturbations progressively
cause the loss of global stability, triggering large-scale avalanches.

The example that inspired us is depicted in Fig. 1: vendors in the market
stack the fruit trying to make most of it accessible to whoever wants to take
it, which happens by placing the container with an optimal inclination. If the
inclination is too slight, some fruit will be more difficult to reach. For a
higher inclination, it will be easier to pick the fruit, but the stability of
the stack could be compromised. Especially for this last case, one may ask how
many fruits one can extract without producing the collapse of the pile (see
Fig. 1). This simple system, familiar to everyone, presents two novelties
compared to the literature on granular materials. First, fruits are generally
displayed in hexagonal or square arrangements, which is possible because of
their low size-dispersion. However, this perfectly crystallized configuration
remains poorly studied. The second remarkable novelty corresponds to the way
of triggering the avalanche. The generation of vacancies by fruit extractions
corresponds to a new precursor of avalanches, producing both local
rearrangements and a progressive decline in the coordination number, which
weakens the contact structure globally. Intuitively, such a system should be
ruled by the contact structure between neighbors, determining forces
distribution and system’s rigidity, and, more broadly, by the system’s packing
fraction.

Figure 1:  Avalanches triggered by random extractions. (a) Arrangement of
fruits in the market (photo credits in [8]). The bottom panels show a
numerical simulation for a system with 90 spheres, randomly colored. (b)
Presents a global view of the system, tilted with an angle $\theta$. In (c),
the first sphere is extracted and the second in (d). In (e), an avalanche
takes place after the extraction of 12 spheres. The inset in (c) highlights
the fence bounding the system.

A disordered monolayer of particles placed on a tilted plane was considered
before by Dorbolo [9]: the slope of the plane is increased quasi-statically
until a critical angle, where the stack becomes unstable and produces an
avalanche. In some cases, only the lower part of the particles falls, while
the rest remains in place. This observation represents a manifestation of the
Jansen effect, revealing the complex network of contacts between spheres [9].
Contact networks in granular assemblies indeed show filamentary structures or
force chains ruling the stability of the media [10, 11, 12, 13] and physically
connecting micro and macro scales. For instance, a local failure could be
extended in a long-ranged response through a force chain without producing
significant changes in particles’ positions or orientations [14, 6].
Therefore, exploring the so-called mesoscale variables in granular avalanches
is essential. Accordingly, we analyze the coordination number and the elastic
energy through contacts, which is related to the intensity of force chains.

Finally, it should be noticed that in compaction processes, particularly in
monolayers, an increase in the packing fraction is usually observed [15, 16,
17]. However, suppose the compaction process triggers a flowing state. In that
case, it generally does through a slight decrease of the packing fraction just
before the avalanche, related to Reynolds dilatancy [18]. Indeed, Levy dit
Vehel and collaborators [19] subjected a monolayer to shear by applying torque
to a cylindrical film configuration, stressing the relationship between the
magnitude of the global dilation and the intensity of catastrophic events.
Consequently, it seems reasonable to examine dilation to forecast sudden
events like avalanches.

This article is organized as follows. Our model is presented in section II,
together with its numerical implementation. Our results are presented in
section III, particularly a phase diagram for avalanches and the study of the
evolution of the coordination number and elastic energy in the system. Then,
the regime where avalanches are ruled by the decrease of elastic energy is
discussed, followed by the conclusions in Section IV.

##  II  Setup and numerical implementation

Figure 2:  Phase diagram for avalanches triggered by random extractions. (a)
Migration rate $r_{m}$ of spheres as a function of the inclination angle.
Symbols represent different system sizes. The horizontal dashed line at 0.8
shows our threshold to define an avalanche. (b) Phase diagram of the three
main regimes observed, as a function of the angle $\theta$ and system size
$f$: No avalanche (crosses); Instant avalanche (triangles) and Avalanche
(circles). The shadowed zone corresponds to a region where avalanches can be
predicted (latter explained in the main text).

The present study was performed on a numerical platform using Molecular
Dynamics [20], which is one of the discrete element methods (DEM). This
technique allowed us to compute spheres’ movements and interactions, giving
detailed information about their contacts.

The setup consists of a rectangular monolayer arrangement of identical
spherical particles of diameter $d=$4\text{\,}\mathrm{cm}$$ and mass
$m=$2.93\text{\,}\mathrm{g}$$, similar to table-tennis balls. These spheres,
subjected to gravity in the vertical direction, are disposed on a tilted
plane, with an angle $\theta$ with respect to the horizontal plane (see Fig.
1.b). Initially, the spheres are crystallized forming a perfect hexagonal
configuration. A small square cross-section fence sustains the arrangement.
This barrier is $0.17d$ height and is located all over the monolayer boundary
(see inset in Fig. 1.c). The smallest system we used corresponds to a
monolayer of $90$ spheres. Due to the hexagonal packing, horizontal rows
alternate between eight and seven spheres in this case (see Fig. 1.b). To
analyze size effects, we scaled the system by a factor $f$ from 1 to 5 and
applied it to the smallest rectangle’s width and length. To initialize the
simulations, after the spheres are placed on the tilted surface, a gentle
horizontal movement is applied to the whole system to settle the spheres in a
stable position. The velocity $v(t)=v_{0}\sin(2\pi t/T)$ is applied to the
system for $2\text{\,}\mathrm{s}$, where
$v_{0}=$0.1\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$ and
$T=$1\text{\,}\mathrm{s}$$. During this initial process, the relative movement
between spheres is negligible and the hexagonal configuration is preserved.

After initialization, a random extraction process starts, taking the spheres
one at a time from inside the rectangle that contains and supports the spheres
111For each extraction process, all the spheres in the system are labeled, and
one of them is chosen according to a uniform random distribution.. A velocity
perpendicular to the inclined plane is imposed to each particle selected to be
removed (see Figs. 1.c and 1.d). The interval between extractions was three
seconds, enough time to allow spheres to relocate to their new stable position
or escape outside the fence. By doing so, the quasi-static evolution of the
system is better described by the number of spheres extracted rather than
time.

The discrete element simulations are carried out using the ViscElMat module of
the YADE platform [22]. The Molecular Dynamics method makes use of particle-
particle and particle-plane contact forces when there is an overlap $\delta$
between the elements. The model for the normal contact force $F_{n}$ is
viscoelastic in its linear form: $F_{n}=k_{n}\delta+c_{n}\dot{\delta}$, where
$k_{n}$ and $c_{n}$ are constants. The tangential contact force $F_{t}$ has
the same form with constants $k_{t}$ and $c_{t}$, added to Coulomb friction
with the same static and dynamic friction coefficient $\mu$. The constant
$k_{n}$ for the normal force is set to
$10.000\text{\,}\mathrm{N}\mathrm{/}\mathrm{m}$ to obtain an overlap
$\delta<10^{-3}d$ in all the simulations performed, condition to remain in the
hard particle limit. The first constant in the tangential force is fixed to
$k_{t}=0.5k_{n}$, as used in numerous DEM works (see for instance [23]). The
dissipative terms are fixed to obtain a restitution coefficient $e=0.7$, while
the Coulomb friction coefficient is set to $\tan($$)$.

##  III  Results

###  III.1  Phase diagram for avalanches

To decide whether an avalanche occurs or not, we consider two quantities: 1)
the number of spheres inside the system $N_{in}$; and 2) the spheres that
migrate outside the system $N_{out}$ immediately after a given extraction.
Both quantities are a function of the extraction number $n$, which represents
the umpteenth sphere extracted from the system. The maximum value of
$N_{out}$, which occurs at $n=n_{max}$, is the parameter selected to define
the migration rate $r_{m}$ as:

$r_{m}=\frac{N_{out}(n_{max})}{N_{in}(n_{max})}.$ (1)

Therefore, $r_{m}$ is a measurement of the size of the biggest departure of
spheres compared to the size of the system at that moment. We consider that
$r_{m}>0.8$ constitutes an avalanche. Figure 2a presents the rate $r_{m}$ as a
function of $\theta$ for three system sizes ($f=1,3,5$), where we draw a
dashed line for the avalanche criterion. Figure 2a shows that for a fixed
angle $\theta$, in general, bigger systems imply bigger spheres migrations
(bigger $r_{m}$ values). Following the avalanche criterion, a phase diagram
was built from data of several simulations performed as a function of $\theta$
and $f$ (see Fig. 2b). Three phases are distinguished: 1) no avalanche zone,
where $r_{m}\leq 0.8$ ($\times$ markers); 2) avalanches triggered by
extraction of spheres, where $0.8<r_{m}\leq 1$ and $n_{max}>0$ (circle
markers); and 3) instant avalanche zone, where $r_{m}=1$ and $n_{max}=0$
(triangle markers). The angle $None$ represents a constant inferior limit for
avalanches. Also, figure 2b shows that in larger systems, instant avalanches
start at smaller angles. This is simply because the boundaries play a more
significant role for smaller systems. The shadowed area in this figure
represents a sub-zone of predictable avalanches, to be discussed in section
III.4.

Figure 3:  Spatial distribution of coordination number and elastic energy. (a)
Coordination number field $z$, and (b) elastic energy field $E$. Initially,
the system contains 2370 spheres ($f$ = 5 as defined in the main text) and the
tilting angle is fixed to $\theta=$$$. Five states are presented by means of
snapshots obtained after the extraction of $n$ particles. From left to right,
are shown $n$ = 0, 60, 120, 180 and 190. In the last step, the collapse is
taking place.

###  III.2  Evolution of coordination number and elastic energy

We focus now on the coordination number $z$ and elastic energy of the contacts
$E$, obtained from simulations. The coordination number was computed by
counting lateral contacts sphere-sphere and sphere-fence (i.e.: contacts with
the bottom plane are excluded), which gives 6 for spheres in the bulk of a 2D,
densely packed hexagonal configuration. On the other hand, the elastic energy
was computed in each contact as $1/2k_{n}\delta^{2}$. We present the fields
$z$ and $E$ on the entire system (Fig. 3), and their average values (Fig. 4)
to obtain a broader picture of their evolution during the extraction process.

Figure 3a shows the coordination number field $z$ on the spheres of the system
for $f=5$ and $\theta=$$$. This figure corresponds to a sequence for different
extraction numbers ($n$ values) from the system’s initial state until it
collapses. For the initial condition ($n=0$), there is a considerable amount
of particles with less than 6 contacts. However, this loss of contacts does
not indicate a departure from the hexagonal configuration. Instead, fewer
contacts are at the boundaries, and some were lost during system
initialization. Also, for $n=0$, Fig. 3.b shows that the elastic energy $E$ is
not homogeneous but is more concentrated in a small central core, which
implies a non-hydrostatic pressure distribution. This non-hydrostatic
distribution is consistent with the Janssen effect, where the strongest
pressure is not located at the bottom of the medium but in a central place.
This effect has been reported for confined granular monolayers [9, 16]. Then,
for $n=60$ and $120$ the coordination number decreases in magnitude across the
system by forming channels of lower contacts. The elastic energy also
decreases, particularly in the central core. At $n=180$ an empty zone appears
in the upper-left corner of the system produced by a small migration of
particles. Note that this migration does not constitute an avalanche because
the avalanche criterion $r_{m}>0.8$ is not satisfied. A consequence of this
departure is an increment of the elastic energy just under the empty zone,
producing diagonal segments of strong elastic energy, which increases the
system’s strength locally. For $n=190$ the system collapses and all particles
depart outside the fence of the system (note that in Fig. 3 at $n=190$
particles are moving down). At this point, most of the contacts are lost (low
$z$ values), and the elastic energy shows a non-homogeneous behavior with
zones of strong collisions inside of a pressure-relaxed field.

Figure 4:  Evolution of averaged quantities. (a) Average coordination number
of the system $\bar{z}$ as a function of the extraction number $n$, for
different scale factor $f$ at angle $\theta=20^{o}$. The inner figure shows
the number of coordination as a function of the extraction number normalized
by the initial particle number. Segmented lines represent the model given in
equation 6. (b) Average elastic energy $\bar{E}$. The segmented and dashed
lines represent the model given in equation 11 for $\lambda=1$ and
$\lambda=1.5$, respectively. $\lambda=1.5$ corresponds to the average fitting
value of all evolution curves.

Figure 4a corresponds to the average coordination number of the system
$\bar{z}$ as a function of the extraction number $n$, for different scale
factor $f$, and considering a fixed angle $\theta=$$$. Regardless of $f$, the
coordination number starts near $5.25$, less than 6 for a hexagonal
configuration, as discussed earlier. The variable $\bar{z}$ decay monotonously
with $n$ until a final sharp drop when the avalanche is produced. All curves
follow a trend of linear decay before the avalanche, except for two punctual
intermediate drops for $f=2$ and $3$ associated to internal rearrangements
where the hexagonal configuration is lost locally. The dashed lines correspond
to Eqn. (3), coming from a simple statistical model to be described next. The
inset in Fig. 4a shows $\bar{z}$ as a function of the vacancy rate $n/N_{0}$,
where $N_{0}$ corresponds to the initial number of particles in the system. In
this inset, all curves collapse near the model, reinforcing the idea of a
statistical rule dominating the evolution of $\bar{z}$. In the end, when the
avalanche occurs, the granular medium flows and becomes less dense, which
causes a fall in all $\bar{z}$ curves. From the inset, we can notice that in
all cases, the vacancy rate is similar at collapse and near a 9%.

On the other hand, Fig. 4b shows the evolution of the average elastic energy
$\bar{E}$ obtained by summing the normal elastic energy of each contact
($1/2k_{n}\delta^{2}$), divided by the number of spheres inside the system
$N_{in}$. In this figure, it is clear that a larger size of the system ($f$
factor) implies a larger starting point of the mean elastic energy ($n=0$),
which is produced by the increment of the initial gravitational potential
energy. $\bar{E}$ also follows a linear decay trend until the collapse of the
system, where the energy drastically goes up and down. At $n=170$ for $f=5$ it
is possible to see how the average energy goes up. This effect can be
explained by the strong diagonal segments of the elastic energy formed under a
small departure of particles shown in Fig. 3 for $n=180$. The average elastic
energy evolution was also modeled by a statistical rule, which is explained in
the next point.

###  III.3  Statistical rules for $\bar{z}$ and $\bar{E}$

To formulate a simple statistical model for the average coordination number,
we will consider that the lost contacts after an extraction $n$ are
homogeneously distributed between the $N^{n}_{in}$ spheres remaining in the
system. Note that this hypothesis implies an equal and, in general, no integer
number of contacts for all particles. Denote the total number of contacts in
the system by $C_{n}$, then consider that each contact acts in two spheres:

$\displaystyle\bar{z}_{n}$ $\displaystyle=$
$\displaystyle\frac{2C_{n}}{N^{n}_{in}}.$ (2)

Assuming that the contacts in the system remain the same aside from the zones
of extraction, the contacts $C_{n}$ are equal to the contacts for the
extraction $n-1$, associated to the $N^{n}_{in}$ particles remaining in the
system (which are $\bar{z}_{n-1}N_{in}^{n}/2$), less the lost contacts around
the new vacancy ($\bar{z}_{n-1}$):

$\displaystyle\bar{z}_{n}$ $\displaystyle=$
$\displaystyle\frac{2\left({\bar{z}_{n-1}N_{in}^{n}/2-\bar{z}_{n-1}}\right)}{N^{n}_{in}}$
(3) $\displaystyle=$
$\displaystyle\bar{z}_{n-1}\left({1-\frac{2}{N^{n}_{in}}}\right).$

The Eqn. (3) can be written as a function of the initial average coordination
number $\bar{z}_{0}$ and the initial number of spheres inside the system
$N^{0}_{in}=N_{0}$:

$\displaystyle\bar{z}_{n}$ $\displaystyle=$
$\displaystyle\left\\{{\begin{tabular}[]{l l}$\bar{z}_{0}$&, $n=0$\\\
$\bar{z}_{0}\prod_{i=1}^{n}\left({1-\frac{2}{N_{0}-i}}\right)$&,
$n=1,2,3,...$,\end{tabular}}\right.$ (6)

where the factor $\frac{2}{N_{0}-i}$ represents the decreasing fraction of the
average coordination number when the $i$-th particle is extracted from the
system.

We included the results of this model as dashed lines in Fig. 4a. Strikingly,
the statistical rule is satisfied most of the time, implying that the general
hexagonal configuration is mostly preserved, except in the vacancies
themselves.

A similar model for the average elastic energy can be proposed if we consider
that the elastic energy is proportional to the coordination number. This means
that when we extract a particle, a certain amount of the potential elastic
energy $\Delta E_{e}$ is cut down proportionally to the contacts deleted,
then:

$\displaystyle\bar{E}_{n}$ $\displaystyle=$
$\displaystyle\bar{E}_{0}\frac{\bar{z}_{n}}{\bar{z}_{0}}.$ (7)

However, it is reasonable to consider other dissipation sources during the
particle’s extraction process due to non-conservative forces—for example,
viscous damping and Coulomb friction [24, 18]. Associating those extra losses
to the quantity $\Delta E_{d}$, we can define the nondimensional factor

$\displaystyle\lambda$ $\displaystyle=$ $\displaystyle\frac{\Delta
E_{e}+\Delta E_{d}}{\Delta E_{e}},$ (8)

and introduce it in the model, with $\lambda\geq 1$. Thereby, the elastic
energy evolution is written as:

$\displaystyle\bar{E}_{n}$ $\displaystyle=$
$\displaystyle\left\\{{\begin{tabular}[]{l l}$\bar{E}_{0}$&, $n=0$\\\
$\bar{E}_{0}\prod_{i=1}^{n}\left({1-\frac{2\lambda}{N_{0}-i}}\right)$&,
$n=1,2,3,...$\end{tabular}}\right.$ (11)

The evolution model for the average elastic energy of the system is plotted in
Fig. 4b. Two $\lambda$ values were used: 1) $\lambda=1$, which represents that
the only mechanism for losing elastic energy was the loss of contacts when we
remove a particle; 2) $\lambda=1.5$, obtained from averaging the fitting
values from all the energy curves presented in Fig. 4b. Note in Fig. 4 that
the elastic energy before the avalanche presents a clear tendency to increase
linearly for larger systems, while the coordination number before collapse is
noisy. This observation impulsed us to model the elastic energy (instead of
just the coordination number) to predict the critical amount of particles to
be extracted before avalanches.

Figure 5: (a) Critical average elastic energy as a function initial elastic
energy. The segmented line shows a linear fit with an energy ratio of
$\bar{E}_{c}/\bar{E}_{0}\approx 0.77$. (b) Critical extraction number as a
function the initial number of spheres in the system. The segmented line shows
the model prediction according Eqn. (12) using $\lambda=1.5$. Angles (18, 20,
22) are considered, which correspond to range where the model here presented
is valid.

###  III.4  Prediction for critical extraction number and packing fraction

The model given by Eqn. (11) implies a continuous weakening process, where the
extraction of the $n$-th particle produces a drop in the average elastic
energy with respect to the step before of $2\lambda/(N_{0}-n)\bar{E}_{n-1}$.
This process ends at the extraction number $n_{c}$, when the elastic energy
reaches a critical value $\bar{E}_{c}$, for which the avalanche is triggered.
From Eqn. (11), considering $i\ll N_{0}$, we obtain:

$\displaystyle n_{c}$ $\displaystyle=$
$\displaystyle\frac{\ln\left({\frac{\bar{E}_{c}}{\bar{E}_{0}}}\right)}{\ln\left({1-\frac{2\lambda}{N_{0}}}\right)}.$
(12)

In order to predict the critical extraction number $n_{c}$, Eqn. (12) needs a
relationship between the initial and final value of the average elastic
energy, respectively $\bar{E}_{c}$ and $\bar{E}_{0}$. After initialization at
$n=0$, $\bar{E}_{0}$ is collected. On the other end, $\bar{E}_{c}$ is picked
up just before the avalanche, where the mean elastic energy becomes ill-
defined (see vertical lines in Fig. 4b). Following this procedure and
averaging the results of three simulation at each point, we found a linear
relation between $\bar{E}_{c}$ and $\bar{E}_{0}$ for a range
$\theta\in[$$,$$]$, as shown in Fig. 5a. Note that these points where
$\bar{E}_{c}$ and $\bar{E}_{0}$ are linearly related are shadowed in the phase
diagram (see Fig. 2b). In Fig. 5a the dashed line represents the fitting
function:

$\bar{E}_{c}=\bar{E}_{0y}+0.77\bar{E}_{0},$ (13)

with $\bar{E}_{0y}=4.5\cdot 10^{-8}\text{ J}$, giving approximately
$\bar{E}_{c}/\bar{E}_{0}~{}\approx~{}0.77$. This ratio and the value
$\lambda=1.5$ (fitted from simulations in Fig. 4b) are included in Eqn. (12),
allowing to make a prediction for $n_{c}$ as a function of $N_{0}$, shown as a
dashed line in Fig. 5b. The agreement between the prediction of $n_{c}$ and
the numerical experiment is remarkable. Also, the numerical experiments
confirm that, for a restricted range of angle ($\theta\in[$$,$$]$), the
critical extraction number is $\theta$ independent, as predicted by our model.
The independence on $\theta$ is striking yet natural: if the tilting angle
increases, both the initial elastic energy and the critical energy at collapse
increase similarly.

For this $\theta$ range, the linearity of $n_{c}$ implies that avalanches
triggered by vacancies occur for a constant packing fraction $\phi_{c}$:

$\displaystyle\phi_{c}$ $\displaystyle=$
$\displaystyle\phi_{HE}\left({1-\frac{n_{c}}{N_{0}}}\right),$ (14)

where $\phi_{HE}=\pi/(2\sqrt{3})$ represents the packing fraction of the
hexagonal configuration for a 2D monolayer of disks. Now considering the
$n_{c}/N_{0}$ value from the slope of the model relation in Fig. 5b (which
fits to $n_{c}=0.087N_{0}$ for the range of $N_{0}$ considered), we arrive to
$\phi_{c}=0.83$, which in turn corresponds to $\bar{z}_{c}\approx 4.2$ (see
Fig. 4.a), and an accumulated loss of elastic energy around $23\%$ just before
the avalanche. A possible interpretation of $\phi_{c}$ is that subtracting
particles is analogous to dilating the system continuously. In this context,
$\phi_{c}$ corresponds to a threshold for the medium to flow.

##  IV  Conclusions

We conducted discrete element simulations to study avalanches triggered by
vacancies in a granular medium formed by a monolayer of crystallized spheres
on a tilted plane. By varying the inclination angle of the plane and the size
of the system, a phase diagram was built where three zones are distinguished:
no avalanche, avalanches triggered by extractions, and instant avalanches. We
observed that, as particles are extracted, the average coordination number
($\bar{z}$) decreases, weakening the structure of the system, thus causing the
average elastic energy ($\bar{E}$) also to decrease. On the other hand, we
showed there is a region of the phase space ($\theta\in[$$,$$]$) where the
number of extractions before an avalanche can be predicted. We proposed a
model based on the evolution of the average coordination number and the
elastic energy to predict the critical extraction number for this zone. We
found that the model agrees with the simulations for systems of different
sizes. The critical extraction number predicted is associated with a unique
critical packing fraction of $\phi_{c}=0.83$, associated with a loss of
elastic energy of $23\%$. Our numerical simulations and the proposed model
shed some light on a problem that is difficult to tackle experimentally. It
remains to be seen how our results could be extended to other situations, such
as a non-crystallized monolayer or a 3D configuration. As a final thought, we
would like to emphasize the value of considering everyday phenomena to inspire
new pathways to learn about complex problems like avalanches.

###### Acknowledgements.

The authors thank Nader Droguett, Gabriel Maureira, and Thomas Olivares for
performing preliminary experiments that inspired this study. The authors also
thank project ANID ING2030 16ENI2-71940 (E.R.), as well as Fondecyt Grants:
11190900 (V.S.), 11200464 (G.C.) and 11191106 (P.G.) for financially
supporting this research.

## References

  * Daerr and Douady [1999] A. Daerr and S. Douady, Two types of avalanche behaviour in granular media, Nature 399, 241 (1999).
  * Lastakowski  _et al._ [2015] H. Lastakowski, J.-C. Géminard, and V. Vidal, Granular friction: Triggering large events with small vibrations, Scientific Reports 5, 1 (2015).
  * Salinas  _et al._ [2021] V. Salinas, C. Quiñinao, S. González, and G. Castillo, Triggering avalanches by transverse perturbations in a rotating drum, Scientific Reports 11, 1 (2021).
  * Daniels and Hayman [2008] K. E. Daniels and N. W. Hayman, Force chains in seismogenic faults visualized with photoelastic granular shear experiments, Journal of Geophysical Research 113, B11411 (2008).
  * Le Bouil  _et al._ [2014] A. Le Bouil, A. Amon, S. McNamara, and J. Crassous, Emergence of Cooperativity in Plasticity of Soft Glassy Materials, Physical Review Letters 112, 246001 (2014).
  * Barés  _et al._ [2017] J. Barés, D. Wang, D. Wang, T. Bertrand, C. S. O’Hern, and R. P. Behringer, Local and global avalanches in a two-dimensional sheared granular medium, Physical Review E 96, 052902 (2017).
  * Lherminier  _et al._ [2019] S. Lherminier, R. Planet, V. Levy dit Vehel, G. Simon, L. Vanel, K. J. Måløy, and O. Ramos, Continuously Sheared Granular Matter Reproduces in Detail Seismicity Laws, Physical Review Letters 122, 218501 (2019).
  * hermanoportes [2016] hermanoportes, Market, juice, fruits - free photo on pixabay. pixabay license (2016), [Accessed June 2, 2022].
  * Dorbolo [2005] S. Dorbolo, Stability limit of a granular monolayer, The European Physical Journal E 17, 77 (2005).
  * Liu  _et al._ [1995] C. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, S. Majumdar, O. Nayaran, and T. A. Witten, Force fluctuations in bead packs, Science 269, 513 (1995).
  * Jaeger  _et al._ [1996] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Granular solids, liquids, and gases, Reviews of Modern Physics 68, 1259 (1996).
  * Majmudar and Behringer [2005] T. S. Majmudar and R. P. Behringer, Contact force measurements and stress-induced anisotropy in granular materials, Nature 435, 1079 (2005).
  * Wensrich  _et al._ [2021] C. M. Wensrich, E. H. Kisi, V. Luzin, A. Rawson, and O. Kirstein, Evolution of a contact force network in a 2D granular assembly: an examination using neutron diffraction, Granular Matter 23, 70 (2021).
  * Wyart [2005] M. Wyart, On the rigidity of amorphous solids, Annales de Physique 30, 1 (2005).
  * Lumay and Vandewalle [2005] G. Lumay and N. Vandewalle, Experimental Study of Granular Compaction Dynamics at Different Scales: Grain Mobility, Hexagonal Domains, and Packing Fraction, Physical Review Letters 95, 028002 (2005).
  * Scheller  _et al._ [2006] T. Scheller, C. Huss, G. Lumay, N. Vandewalle, and S. Dorbolo, Precursors to avalanches in a granular monolayer, Physical Review E 74, 031311 (2006).
  * Dorbolo  _et al._ [2011] S. Dorbolo, T. Scheller, F. Ludewig, G. Lumay, and N. Vandewalle, Influence of a reduced gravity on the volume fraction of a monolayer of spherical grains, Physical Review E 84, 041305 (2011).
  * Andreotti  _et al._ [2013] B. Andreotti, Y. Forterre, and O. Pouliquen,  _Granular media: between fluid and solid_ (Cambridge University Press, 2013).
  * Levy dit Vehel  _et al._ [2021] V. Levy dit Vehel, T. Hatano, L. Vanel, K. J. Måløy, and O. Ramos, Dilation as a precursor in a continuous granular fault, EPJ Web of Conferences 249, 15006 (2021).
  * Cundall and Strack [1979] P. A. Cundall and O. D. Strack, A discrete numerical model for granular assemblies, geotechnique 29, 47 (1979).
  * Note [1] For each extraction process, all the spheres in the system are labeled, and one of them is chosen according to a uniform random distribution.
  * Šmilauer et al. [2021] V. Šmilauer et al.,  _Yade Documentation_ , The Yade Project, 3rd ed. (2021), doi:10.5281/zenodo.5705394 (http://yade-dem.org/doc/).
  * Rojas  _et al._ [2019] E. Rojas, V. Vergara, and R. Soto, Case study: Discrete element modeling of wear in mining hoppers, Wear 430-431, 120 (2019).
  * Luding [2008] S. Luding, Introduction to discrete element methods: basic of contact force models and how to perform the micro-macro transition to continuum theory, European journal of environmental and civil engineering 12, 785 (2008).