# Understanding Magnetism in Double Double Perovskites: A Complex Multiple
Magnetic Sublattice System

Anita Halder1,2 Shreya Das1 Prabuddha Sanyal3 Tanusri Saha-Dasgupta1
<EMAIL_ADDRESS>1Department of Condensed Matter Physics and Material
Sciences, S.N. Bose National Centre for Basic Sciences, JD Block, Sector-III,
Salt Lake City, Kolkata 700 106, India  2 School of Physics, Trinity College
Dublin, Dublin, Ireland.  3 Maulana Abul Kalam Azad University of Technology,
Kolkata, India.

###### pacs:

75.50.-y,71.20.-b,75.10.Dg

Understanding magnetism in multiple magnetic sublattice system, driven by the
interplay of varied nature of magnetic exchanges, is on one hand challenging
and on other hand intriguing. Motivated by the recent synthesis of
AA${}^{{}^{\prime}}$BB${}^{{}^{\prime}}$O6 double double perovskites with
multiple magnetic ions both at A- and B-sites, we investigate the mechanism of
magnetic behavior in these interesting class of compounds. We find that the
magnetism in such multiple sublattice compounds is governed by the interplay
and delicate balance between two distinct mechanisms, a) kinetic energy-driven
multiple sublattice double exchange mechanism and b) the conventional super-
exchange mechanism. The derived spin Hamiltonian based on first-principles
calculations is solved by classical Monte Carlo technique which reproduces the
observed magnetic properties. Finally, the influence of off-stoichiometry, as
in experimental samples, is discussed. Some of these double double perovskite
compounds are found to possess large total magnetic moment and also are found
to be half-metallic, which raises the hope of future applications of these
large magnetic moment half-metallic oxides in spintronics and memory devices.

Perovskite structured ABO3 transition metal oxides remained the holy grail of
condensed matter physics due to wide range of fascinating properties exhibited
by them, which includes properties like high temperature superconductivity,
colossal magneto-resistance, half-metallicity etc.C. N. R Rao(1989) ; AS
Bhalla(2000)  With an aim to tailor properties further, one of the common
route is cation substitution. Substitution and 1:1 ordering of cations in B
sublattices in rock-salt arrangement give rise to A2BB${}^{{}^{\prime}}$O6
double perovskites.King (2010) ; Tanusri(2013) ; Sami Vasala(2015) ;
Tanusri(2020)  The topic of magnetism in transition metal oxides with two
magnetic ions as in double perovskite structure has received significant
attention.D.D. Sarma(2000) ; K.-I. Kobayashi(1998) ; Hena Das(2008) ; DP1 ;
DP2 ; DP3 ; DP4 ; Prabuddha Sanyal(2009) ; kato(2007) ; Krockenberger(2007) ;
Nyrissa S. Rogado(2005) ; Hena Das(2009) ; Prabuddha Sanyal(2017) ; Hena
Das(2011) ; Anita(2019) ; kartik(2015)

In this backdrop, it is an interesting issue to ask, what happens if the
compounds involve even larger number of magnetic ions, i.e. more than two
magnetic ions, as in double double perovskites. Due to the structural and
compositional flexibility of the structure, perovskites can accommodate almost
all of the elements of the periodic table, and also can support various
possible coordination. In recent time double double perovskites of general
formula
AA${}^{{}^{\prime}}_{0.5}$A${}^{{}^{\prime\prime}}_{0.5}$BB${}^{{}^{\prime}}$O6
have been synthesized using high pressure and temperature,E. Solana-
Madruga(2016) ; McNally (2017) ; E. Solana-Madruga(2018)  combining columnar
ordering at A sublattice and rock-salt ordering at B sublattice with five
independent cation sites, A, A${}^{{}^{\prime}}$, A${}^{{}^{\prime\prime}}$, B
and B${}^{{}^{\prime}}$, hosting rare-earth or alkaline-earth ion at A site,
3d transition metals at A${}^{{}^{\prime}}$, A`` and B sites, and 5d
transition metal at B${}^{{}^{\prime}}$ site. Use of high pressure is able to
stabilize small magnetic transition-metal ions such as Mn2+ at the A-sites of
perovskites in place of large, nonmagnetic cations like Ca2+ and Sr2+, with
reduced coordination of tetrahedral (4) and square planar (4) instead of usual
dodecahedral (12) coordination of A site.A. J. Dos santos-Garca(2015)  This
introduces a source of magnetism at A-site, which in turn drives the interplay
of magnetism between multiple sublattices, resulting in highly enriched
magnetic properties.

One would naively expect presence of multiple magnetic ions with multiple
magnetic exchange would lead to a frustrating situation and spin glass like
ground state. Contrary to this expectation, recently synthesized double double
perovskites CaMnMReO6 (M=Ni,Co) are found to be magnetically ordered, Elena
Solana-Madruga(2019)  showing ferromagnetic ordering in CaMnNiReO6 with
parallel alignment of spins, which is found to change to ferrimagnetic
ordering when Ni is replaced by Co, the neighboring element in periodic table.
We note the net moment of such multi-component ferromagnetic system is very
high, paving the way to design large moment magnetic oxides. This makes the
situation rather curious in the sense, what makes the three or more magnetic
sublattice system CaMnNiReO6 ferromagnetic, and why replacement of Ni by Co,
the neighboring element in periodic table, makes it ferrimagnetic. What is the
driving mechanism of magnetism in such multi sublattice magnetic system ?
Understanding of such a complex, multiple magnetic sublattice system is
expected to bring out rich physics, which would help future designing of such
oxides.

Motivated by these developments, we present here a first-principles density
functional theory (DFT) based study of these compounds which takes into
account the structural and chemical details in an accurate way, followed by
construction of DFT-derived spin Hamiltonian, which is solved with Monte Carlo
(MC) simulation. Our study uncovers a novel exchange mechanism to be operative
in these compounds, which turn out to be a combination of multi-sublattice
hybridization or kinetic energy-driven double-exchange mechanism, and the more
conventional super-exchange mechanism, the nature of ground state magnetic
order being decided by the competition of these two. While for CaMnNiReO6, the
multi-sublattice hybridization-driven double-exchange mechanism wins over the
super-exchange mechanism stabilizing the long range ferromagnetic state, the
replacement of Ni by Co, increases the core spin value at B sublattice from
S=1 to S=3/2, thus toggling the balance between two exchange mechanisms,
favoring the long range ferrimagnetic behavior. The spin-Hamiltonian with
parameters derived in a first-principles manner provide good description of
measured magnetic properties.Elena Solana-Madruga(2019)  Introduction of off-
stoichiometry is found to maintain the magnetic ground state, encouraging
exploration of many more candidates in such multiple magnetic sublattice
double perovskites. Interestingly, the large magnetic moment compounds,
arising due to long range ferromagnetic ordering between multiple magnetic
sublattices, also turned out to be half-metallic having important implication
for spintronic applications.

##  I Results

###  I.1 Crystal Structure

Fig. 1 shows the four formula unit tetragonal, P42/n crystal structure of
stoichiometric CaMnNiReO6 (CMNRO).Elena Solana-Madruga(2019)  CaMnCoReO6
(CMCRO) compound is isostructural to CaMnNiReO6. The structure consists of
four magnetic sublattices, tetrahedrally coordinated 3d transition metal (TM)
Mn1 at A${}^{{}^{\prime}}$ site, square planar coordinated 3d transition metal
Mn2 at A`` site, octahedrally coordinated 3d transition metal Ni/Co at B site
and octahedrally coordinated 5d transition metal Re at B${}^{{}^{\prime}}$
site, making it a 3d-5d TM magnetic system. Mn1 is connected to two nearest
neighbour (NN) Mn2 sites through Mn1-O-O-Mn2 superexchange paths, while it is
connected to 4 NN Ni(Co)/Re through Mn1-O-Ni(Co)/Re super-exchange paths.
Ni(Co) and Re are connected to each other through corner shared Ni(Co)-O-Re,
of bond angles 141-152o.

###  I.2 Electronic Structure

We analyze the electronic structure of the studied compounds in terms of spin-
polarized density of states, and its projection to orbital characters which
provide us the information on charge and spin states of the transition metal
ions. The GGA+$U$ density of states (DOS), with choice of $U$ = 5 (2) eV and
$J_{H}$ = 0.9 (0.4) eV at 3d TM (Re) sites, for CMNRO and CMCRO are shown in
top and bottom panels of Fig. 2, respectively. In agreement with experimental
findings, the ground state of CMNRO is found to be ferromagnetic, with moments
at three 3d TM sublattices Mn1, Mn2 and Ni sites aligned in parallel
direction, while the moment at Re site is found to be aligned opposite to the
moments at Mn1, Mn2 and Ni sites. The calculated magnetic moments at Mn (Mn1
and Mn2), Ni and Re are found to be 4.5 $\mu_{B}$, 1.6 $\mu_{B}$ and 0.5
$\mu_{B}$ respectively, with a large total moment of 24 $\mu_{B}$ in the unit
cell. On the contrary the ground state of CMCRO is found to be ferrimagnetic,
as observed experimentally, with moments of Mn1, and Mn2 aligned in
antiparallel direction, and Co moment pointing in the direction of Mn1. The Re
moment is found to be antiparallel to Mn1 and Co, with calculated moments
values of 4.5 $\mu_{B}$ (Mn1 and Mn2), 2.6 $\mu_{B}$ (Co) and 0.5 $\mu_{B}$
(Re) and total moment of 8 $\mu_{B}$ in the unit cell. The calculated moments
are in conformity with nominal 2+ valence of Mn1 and Mn2 with high spin d5
occupancy, 2+ valence of Ni/Co with high spin (HS) d8/d7 occupancy and 6+
valence of Re with d1 occupancy. Following this, in DOS of CMNRO we find Mn1
and Mn2 states are filled in the majority spin channel and empty in the
minority channel. Ni eg DOS in CMNRO get filled in majority spin channel and
empty in minority, while Ni t2g states are filled in both spin channels. The
partially filled Re t2g states in CMNRO, with one electron in the minority
spin channel and strongly hybridized with Mn1/Mn2 d and Ni eg states, crosses
the Fermi level making the solution metallic in minority spin channel and
gaped in the majority spin channel. This half metallic solution persists in
CMCRO, though Mn1 and Mn2 d states now become filled and empty, respectively
in two opposite spin channels and Co t2g becomes partly empty. This points to
possibility of achieving spin-dependent nature of the carrier scattering in
these compounds, with a large spin value, which would allow for the resistance
of these large moment compounds to be strongly influenced by the low magnetic
field.

Effect of spin-orbit coupling (SOC) was checked, which is expected to be
appreciable for 5d TM element, Re. The qualitative results are found to remain
unchanged upon inclusion of SOC, apart from an appreciable orbital moment of
$\sim$ 0.15 $\mu_{B}$ that develops at Re site, antiparallel to its spin
moment.

###  I.3 Mechanism of Magnetism

In order to shed light on the mechanism of magnetism in this interesting class
of compounds, we derive the low energy spin Hamiltonian out of DFT inputs. For
this purpose, we perform muffin tin orbital based downfolding calculationsO.
K. Andersen(2000)  that integrate out degrees of freedom which are not of
interest in an energy selective manner. Wannier representation of the
downfolded Hamiltonian provides the estimates of the onsite energies and the
hopping interactions between the orbitals retained in the basis during the
process of downfolding. In the first step of downfolding calculations, we
retain the Mn1, Mn2 d states, Ni eg/ Co d states and Re t2g in the basis and
integrate out the rest. In second step, Mn1, Mn2 and Ni/Co degrees of freedom
are downfolded retaining only the Re t2g degrees of freedom in the basis. The
latter massive downfolding provides the estimates of the Re t2g onsite
energies renormalized by the hybridization from Mn1, Mn2 and Ni/Co states.
Thus the onsite matrix elements of the real space Hamiltonian defined in the
first and second step of downfolding calculations, give the energy level
positions before and after switching on the hybridization between
Mn1/Mn2/Ni(Co) and Re states, respectively. Results of two step downfolding
calculations for CMNRO and CMCRO are presented in top and bottom panels of
Fig. 3, respectively. Mn1 -d, Mn2-d, Ni eg (Co -d ) and Re t2g states are both
crystal field split and exchange split. In distorted tetrahedral coordination,
Mn1 d states are split into 1-1-1-2 fold degeneracies, while Mn2 d states in
square planar coordination are split into 2-2-1 fold degeneracies. The
trigonal distortion in ReO6 octahedra splits its t2g states in 1-2 fold
degeneracies.

Examination of Fig. 3 reveals several interesting aspects, key to construct
the low energy spin Hamiltonian. First of all, in Mn1-Mn2-Ni(Co)-Re basis, Re
t2g states are essentially non-magnetic with negligible exchange splitting.
Secondly, the Re t2g states lie within the exchange split states of Mn1 d, Mn2
d and Ni eg/Co d. Third and most importantly, upon switching on the
hybridization between Mn1 d/Mn2 d/Ni eg(Co d) and Re t2g, captured through
massive downfolding procedure, an exchange splitting of 0.6-0.8 eV is induced
among Re t2g states, with the direction of spin splitting opposite to that of
Mn1 or Mn2 or Ni/Co. This essentially establishes the hybridization-driven
multi sublattice double-exchange process to be operative, in which a negative
spin splitting in the essentially non-magnetic site is induced through
hybridization between the localized spin and itinerant electrons.D.D.
Sarma(2000) ; K.-I. Kobayashi(1998)  In the present context, this can be
captured in terms of a 3+1 sublattice Kondo Lattice model, consisting of a) a
large core spin at the Mn1, Mn2 and Ni(Co) sites, b) strong coupling on the
Mn1/Mn2/ Ni (Co) site between the core spin and the itinerant electron,
strongly preferring one spin polarization of the itinerant electron, and c)
delocalization of the itinerant electron on the Mn1-Mn2-Ni(Co)-Re network, in
the similar spirit as in Prabuddha Sanyal(2009) ,

$\displaystyle H_{DE}$ $\displaystyle=$
$\displaystyle\epsilon_{B}\sum_{i\sigma}b_{i\sigma}^{\dagger}b_{i\sigma}+\epsilon_{Mn1}\sum_{i}m^{1\dagger}_{i\sigma}m^{1}_{i\sigma}$
$\displaystyle+$
$\displaystyle\epsilon_{Mn2}\sum_{i}m^{2\dagger}_{i\sigma}m^{2}_{i\sigma}+\epsilon_{Re}\sum_{i}r^{\dagger}_{i\sigma}r_{i\sigma}$
$\displaystyle+$ $\displaystyle
t_{B-Re}\sum_{<ij>}(b_{i\sigma}^{\dagger}r_{j\sigma}+h.c.)$ $\displaystyle+$
$\displaystyle t_{Mn1-Re}\sum_{<ij>}(m_{i\sigma}^{1\dagger}r_{j\sigma}+h.c.)$
$\displaystyle+$ $\displaystyle
t_{Mn2-Re}\sum_{<ij>}(m_{i\sigma}^{2\dagger}r_{j\sigma}+h.c.)$
$\displaystyle+$ $\displaystyle J_{B}\sum_{i\in B}\vec{S}_{i}^{B}\cdot
b_{i\alpha}^{\dagger}\vec{\sigma}_{\alpha\beta}b_{i\beta}$ $\displaystyle+$
$\displaystyle J_{Mn1}\sum_{i\in A^{\prime}}\vec{S}_{i}^{Mn1}\cdot
m_{i\alpha}^{1\dagger}\vec{\sigma}_{\alpha\beta}m_{i\beta}^{1}$
$\displaystyle+$ $\displaystyle J_{Mn2}\sum_{i\in
A^{\prime\prime}}\vec{S}_{i}^{Mn2}\cdot
m_{i\alpha}^{2\dagger}\vec{\sigma}_{\alpha\beta}m_{i\beta}^{2}$

The $m$’s refer to the Mn sites and the $b$’s to the B (Ni/Co) sites.
$t_{B-Re}$, $t_{Mn1-Re}$, $t_{Mn2-Re}$ represent the nearest neighbor B-Re,
Mn1-Re, Mn2-Re hoppings respectively, with onsite elements $\epsilon_{B}$,
$\epsilon_{Mn1}$, $\epsilon_{Mn2}$ and $\epsilon_{Re}$. The ${\bf S}_{i}$ are
‘classical’ (large $S$) core spins at the Mn1/Mn2/B sites, coupled to the
itinerant Re electrons through a coupling $J$, when the Re electron hops onto
the respective sublattice.

It is to be noted that in $H_{DE}$, the Kondo coupling parameter $J$ is
present only at the magnetic Mn1, Mn2 and B sites, which possess a large core
spin (S=5/2 at Mn1 and Mn2, and S=1 for Ni and S=3/2 for Co), with which the
itinerant Re electron interacts when it hops onto these magnetic sites. The
$J/W$ ratio between the Kondo exchange coupling, $J$ and the bandwidth, $W$ is
thus relevant only on the magnetic Mn1, Mn2 and B (Ni/Co) sites. Following the
DFT inputs, the calculated $J/W$ ratios for the Mn1, Mn2 and Ni sites in CMNRO
are found to be 3.77, 2.06 and 2.7 respectively, while the ratios for the Mn1,
Mn2 and Co sites in CMCRO are given by 3.529, 1.87 and 3.5 respectively. Thus
the exchange coupling $J$ is appreciably larger than the bandwidth $W$ for all
the magnetic sites. Moreover, since the bandwidth, $W\approx z\times t$ where
$t$’s are the hopping parameters appearing in the Hamiltonian, and $z$ is the
number of neighbors, the $J/t$ ratios are even larger, of the order of 8 or
10, justifying use of $J\rightarrow$ $\infty$ model. The $J\rightarrow$
$\infty$ limit of double exchange models was used by Anderson and
HasegawaAnderson-hasegawa  and later studied by P. De GennesGennes  in the
context of perovskites. This limit was studied in the context of double
perovskites in Ref.DP3, . $J\rightarrow$ $\infty$ approximation has been used
for other magnetic perovskite and double perovskite compounds in the
literature for similar $J/t$ ratios.manganites ; millis ; Prabuddha
Sanyal(2009) ; DP3

Invoking the $J\rightarrow$ $\infty$ approximation, one can derive the
effective spin Hamiltonian for CMNRO in terms of the core spins at Mn1
(S=5/2), Mn2 (S=5/2) and Ni (S=1) site as given in the following. The details
of the derivation can be found in the supplementary information (SI).

$\displaystyle H^{{}^{\prime}}_{DE}$ $\displaystyle=$ $\displaystyle
4D_{Mn1-Mn2}\sum_{<ij>i\in A^{\prime},j\in
A^{\prime\prime}}\sqrt{\frac{1+\mathbf{S}^{Mn1}_{i}\cdot\mathbf{S}^{Mn2}_{j}}{2}}$
(2) $\displaystyle+$ $\displaystyle 8D_{Mn1-Ni}\sum_{<ij>i\in A^{\prime},j\in
B}\sqrt{\frac{1+\mathbf{S}^{Mn1}_{i}\cdot\mathbf{S}^{Ni}_{j}}{2}}$
$\displaystyle+$ $\displaystyle 8D_{Mn2-Ni}\sum_{<ij>i\in
A^{\prime\prime},j\in
B}\sqrt{\frac{1+\mathbf{S}^{Mn2}_{i}\cdot\mathbf{S}^{Ni}_{j}}{2}}$

A similar Hamiltonian can be written for the Co compound, with
$\mathbf{S}^{Ni}$ (S=1) replaced by $\mathbf{S}^{Co}$ (S=3/2), where the
coupling constants are $D_{Mn1-Mn2}$, $D_{Mn1-Co}$ and $D_{Mn2-Co}$.

The above described multi sublattice double-exchange Hamiltonian although is
capable of describing the ferromagnetic state of CMNRO, does not account for
the fact that replacement of Ni by Co in CMCRO changes ferromagnetic state to
ferrimagnetic state. This suggests that together with HDE another source of
magnetism needs to be considered. Indeed, there exists another source of
magnetism, namely the super-exchange between the half-filled Mn1-d, Mn2-d, Ni
eg, and high spin Co (t2g+eg) states. The Goodenough-Kanamori ruleGoodenough
states that superexchange interactions are antiferromagnetic where the virtual
electron transfer is between overlapping orbitals that are each half-filled,
but they are ferromagnetic where the virtual electron transfer is from a half-
filled to an empty orbital or from a filled to a half-filled orbital.
Following this, the super-exchange contributions are all antiferromagnetic in
nature (cf Fig. 4) with its strength defined by hopping integrals ($t$) and
onsite energy differences ($\Delta$), $J$ $\propto$
$\sum_{m,m^{{}^{\prime}}}t_{m,m^{{}^{\prime}}}^{2}/$ (U+
$\Delta_{m,m^{{}^{\prime}}}$), where $m$ and $m^{{}^{\prime}}$ are the
orbitals at site $i$ (Mn1/Mn2/Ni (Co)) and $j$ (Mn1/Mn2/Ni (Co)).

The net Hamiltonian for CMNRO adding the contribution of double exchange and
super-exchange is thus given by,

$\displaystyle H$ $\displaystyle=$ $\displaystyle 4D_{Mn1-Mn2}\sum_{<ij>i\in
A^{\prime},j\in
A^{\prime\prime}}\sqrt{\frac{1+\mathbf{S}^{Mn1}_{i}\cdot\mathbf{S}^{Mn2}_{j}}{2}}$
(3) $\displaystyle+$ $\displaystyle 8D_{Mn1-Ni}\sum_{<ij>i\in A^{\prime},j\in
B}\sqrt{\frac{1+\mathbf{S}^{Mn1}_{i}\cdot\mathbf{S}^{Ni}_{j}}{2}}$
$\displaystyle+$ $\displaystyle 8D_{Mn2-Ni}\sum_{<ij>i\in
A^{\prime\prime},j\in
B}\sqrt{\frac{1+\mathbf{S}^{Mn2}_{i}\cdot\mathbf{S}^{Ni}_{j}}{2}}$
$\displaystyle+$ $\displaystyle 4J_{Mn1-Mn2}\sum_{<ij>i\in A^{\prime},j\in
A^{\prime\prime}}\mathbf{S}^{Mn1}_{i}\cdot\mathbf{S}^{Mn2}_{j}$
$\displaystyle+$ $\displaystyle 8J_{Mn1-Ni}\sum_{<ij>i\in A^{\prime},j\in
B}\mathbf{S}^{Mn1}_{i}\cdot\mathbf{S}^{Ni}_{j}$ $\displaystyle+$
$\displaystyle 8J_{Mn2-Ni}\sum_{<ij>i\in A^{\prime\prime},j\in
B}\mathbf{S}^{Mn2}_{i}\cdot\mathbf{S}^{Ni}_{j}$

and a similar one for CMCRO.

In order to estimate the various coupling constants, $D_{Mn1-Mn2}$,
$D_{Mn1-Ni/Co}$, $D_{Mn2-Ni/Co}$ and $J_{Mn1-Mn2}$, $J_{Mn1-Ni/Co}$ and
$J_{Mn2-Ni/Co}$, we apply a two step process. In the first step, we apply
downfolding procedure O. K. Andersen(2000)  to construct a spin unpolarized
Mn1-Mn2-Ni(Co) Hamiltonian defined in effective Mn1-d, Mn2-d, Ni eg (Co d)
basis. The real space representative of this Hamiltonian provides the estimate
of onsite matrix elements of Mn1-d, Mn2-d, Ni eg (Co d) and the hopping
interactions (assumed to be nearest neighbour) between them. Following this,
$J_{Mn1-Mn2}$, $J_{Mn1-Ni/Co}$ and $J_{Mn2-Ni/Co}$ were estimated using the
super-exchange formula, $J$ = $\varSigma_{m,m^{{}^{\prime}}}$ 2 $\times$
$t_{m,m^{{}^{\prime}}}^{2}/$ (U+ $\Delta_{m,m^{{}^{\prime}}}$). In the second
step, the total energies for different possible spin configurations at Mn1,
Mn2 and Ni(Co) sites are calculated, and mapped on to the spin Hamiltonian
given in Eqn. 2. Putting the values of $J_{Mn1-Mn2}$, $J_{Mn1-Ni/Co}$ and
$J_{Mn2-Ni/Co}$ obtained from super-exchange formula, the estimates of
$D_{Mn1-Mn2}$, $D_{Mn1-Ni/Co}$ and $D_{Mn2-Ni/Co}$ are obtained. The estimated
values of $D$’s, and $J$’s, for the two compounds are given in Table I.

For CMNRO, we find that effective Mn1-Mn2, Mn1-Ni and Mn2-Ni interactions are
all negatively signed, i.e ferromagnetic, in conformity with the ferromagnetic
ground state found in the experiment, as well as in DFT total energy
calculations. Similarly for CMCRO, we find Mn1-Mn2 and Mn2-Co effective
interactions are positively signed i.e antiferromagnetic, while Mn1-Co
interaction is ferromagnetic in conformity with its ferrimagnetic ground
state.

Inspecting Table I, we further find while the strength of Mn1-Mn2 super-
exchange ($J_{Mn1-Mn2}$) remains similar between the two compounds, the
strength of Mn1/Mn2 - B super-exchange is greatly enhanced in CMCRO compared
to CMNRO, $J_{Mn1-Ni/Co}$ being enhanced by a factor of 1.6 and
$J_{Mn2-Ni/Co}$ being enhanced by a factor of 3.2. This is expected due to the
fact that while for Ni two unpaired eg electrons participate in the super-
exchange process, for Co, three unpaired electrons, two belonging to eg
manifold and one belonging to t2g contribute. At the same time, we find a
significant weakening of the hybridization-driven exchange between Mn1-Mn2,
reduced by two orders of magnitude compared to Ni compound. These important
changes turn the net interaction to be antiferromagnetic between Mn1 and Mn2,
and that between Mn2 and Co, compared to all ferro interaction for Ni
compound.

###  I.4 Monte Carlo Study of the Spin Hamiltonian

In order to evaluate the finite temperature properties of the defined spin
Hamiltonian, described by Eqn. 3, we perform Monte Carlo simulations. The
total energy of a particular spin configuration can be obtained from the spin
Hamiltonian by plugging in input parameters $D$’s and $J$’s, as listed in
Table I. The spin configurations at Mn1, Mn2 and Ni/Co sites are generated
through Metropolis algorithm in a 3$\times$3$\times$3 unit cell simulation box
with periodic boundary condition. Starting from an initial temperature of 400
K (1000 K) for CMNRO (CMCRO) the simulation temperature is stepped down to T =
1 K with an interval of 2 K. Hundred thousand Monte Carlo steps are employed
to ensure a large sample space, while the physical quantity like magnetization
is calculated by averaging over last 10,000 Monte Carlo steps. The
magnetizations plotted as a function of temperature for CMNRO and CMCRO are
shown in Fig. 5.

For CMNRO compound, our Monte Carlo simulation correctly reproduces the ground
state of this compound where the spins of Mn1. Mn2 and Ni are all found to be
aligned in parallel to each other (cf top, left panel, Fig. 5). We note that
the total moment at low temperature is found out to be 28$\mu_{B}$/unit cell,
corresponding to the sum of the nominal moment 5$\mu_{B}$ for 2 Mn1 and 2 Mn2
with the nominal moment 2$\mu_{B}$ for 4 Ni sites. Transition temperature (Tc)
may be obtained from the inflection point of the derivative of magnetization
versus temperature curve, as shown in the top panel, Fig 5. The Tc is found to
be 142 K which is close to experimentally reported value of 158 K.Elena
Solana-Madruga(2019)  For CMCRO compound, the ferrimagnetic ground state is
also correctly captured. In this case, Mn2 is found to be antiparallel to Mn1
and Co moment, giving rise to the total moment of 12 $\mu_{B}$/unit cell,
arising from magnetic moment of 3$\mu_{B}$ in 4 Co sites and cancellation of
moments at Mn1 and Mn2 sites. The transition in the case of CMNRO is
noticeably sharper compared to CMCRO, as reflected in the narrower width of
the inverse peak in dM/dT curve. Additionally in CMCRO, a shoulder is observed
in the left of the inverse peak which is completely absent in CMNRO.

By repeating the calculation with larger simulation cell size of
4$\times$4$\times$4 (see SI), we find the peak and shoulder structure in CMCRO
is robust, which arises due to the competition between effective FM Mn1-Co
interaction and the two effective AFM Mn1-Mn2 and Mn1-Co interactions. To
demonstrate the effect of competing nature of magnetic interactions in dM/dT
curve of Co compound, we further present the dM/dT curve for varying
$D_{Mn1-Co}$ value in the inset of bottom right panel of Fig. 5. As found,
upon reducing $D_{Mn1-Co}$ from DFT estimated value of -143.9 meV to -135.9
meV, i.e. weakening ferro interaction, the high temperature peak is shifted to
lower temperature, along with redistribution of weight between the peak and
the shoulder, converting the shoulder to a peak. Thus the shoulder feature is
reminiscent of second peak which is not resolved when the strength of ferro
and antiferro interactions are comparable. This implies the high temperature
peak feature arises from the ferro interaction, with the lower temperature
feature arising due to antiferro interaction. The experimental studyElena
Solana-Madruga(2019)  on CaMnReCoO6 reports only magnetic susceptibility, and
do not report dM/dT. However reported experimental dM/dT data for multiple
magnetic ion containing Nd2NiMnO6 double perovskite does exhibit such two
feature structure.das  Such two feature dM/dT curve is also seen for
ferrimagnetic compound NiCr2O4.nicro  We notice with choice of DFT estimated
value of $D_{Mn1-Co}$, the second feature appears around 200 K, very close to
experimentally reported TC of 188 K.Elena Solana-Madruga(2019)

###  I.5 Effect of Off-stoichiometry

The discussion above involves stoichiometric compounds, while the experimental
samples of CMNRO and CMCRO are reported to be off-stoichiometric.Elena Solana-
Madruga(2019)  The experimental samples show high degree of B-site cation
ordering with nominal antisite disorder of 3.4 and 2.5$\%$ for Co and Ni
compounds respectively. This has been attributed to high degree of charge
contrast between (Mn/Co/Ni)2+ and Re6+. However, for CMCRO, while there is
96$\%$ Co at the octahedral B site, 30-40$\%$ of Co was reported to substitute
Mn at the A-sites, leading to an overall Co-rich composition of
CaMn0.7Co1.3ReO6 as opposed to the stoichiometric formula of CaMnCoReO6.
Similarly, for CMNRO, there is an overall Ni-poor composition of
CaMn1.2Ni0.8ReO6 in the experimental sample with some of the Mn atoms
occupying Ni sites in B sublattice. We therefore, need to check whether the
above theoretical understanding also holds good for the off-stoichiometric
compounds.

In order to mimic these off-stoichiometric situations, we replace one out of
the four Ni atoms in the unit cell by Mn, giving rise to Ni poor composition
CaMn1.25Ni0.75ReO6, close to experimental composition. Since all four Ni sites
are equivalent in the unit cell, any one chosen out of four possible sites,
give rise to same results. Similarly, for CMCRO, an extra Co atom replacing
one of the four Mn atoms in the unit cell is introduced, giving rise to
composition CaMn0.75Co1.25ReO6. Total energy calculations show that Co prefers
to occupy the square planar Mn site (Mn2) over Mn1. Interestingly it is found
that for CMNRO even in presence of off-stoichiometry the ground state remains
ferromagnetic with Mn1, Mn2, Mn@Ni and Ni spins aligned in parallel. This
highlights the dominant role of hybridization-driven magnetism, as opposed to
super-exchange driven mechanism, which depends primarily on the positioning of
energy levels, and not on the exchange pathways. Similarly, for CMCRO, even in
presence of off-stoichiometry, Mn1 and Mn2 spins continue to remain
antiparallel, while the Co spins (both at A`` site and B site) are found to be
oppositely aligned to Mn2. The magnetic ground states are thus found to be
robust, and remain unaltered even in presence of off-stoichiometry, as also
found experimentally.

The computed Mn1-Mn2, Mn1-Ni(Co), Mn2-Ni(Co) exchanges for CaMn1.25Ni0.75ReO6
and CaMn0.75Co1.25ReO6 are found not to change significantly compared to their
stoichiometric counterparts ($\sim$ 3-10$\%$) (see Table I), although
introduction of off-stoichiometry introduces few additional interactions like
in CMNRO, Mn@Ni-Mn1, Mn@Ni-Mn2 replacing some of Ni-Mn1, Ni-Mn2 interactions,
respectively, and in CMCRO, Co@Mn2-Mn1, Co@Mn2-Co, replacing some of Mn2-Mn1,
Mn2-Co interaction, respectively. Computation of these additional interactions
show the signs of effective interactions corresponding to these additional
interactions are the same as those of replacing interactions, with values
within 5-7 $\%$.

This suggests the magnetic transition temperature to be not altered
drastically by the off-stoichiometry effect. To check this explicitly, we
carried out Monte Carlo study of Ni-poor and Co-rich compounds as well. Within
the 3$\times$3$\times$3 unit cell simulation box size, for the off-
stoichiometric compounds, it is possible to have different possible atomic
configurations of Mn@Ni and Co@Mn2. Total energy calculations show that extra
Mn (Co) atoms at Ni (Mn2) sites prefer to be uniformly distributed rather than
being clustered. Considering uniform distribution of extra Mn (Co) atoms in
3$\times$3$\times$3 unit simulation cell, this leads to 152 different
configurations. To take this into account, the MC results are averaged over
atomic configurations.

In conformity with DFT results the ground state is found to be ferromagnetic
for CMNRO and ferrimagnetic for CMCRO. The variation of moment with
temperature results in two cases is shown Fig. 6. In presence of off-
stoichiometry, the saturation moment for CMNRO and CMCRO becomes 31
$\mu_{B}$/unit cell and 20 $\mu_{B}$/unit cell, respectively. The dM/dT curves
presented in insets, show similarity with that found for stoichiometric
compounds, confirming the transition temperature not be significantly effected
by off stoichiometry.

##  II Summary and Outlook

In this communication, we study the magnetism in systems containing multiple
magnetic sublattices. The study has been motivated by synthesis of double
double perovskite compounds of general formula,
AA${}^{{}^{\prime}}_{0.5}$A${}^{{}^{\prime\prime}}_{0.5}$BB${}^{{}^{\prime}}$O6,
having transition metal magnetic ions in both A and B sites. The key findings
of our study are summarized in the following,

$\bullet$ Our theoretical analysis combining first-principles and model
Hamiltonian approaches, uncovers the microscopic origin of counter-intuitive
long range ordered magnetism in double double perovskite compounds containing
3d magnetic ions at A and B sites, and 5d magnetic ions at B${}^{{}^{\prime}}$
sites, which turn out to be an interplay of hybridization-driven multi-
sublattice double exchange and super-exchange mechanism of magnetism.

$\bullet$ This interplay relies on the positioning of the $d$ energy levels as
well as the filling. This in turn, triggers ferromagnetic long range order in
CMNRO compound containing two different Mn ions at A sites, and Ni and Re ions
at B and B${}^{{}^{\prime}}$ sites, while replacement of Ni by Co at B site,
decreasing the filling by one in CMCRO stabilizes ferrimagnetism, accounting
for the experimental observations.Elena Solana-Madruga(2019)

$\bullet$ The spin Hamiltonian, capturing the interplay of hybridization-
driven multi-sublattice double exchange and super-exchange mechanism of
magnetism is parameterized in terms of three exchange constants $D_{Mn1-Mn2}$,
$D_{Mn1-Ni/Co}$, $D_{Mn2-Ni/Co}$ for the hybridization-driven multi-sublattice
double exchange and another three exchange constants $J_{Mn1-Mn2}$,
$J_{Mn1-Ni/Co}$, $J_{Mn2-Ni/Co}$ for the super-exchange, with the parameters
derived from first-principles estimated hopping interactions, onsite energies
and total energies of different spin configurations.

$\bullet$ The computed temperature dependent magnetization by Monte Carlo
reproduces the measured magnetic transition temperature of CMNRO with
reasonable accuracy. For CMCRO, it is found that the competition between ferro
and antiferro nature of effective interactions, manifests as two hump
structure of dM/dT, which should be probed further experimentally.

$\bullet$ The calculations are further extended to off stoichiometric
composition of Ni-poor and Co-rich compounds, in order to mimic the
experimental situation. The magnetic properties are found to be retained even
in presence of off-stoichiometry, since the hybridization-driven multi-
sublattice double exchange, a dominant contributor in exchange mechanism of
CMNRO and CMCRO, relies on energy level positioning, rather than on exchange
pathways.

The proposed theory of magnetism being general in nature, should be applicable
to multi sublattice mixed 3d-4d/5d transition metal systems, where one of the
transition metal element is a large band width 4d or 5d element with exchange
splitting significantly smaller than the band width. With appropriate choice
of 3d and 4d/5d elements, this opens up the possibility of stabilization of a
large moment ferromagnetic state. Our first-principles calculation shows this
ferromagnetic state with high value of magnetization is furthermore half-
metallic which should be an attractive possibility for spintronics
applications.

##  III Methods

The first-principles DFT calculations are carried out using the plane-wave
pseudo-potential method implemented within the Vienna Ab-initio Simulation
Package (VASP).G. Kresse(1996)  The exchange-correlation functional is
considered within the generalized gradient approximation (GGA).PBE  The
projector-augmented wave (PAW) potentialsP. E. Blochl(1994)  are used and the
wave functions are expanded in the plane-wave basis with a kinetic-energy cut-
off of 600 eV. Reciprocal-space integration is carried out with a $k$-space
mesh of 6 $\times$ 6 $\times$ 6\. The exchange-correlation beyond GGA is
treated within GGA+$U$ approach with local Coulomb interaction parameterized
in terms of Hubbard $U$ and Hund’s coupling $J_{H}$ within the multi-band
formulation.S. L. Dudarev(1998)  For double-counting correction, fully
localized limit (FLL) of double-counting is considered since the around mean
field (AMF) double-counting is found to give magnetic states a significantly
larger energy penalty than that by the FLL counterpart.doublecounting  The
parameters of GGA+$U$ calculations are chosen as $U$ = 5 eV and $J_{H}$ = 0.9
eV for Ni/Co as appropriate for 3d TM atoms, and U = 2 eV and $J_{H}$ = 0.4 eV
for Re as appropriate for 5d TM atoms.U-value  The U values are varied over
1-2 eV and the qualitative results are found to remain unchanged.

In order to extract a few-band tight-binding (TB) Hamiltonian out of the full
DFT calculation, $N$-th order muffin tin orbital $N$MTO) calculations are
carried out.O. K. Andersen(2000)  A prominent feature of this method is the
downfolding scheme. Starting from a full DFT calculation, it defines a few-
orbital Hamiltonian in an energy-selected, effective Wannier function basis,
by integrating out the degrees of freedom that are not of interest. The
$N$-MTO technique relies on the self-consistent potential parameters obtained
out of linear muffin-tin orbital (LMTO)lmto  calculations.

The magnetization data, obtained from the Monte Carlo simulation of the spin
Hamiltonian, are calculated on a N$\times$N$\times$N unit cell of Mn and Ni/Co
atoms. Finite size effect has been checked. The results presented in the
manuscript are obtained from 3$\times$3$\times$3 lattice simulations. The
periodic boundary conditions are applied during the simulation. The magnetic
transition temperatures are estimated from these calculations.

##  IV Acknowledgments

The authors acknowledge the support of DST Nano-mission for the computational
facility used in this study. T.S-D acknowledges J.C.Followship (grant no.
JCB/2020/000004) for research support.

##  V Author contributions statement

A.H., S.D., P.S. did the theoretical calculations. A.H. and S.D. have equal
contributions. The figures were made by A.H. and S.D. The results were
analyzed by T.S.-D, A.H., S.D., and P.S. T.S.-D wrote the manuscript. A.H.,
S.D., P. S. and T.S.-D. finalized the manuscript.

##  VI Competing financial interests

The authors declare no competing financial interests.

## References

  * (1) C. N. R. Rao Annu. Rev. Phys. Chem., 40, 291 (1998).
  * (2) A. S. Bhall, R. Guo, R. Roy,Mater Res Innov. 4, 3 (2000).
  * (3) G. King, P. M. Woodward, J. Mater. Chem., 20, 5785 (2010).
  * (4) T. Saha-Dasgupta, J. Superconductivity and Novel Mag., 26, 1991 (2013).
  * (5) S. Vasala, M. Karppinen, Progress in Solid State Chemistry, 43, 1 (2015).
  * (6) T. Saha-Dasgupta, Materials Research Express, 7, 1 (2020).
  * (7) D. D. Sarma, P. Mahadevan, T. Saha Dasgupta, S. Ray, A. Kumar, Phys. Rev. Lett., 85, 2549 (2000).
  * (8) K.-I. Kobayashi, T. Kimura, H. Sawada, K. Terakura, Y. Tokura, Nature (London), 395, 677 (1998).
  * (9) H. Das, U. V. Waghmare, T. Saha-Dasgupta, D. D. Sarma, Phys. Rev. Lett., 100, 186402 (2008).
  * (10) A. Chattopadhyay and A. J. Millis, Phys. Rev. B 64, 024424 (2001).
  * (11) L. Brey, M. J. Calderon, S. Das Sarma and F. Guinea, Phys. Rev. B, 74, 094429 (2006).
  * (12) J.L.Alonso,L.A. Fernandez, F. Guinea, F. Lesmes, and V Martin-Mayor, Phys. Rev. B, 67, 214423 (2003).
  * (13) J.B. Phillip, P. Majewski, L. Alff, A. Erb, R. Gross, T. Graf, M.S.Brandt, J. Simon, T. Walther, W. Mader, D. Topwal, and D.D. Sarma, Phys. Rev.B, 68, 144431 (2003).
  * (14) P. Sanyal, P. Majumdar, Phys. Rev. B, 80, 054411 (2009).
  * (15) H. Kato et. al. App. Phys. Lett, 81, 328 (2002).
  * (16) Y. Krockenberger et. al., Phys. Rev. B, 75, 020404 (2007).
  * (17) J. Li, A. W. Sleight, M. A. Subramanian, Adv. Mater., 17, 2225 (2005).
  * (18) H. Das, U. V. Waghmare, T. Saha-Dasgupta, D. D. Sarma, Phys. Rev. B, 79, 144403 (2009).
  * (19) P. Sanyal, Phys. Rev. B, 96, 214407 (2017).
  * (20) H. Das, P. Sanyal, T. Saha-Dasgupta, D. D. Sarma, Phys. Rev. B, 83, 104418 (2011).
  * (21) A. Halder, P. Sanyal, T. Saha-Dasgupta, Phys. Rev. B, 99, 020402 (R) (2019).
  * (22) K Samanta, P Sanyal, T Saha-Dasgupta Sci. Rep. 5, 15010 (2015).
  * (23) E. Solana-Madruga, A. M. Arévalo-Lóṕez, A. J. Dos santos-Garcıá, E. Urones-Garrote, D. Avila-Brande, R. Sáéz-Puche, J. P. Attfield, Angew. Chem., Int. Ed., 55, 9340 (2016).
  * (24) G. M. McNally, A. M. Arévalo-López, P. Kearins, F. Orlandi, P. Manuel, J. P. Attfield, Chem. Mater., 29, 8870 (2017).
  * (25) E. Solana-Madruga, A. M. Arévalo-López, A. J. Dos santos-Garcıá, C. Ritter, C. Cascales, R. Sáez-Puche, J. P. Attfield, Phys. Rev. B, 97, 134408 (2018).
  * (26) E. Solana-Madruga, A. J. Dos santos-Garca, A. M. Arvalo Lýpez, D. Avila-Brande, C. Ritter, J. P. Attfield, J. Saez-Puche, R. Dalton Trans, 44, 20441 (2015).
  * (27) E. Solana-Madruga, Y. Sun, A. M. Arévalo-López, J. P. Attfield, Chem. Commun., 55, 2605 (2019).
  * (28) O. K. Andersen, T. Saha-Dasgupta, Phys. Rev. B, 62, R16219 (2000).
  * (29) A. J. Millis, B. I. Shraiman and R. Mueller, Phys. Rev. Lett, 77, 175 (1996); T. V. Ramakrishnan, H. R. Krishnamurthy, S. R. Hassan, and G. Venketeswara Pai, Phys. Rev. Lett. 92, 157203 (2004).
  * (30) A. Chattopadhyay and A. J. Millis, Phys. Rev. B, 64, 024424, (2001).
  * (31) P.W. Anderson and H. Hasegawa, Phys. Rev. 100, 675 (1955).
  * (32) P. -G de Gennes, Phys. Rev. 118, 141 (1960).
  * (33) P. W. Anderson, Phys. Rev., 79 350 (1950); J. B. Goodenough, Phys. Rev. 100 564 (1955); J. Kanamori, J. Phys. Chem. Solids. 10 87 (1959).
  * (34) R. Das, P. Yanda, A.Sundaresan and D.D. Sarma, Mater. Res. Express6 116122 (2019).
  * (35) A. Ali, G. Sharma, Y. Singh, arXiv:1811.07836.
  * (36) G. Kresse, J. Furthmller, Computational Materials Science, 6(1), 15 (1996).
  * (37) J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett., 77, 3865 (1996); ibid, 78, 1396(E) (1991).
  * (38) P. E. Blöchl, Phys. Rev. B, 50, 17953 (1994).
  * (39) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, A. P. Sutton, Phys. Rev. B, 57(3), 1505 (1998).
  * (40) E. R. Ylvisaker, W. E. Pickett, and K. Koepernik, Phys. Rev. B 79, 035103 (2009).
  * (41) I. V. Solovyev, P. H. Dederichs, V. I. Anisimov, Phys. Rev. B, 50, 16861 (1994).
  * (42) O. K. Andersen, O. Jepsen, Phys. Rev. Lett., 53, 2571 (1994).

Table 1: Estimates of the coupling constants connecting the core spins for the double-exchange and super-exchange mechanisms operative in CMNRO and CMCRO, as estimated employing the super-exchange formula and total energy calculations of different spin configurations. Negative (positive) signs indicate ferro-(antiferro-) magnetic interactions. CMNRO |  $D$ (meV) |  $JS^{2}$ (meV) | Effective
---|---|---|---
Mn1-Mn2 | -91.7 (-94.7) | 47.5 (48.8) | -44.2 (-45.9)
Mn1-Ni | -123.5 (-117.3) | 54.8 (56.1) | -68.7 (-61.2)
Mn2-Ni | -28.6 (-29.6) | 10.8 (9.7) | -17.8 (-19.9)
Mn@Ni-Mn1 | -100.9 | 37.9 | -63.0
Mn@Ni-Mn2 | -50.3 | 30.1 | -20.2
CMCRO |  $D$ (meV) |  $JS^{2}$ (meV) | Effective
Mn1-Mn2 | 1.9 (2.1) | 48.5 (44.6) | 50.4 ( 46.7)
Mn1-Co | -143.9 (-147.4) | 87.9 (86.1) | -56.0 (-61.3)
Mn2-Co | -18.7 (-19.8) | 41.0 (43.1) | 22.3 (23.3)
Co@Mn2-Mn1 | 2.4 | 47.7 | 50.1
Co@Mn2-Co | -150.3 | 88.1 | -62.2
Figure 1: (Color online) Crystal structure of the stoichiometric CMNRO
compound. CMCRO compound is isostructural to CMNRO. The left panel shows the
three dimensional network of four magnetic ions in the structure with Mn at
tetrahedral site (Mn1), Mn at square planar site (Mn2), Ni and Re atoms marked
as red, blue, yellow and green coloured balls respectively. The right panel
shows the oxygen coordination of the four magnetic ions, tetrahedral for Mn1,
square planar for Mn2, and octahedral for Ni and Re. Figure 2: (Color online)
The GGA+$U$ density of states of CMNRO (top) and CMCRO (bottom) projected to
Mn1 d (black, solid), Mn2 d (black, dashed), Ni eg/ Co d (red, solid) and Re
t2g (green, solid) states. The zero of the energy is fixed at the Fermi
energy. Figure 3: (Color online) The energy level diagram for CMNRO (top) and
CMCRO (bottom) considering Mn1-d/Mn2-d/Ni eg (Co d)/Re t2g in basis
(Hybridization-off) and in the massively downfolded Re t2g only basis
(Hybridization-on). See text for details. Figure 4: Super-exchange
interactions in CMNRO and CMCRO between half-filled $d$ states of Mn1, Mn2,
half-filled $e_{g}$ states of Ni and half-filled $e_{g}$ and one of the
$t_{2g}$ states of Co. The fully filled levels not contributing in super-
exchange are not shown. Figure 5: (Color online) Magnetic properties of CMNRO
(left) and CMCRO (right) obtained from Monte-Carlo simulation. Top, left
panels show the ground state magnetic structures, while the top, right panels
show the derivative of magnetization as a function of temperature, the minimum
corresponding to the transition temperature of the corresponding compound.
Mn1, Mn2, Ni(Co) atoms are represented by red, blue and light pink (yellow)
balls respectively. The lower panels shows the magnetization plotted as a
function of temperature. The inset in lower, right panel shows the shift of
transition temperature (the minima of the curves) for monotonic decrease of
$D_{Mn1-Co}$. For details see text. Figure 6: (Color online) Magnetic
properties of CaMn1.25Ni0.75ReO6 (left) and CaMn0.75Co1.25ReO6 (right)
obtained from Monte-Carlo simulation. Insets show the derivative of
magnetization as a function of temperature.


# From granular collapses to shallow water waves:
A predictive model for tsunami generation

Wladimir Sarlin<EMAIL_ADDRESS>Université Paris-
Saclay, CNRS, Laboratoire FAST, F-91405 Orsay, France  Cyprien Morize
<EMAIL_ADDRESS>Université Paris-Saclay, CNRS,
Laboratoire FAST, F-91405 Orsay, France  Alban Sauret  University of
California, Santa Barbara, Department of Mechanical Engineering, CA 93106, USA
Philippe Gondret  Université Paris-Saclay, CNRS, Laboratoire FAST, F-91405
Orsay, France

###### Abstract

In this article, we present a predictive model for the amplitude of impulse
waves generated by the collapse of a granular column into a water layer. The
model, which combines the spreading dynamics of the grains and the wave
hydrodynamics in shallow water, is successfully compared to a large dataset of
laboratory experiments, and captures the influence of the initial parameters
while giving an accurate prediction. Furthermore, the role played on the wave
generation by two key dimensionless numbers, i.e., the global Froude number
and the relative volume of the immersed deposit, is rationalized. These
results provide a simplified, yet comprehensive, physical description of the
generation of tsunami waves engendered by large-scale subaerial landslides,
rockfalls, or cliff collapses in a shallow water.

Introduction. Tsunamis are among the most destructive natural disasters for
human coastal settlements. While events generated by earthquakes have been
extensively studied [1, 2, 3, 4], several past or potential occurrences of
high amplitude waves arising from large-scale landslides have also been
reported in past decades [5, 6, 7, 8, 9], which constitutes a grand challenge
in environmental fluid mechanics [10]. The 1958 Lituya Bay tsunami, featuring
the highest recorded wave runup of 524 m [6], is reminiscent of the importance
of understanding the physics underlying such events, for reliable hazard
assessments.

A relevant approach is to experimentally model landslides and pyroclastic
flows using a granular material [11, 6, 12, 13, 14, 15, 16, 17, 18]. Although
the finding of constitutive laws for granular media remains challenging, and
is still attracting many research activities [19, 20, 21, 22, 23], canonical
experimental configurations have been developed to study geophysical flows,
such as the granular collapse experiment [24, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34, 35, 36]. In this situation, a column of dry grains, suddenly released,
falls vertically while spreading horizontally under the effect of gravity. The
resulting final deposit height $H_{\infty}$ and runout distance $\Delta
L_{\infty}=L_{\infty}-L_{0}$ both depend on the aspect ratio $a$ of the
column, i.e., the ratio of its initial height $H_{0}$ to width $L_{0}$,
through nontrivial scaling laws [26, 27]. The aspect ratio also influences the
characteristic timescale of the horizontal spreading [32, 36]. In the case of
an immersed granular collapse, similar scaling laws are obtained for the final
morphology of the deposit [37, 38, 39]. An interesting feature of the granular
collapse experiment was its ability to compare successfully with large-scale
geophysical events [11].

This experimental configuration has also been used recently to generate
impulse waves of geophysical interest, by releasing the grains into a water
layer of depth $h_{0}$ [40, 41, 42, 43, 44, 45]. Several dimensionless
parameters have been observed to be important for the wave generation. The
most intuitive is the Froude number, that can be defined in two ways: globally
as $\mathrm{Fr}_{0}=\sqrt{H_{0}/h_{0}}$ [41, 40, 45], which compares the
typical vertical free-fall velocity $\sqrt{gH_{0}}$ of the grains to the
velocity $\sqrt{gh_{0}}$ of linear gravity waves in shallow water, or locally
as $\mathrm{Fr}_{f}=v_{m}/\sqrt{gh_{0}}$, with $v_{m}$ the maximum horizontal
velocity of the granular front at the water surface [42, 44]. In particular,
the wave amplitude was related analytically to $\mathrm{Fr}_{f}$ through
weakly nonlinear scalings [44], or with the linear approximation
$A_{m}/h_{0}\simeq 1.2~{}\mathrm{Fr}_{f}$ [46, 42]. In addition, the wave
amplitude $A_{m}$ was found to be linked to the final immersed volume of
grains $\Delta V_{\infty}$, which can be measured in the field, and allows
back-calculations of former events [43]. However, so far a unifying model
relating the wave amplitude to the initial geometry of the column and the
water depth remains elusive. Such an attempt was made recently through the
study of the collapse of a Newtonian viscous fluid into a water layer [47].
Nevertheless, this configuration is closer to a dam-break problem and does not
account for the granular nature of the landslide. As a result, the prediction
of the wave amplitude following a granular collapse currently relies on
empirical fits, without relevant physical model.

In the present study, a predictive model for the wave amplitude in shallow
water conditions is introduced, by coupling the dynamics of the granular
collapse highlighted in Ref. [36] with the wave hydrodynamics elucidated in
Ref. [44]. The influence of the initial parameters is revealed by the model,
which therefore captures the role played by the initial Froude number
$\mathrm{Fr}_{0}$ and also explains how the final volume of immersed grains
$\Delta V_{\infty}$ is related to the wave generation. The model is
successfully compared with an extended dataset of experiments.

Figure 1: (a) 3D view of the experimental setup, showing the initial granular
column of height $H_{g}$ and width $L_{0}$, the initial fluid depth $h_{0}$,
the channel width $W$, and the maximum wave amplitude $A_{m}$. (b) 2D view of
the final granular deposit of height $H_{\infty}$, runout length $L_{\infty}$,
immersed volume $\Delta V_{\infty}$, and front position $\ell_{\infty}$ at
$z=h_{0}$. (c)-(d) Impulse waves generated by the collapse of a column with
$H_{g}=39$ cm and $L_{0}=10$ cm: (c) Bore wave for $h_{0}=3$ cm, and (d)
solitary wave for $h_{0}=10$ cm (movies are available in Supplemental Material
[48]).

Experimental setup. The two-dimensional experimental setup, sketched in figure
1(a), consists in a $2\ \mathrm{m}\times 0.30\ \mathrm{m}\times 0.15\
\mathrm{m}$ reservoir, filled up to a height $h_{0}$ with water [42, 44]. A
rectangular dry granular column of initial height $H_{g}$ and width $L_{0}$
rests on a solid step of height $h_{0}$, and is retained by a sliding gate.
The following ranges of initial parameters were investigated: $9\
\mathrm{cm}\leqslant H_{g}\leqslant 50\ \mathrm{cm}$, $2.5\
\mathrm{cm}\leqslant L_{0}\leqslant 20\ \mathrm{cm}$, and $2\
\mathrm{cm}\leqslant h_{0}\leqslant 25\ \mathrm{cm}$. The total height
$H_{0}=H_{g}+h_{0}$ from the bottom plane and the aspect ratio $a=H_{0}/L_{0}$
of the initial column were thus varied in the ranges $11\ \mathrm{cm}\leqslant
H_{0}\leqslant 59\ \mathrm{cm}$ and $1\lesssim a\lesssim 21$, respectively.
This allowed to explore a larger range of initial parameters, especially for
the aspect ratio, than in previous studies [41, 42, 44]. The grains are
monodisperse glass beads of mean diameter $d=5$ mm, density $\rho=2.5\
\rm{g.cm^{-3}}$ with a dense packing fraction $\phi\simeq 0.64$. Using the
same grains for all experiments is motivated by a previous study [42], which
showed that the maximum amplitude of the generated wave does not significantly
depend on the grains’ size and density, at least for millimeter-scale grains
denser than water. At the beginning of the experiment, the gate is quickly
lifted at $V=1\ \rm{m.s^{-1}}$ to ensure no significant influence of the
release process on the collapse dynamics [36]. The solid step restricts fluid
perturbations induced by the withdrawal of the gate [41]. The grains fall
under the effect of gravity into water, leading to the formation of an impulse
wave, as illustrated in figures 1(c)-(d). The video recordings are then
processed to extract the maximum amplitude $A_{m}$ reached by the wave, and
the maximum horizontal velocity $v_{m}=\mathrm{d}\ell/\mathrm{d}t$ of the
advancing granular front at the water surface $z=h_{0}$ [44]. In addition, the
final morphology of the deposit is characterized by measuring the final runout
distance $L_{\infty}$, height $H_{\infty}$, volume $\Delta V_{\infty}$ of
immersed grains, and position $\ell_{\infty}$ of the granular front at
$z=h_{0}$, as illustrated in figure 1(b).

Figure 2: Evolution of the maximum wave amplitude $A_{m}$ with the fluid depth
$h_{0}$, for three different initial granular columns of same width $L_{0}$
but different height $H_{g}$. The horizontal dashed lines indicate the
saturation of $A_{m}$ at large $h_{0}$. Full symbols refer to shallow water
waves, while empty symbols correspond to deep water waves. The shaded area
corresponds to deep water waves, where $A_{m}\lesssim 0.45~{}h_{0}$, and the
black dashed line ( ) represents the observed limit between bore and solitary
waves ($A_{m}\simeq 0.96~{}h_{0}$).

Figure 3: (a) Relative final height $H_{\infty}/L_{0}$ and (b) runout distance
$\Delta L_{\infty}/L_{0}$ reached by the final deposit as a function of the
aspect ratio of the initial column $a=H_{0}/L_{0}$. ( ) Power laws
$H_{\infty}/L_{0}=\alpha a^{n}$ [with $\alpha=1$ and $n=1$ when $a\lesssim
0.93$, and $\alpha=0.95$ and $n=1/3$ when $a\gtrsim 0.93$], and $\Delta
L_{\infty}/L_{0}=\beta a^{m}$ [with $\beta=1.65$ and $m=1$ when $a\lesssim 3$,
and $\beta=2.38$ and $m=2/3$ when $a\gtrsim 3$]. ($\boldsymbol{\circ}$)
Experiments of bore and solitary waves which are not significantly affected by
the presence of the solid step. ($\boldsymbol{\times}$) Bore and solitary
waves with a noticeable influence of the solid step, and nonlinear transition
waves.

Influence of the initial parameters on the generated wave and final deposit.
The evolution of the maximum wave amplitude $A_{m}$ with the water depth
$h_{0}$ is presented in figure 2 for three initial geometries of the column.
For a given column geometry, $A_{m}$ first increases with $h_{0}$ at the
shallowest depths, then decreases for intermediate depths, before saturating
at a constant asymptotic value. When the initial fluid depth $h_{0}$ is small
compared to the column height $H_{g}$, the collective motion of the spreading
grains is mainly horizontal, and the advancing front is similar to a piston
pushing the water, as illustrated in figures 1(c)-(d) (see also Supplemental
Material [48]). This situation leads to the formation of either transient bore
($A_{m}\gtrsim 0.96~{}h_{0}$) or solitary ($0.45~{}h_{0}\lesssim A_{m}\lesssim
0.96~{}h_{0}$) waves [44]. In that case, the generation mechanism, i.e., the
advancing granular front pushing water like a moving piston, is in agreement
with recent numerical simulations which reproduced finely a specific
geophysical case with a complex topography [9], although the initial
conditions are quite different than from our model laboratory configuration.
On the contrary, when $h_{0}$ is of the order or greater than $H_{g}$, the
grains have essentially a vertical motion, closer to a “granular spillway”
situation, in which a finite volume of grains discharges into deep water (see
Supplemental Material [48]). This regime departs from shallow water conditions
and leads to the formation of the so-called nonlinear transition waves for
which $A_{m}\lesssim 0.45~{}h_{0}$ [44]. Increasing the initial height $H_{g}$
of the column, while keeping the other initial parameters constant,
systematically leads to higher values for the wave amplitude. In the
following, the study focuses on shallow water waves ($A_{m}\gtrsim
0.45~{}h_{0}$), whose understanding is of significant importance for risk
assessments as bore waves are potentially very destructive in enclosed fjords
scenarios [8], while at the same time solitary waves can propagate and
distribute part of the collapse energy far from the source.

In the present configuration, the generated wave is a direct consequence of
the granular collapse. Therefore, it is interesting to focus also on the main
features of the final granular deposit, characterized by its final height
$H_{\infty}$ and runout distance $\Delta L_{\infty}=L_{\infty}-L_{0}$. In the
dry case, these two final parameters, rescaled by the initial width of the
column, $H_{\infty}/L_{0}$ and $\Delta L_{\infty}/L_{0}$, are known to be
governed mainly by the initial aspect ratio $a$ of the column, and slightly
influenced by the material and frictional properties of the grains and the
substrate [26, 27, 29, 31]. Figures 3(a) and 3(b) show the evolution of
$H_{\infty}/L_{0}$ and $\Delta L_{\infty}/L_{0}$, respectively, as a function
of $a$. In the present experiments, the initial granular column stands on a
solid step of height $h_{0}$, and one may expect that the final deposit may be
different from the one that would occur without a solid step. The initial
dataset is thus subdivided into two parts: ($\boldsymbol{\circ}$) the
experiments for which the solid step has no significant influence on the final
deposit, i.e., where $H_{0}/h_{0}\geqslant 4$ and $h_{0}\leqslant L_{0}$
according to [49], and ($\boldsymbol{\times}$) those impacted by the presence
of the step. For the first set ($\boldsymbol{\circ}$), the relative final
height $H_{\infty}/L_{0}$ [figure 3(a)], and runout distance $\Delta
L_{\infty}/L_{0}$ [figure 3(b)] are then fitted by piecewise power laws of the
initial aspect ratio, following Lajeunesse et al. (2005) [26]:

$\displaystyle\frac{H_{\infty}}{L_{0}}=\alpha a^{n}=$
$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle a\\\\[8.0pt]
\displaystyle
0.95~{}a^{1/3}\end{array}\right.\begin{array}[]{ll}\displaystyle\ \mathrm{for\
}a\lesssim 0.93,\\\\[8.0pt] \displaystyle\ \mathrm{for\ }a\gtrsim
0.93,\end{array}$ (1) $\displaystyle\frac{\Delta L_{\infty}}{L_{0}}=\beta
a^{m}=$ $\displaystyle\left\\{\begin{array}[]{ll}\displaystyle
1.65~{}a\\\\[8.0pt] \displaystyle
2.38~{}a^{2/3}\end{array}\right.\begin{array}[]{ll}\displaystyle\ \mathrm{for\
}a\lesssim 3,\\\\[8.0pt] \displaystyle\ \mathrm{for\ }a\gtrsim 3,\end{array}$

where the prefactors $\alpha$ and $\beta$ are slightly dependent on the
material and frictional properties (except the trivial value $\alpha=1$ for
low enough aspect ratios) [28, 31, 35], and would also be impacted by the
presence of cohesion [14, 50]. Therefore, using different grain or substrate
properties would require one to adjust these coefficients accordingly. The
exponents of Eqs. (1) differ from 2D to 3D configurations [27, 26], and are a
signature of the gravity-driven dynamics of the collapse [30]. It should be
mentioned that although no experimental data are available here at low aspect
ratio ($a\lesssim 0.93$), the scaling $H_{\infty}/L_{0}=a$ is represented, as
this case corresponds to the trapezoidal final shape of the deposit already
described by previous work [26, 27, 30, 34, 36], where only a fraction of the
initial granular column collapses, so that $H_{\infty}=H_{0}$. The second set
($\boldsymbol{\times}$) of experiments deviates significantly from the
scalings of Eqs. (1) for the final height as well as for the run-out distance.
These data are either impacted by the presence of the solid step or depart
from shallow water conditions (deep water waves), and are therefore not
considered in the following.

Figure 4: (a) Relative maximum amplitude $A_{m}/h_{0}$ of the wave as a
function of the relative final extension of the granular front
$\Delta\ell_{\mathrm{\infty}}/h_{0}$, with ( )
$A_{m}/h_{0}=0.50\sqrt{\Delta\ell_{\mathrm{\infty}}/h_{0}}$. (b) Comparison
between the measurement of $\Delta\ell_{\mathrm{\infty}}$ and the model
prediction $\Delta\ell_{\mathrm{\infty,th}}$. ( )
$\Delta\ell_{\mathrm{\infty}}=0.82\Delta\ell_{\mathrm{\infty,th}}$. (c)
$A_{m}$ as a function of $\sqrt{\Delta\ell_{\mathrm{\infty,th}}h_{0}}$, with (
) Eq. (4). Symbols designate ($\triangleright$) bore and ($\square$) solitary
waves. The shaded area in (a) correspond to the deep water waves region where
$A_{m}\lesssim 0.45~{}h_{0}$.

A predictive model for the wave amplitude. To understand the bell-shaped
curves for $A_{m}(h_{0})$ in shallow water, as illustrated in figure 2, the
local Froude number $\mathrm{Fr}_{f}=v_{m}/\sqrt{gh_{0}}$ needs to be related
to the initial parameters of the granular column as it is known that, at first
order, $A_{m}/h_{0}\simeq 1.2~{}\mathrm{Fr}_{f}$ [46, 42]. The granular front
velocity $v_{m}$ is expected to scale as $\Delta\ell_{\infty}/\tau_{x}$, where
$\Delta\ell_{\infty}=\ell_{\infty}-L_{0}$ is the front position of the final
deposit and $\tau_{x}$ the typical spreading time of the granular mass, both
taken at $z=h_{0}$. In a recent study, $\tau_{x}$ was shown to be proportional
to $\sqrt{\Delta\ell_{\infty}/g}$ [36]. Using this result, $\mathrm{Fr}_{f}$
may be estimated as

$\mathrm{Fr}_{f}=\frac{v_{m}}{\sqrt{gh_{0}}}\propto\frac{\Delta\ell_{\infty}}{\tau_{x}\sqrt{gh_{0}}}\propto\sqrt{\frac{\Delta\ell_{\infty}}{h_{0}}},$
(2)

and therefore $A_{m}/h_{0}$ is also expected to scale as
$\sqrt{\Delta\ell_{\infty}/h_{0}}$, due to the aforementioned approximate
linear relation between $A_{m}/h_{0}$ and $\mathrm{Fr}_{f}$. Figure 4(a)
presents the rescaled wave amplitude $A_{m}/h_{0}$ as a function of
$\Delta\ell_{\infty}/h_{0}$ for ($\triangleright$) bore and ($\square$)
solitary waves. All experimental data corresponding to shallow water waves
collapse onto a master curve of equation $A_{m}/h_{0}\simeq
0.5\sqrt{\Delta\ell_{\infty}/h_{0}}$. Therefore, in this case the moving
granular front acts like a piston pushing the water, whose velocity would be
coupled to its stroke [as highlighted by Eq. (2)].

An estimate of the typical horizontal extension
$\Delta\ell_{\mathrm{\infty,th}}$ can be obtained by assuming a triangular
shape for the final deposit: $\Delta\ell_{\mathrm{\infty,th}}=\Delta
L_{\infty}-h_{0}L_{\infty}/H_{\infty}$. Using Eqs. (1) the final runout
distance at $z=h_{0}$ is thus expected to be given by the relation

$\Delta\ell_{\mathrm{\infty,th}}=\beta a^{m-1}H_{0}-\frac{1+\beta
a^{m}}{\alpha a^{n}}h_{0}.$ (3)

In figure 4(b), the measured front position $\Delta\ell_{\infty}$ of the final
deposit is compared with the prediction of Eq. (3). The relation
$\Delta\ell_{\infty}=0.82~{}\Delta\ell_{\mathrm{\infty,th}}$ captures all the
experiments, with a prefactor revealing a systematic overestimate of
$\Delta\ell_{\infty}$ by the model, as the final deposit is not perfectly
triangular, but presents some curvature instead [26, 27, 30, 36]. From the two
fits of figures 4(a) and 4(b), the wave amplitude should be given by

$A_{m}=0.45\sqrt{\Delta\ell_{\mathrm{\infty,th}}h_{0}}.$ (4)

The prediction given by Eq. (4) fits well the data, as presented in figure
4(c), which reveals that the model captures the physics behind the wave
generation in shallow water. By inserting Eq. (3) in Eq. (4), the wave
amplitude can be expressed as a function of the initial parameters $H_{0}$,
$L_{0}$, and $h_{0}$ to obtain

$A_{m}=0.45\,h_{0}\,\sqrt{\beta\left(\frac{H_{0}}{L_{0}}\right)^{m-1}\frac{H_{0}}{h_{0}}-\frac{1+\beta(H_{0}/L_{0})^{m}}{\alpha(H_{0}/L_{0})^{n}}},$
(5)

where the coefficients $\alpha$, $\beta$, $n$ and $m$ are given in Eqs. (1).

To highlight the influence of the initial parameters on the wave amplitude,
figures 5(a)-(b) report the prediction of $A_{m}$ given by Eq. (5) as a
function of $h_{0}$ for $a=3$ and different values of $H_{0}$ from 20 to 60 cm
[figure 5(a)], and for $H_{0}=40$ cm and different values of $a$ from 1 to 10
[figure 5(b)]. In all cases, we observe that $A_{m}\sim{h_{0}}^{1/2}$ for the
shallowest depths, then reaches a maximum value before decreasing. Note that
this maximum, as well as the critical depth at which it occurs, both increase
linearly with $H_{0}$ when $a$ is constant [see Appendix for more details].
The present model reproduces thus well the bell-shaped part of the curves
obtained in figure 2 which were also observed in previous experimental studies
[16, 42].

Figure 5: (a)-(b) Maximum wave amplitude $A_{m}$, given by Eq. (5), as a
function of the water depth $h_{0}$, for (a) different initial height $H_{0}$
with the same aspect ratio $a=3$, and (b) different $a$ with $H_{0}=40$ cm.
(c)-(d) $A_{m}/h_{0}$ as a function of (c) $\mathrm{Fr}_{0}$ as given by Eq.
(6), and (d) $\Delta V_{\mathrm{\infty,th}}/(W{h_{0}}^{2})$ as given by Eq.
(7), for $H_{0}=40$ cm and different $a$. Shaded areas are the same as in
figure 2.

Figure 6: (a) $A_{m}$ as a function of $h_{0}$, for the experiments of figure
2. (b) $A_{m}/h_{0}$ as a function of the global Froude number
$\mathrm{Fr_{0}}$. (c) Comparison between the experimental and predicted
immersed volumes, with ( ) $\Delta V_{\infty}=0.65\Delta
V_{\mathrm{\infty,th}}$. (d) $A_{m}/h_{0}$ as a function of the relative
volume $\Delta V_{\infty}/(W{h_{0}}^{2})$ of the immersed deposit, with ( )
the empirical law $A_{m}/h_{0}=0.25\left[\Delta
V_{\infty}/(W{h_{0}}^{2})\right]^{0.8}$ from [43]. The colored solid lines in
(a), (b), and (d) correspond to the model predictions for the same initial
granular columns as in the experiments, and the shaded areas are the same as
in figure 2.

Let us now explore the prediction of the model for the wave generation as a
function of two dimensionless numbers, reflecting either the initial state or
the final deposit of the collapse: the global Froude number
$\mathrm{Fr}_{0}=\sqrt{H_{0}/h_{0}}$ [40, 41, 42, 45] and the relative volume
$\Delta V_{\infty}/(W{h_{0}}^{2})$ of the immersed deposit [43]. Using Eq. (5)
straightforwardly leads to

$\frac{A_{m}}{h_{0}}=0.45\sqrt{\beta
a^{m-1}{\mathrm{Fr}_{0}}^{2}-\frac{1+\beta a^{m}}{\alpha a^{n}}}.$ (6)

In addition, within the present model the final immersed deposit would have a
trapezoidal shape, of volume $\Delta V_{\mathrm{\infty,th}}=(\Delta
L_{\infty}+\Delta\ell_{\mathrm{\infty,th}})Wh_{0}/2$. As a result, Eq. (4)
leads to

$\frac{A_{m}}{h_{0}}=0.45\sqrt{2\frac{\Delta
V_{\mathrm{\infty,th}}}{W{h_{0}}^{2}}-\frac{\Delta L_{\infty}}{h_{0}}}.$ (7)

The relative wave amplitude predicted by Eqs. (6) and (7) is presented in
figures 5(c) and 5(d), respectively, for different initial aspect ratios from
1 to 10 and $H_{0}=40$ cm. In both cases, $A_{m}/h_{0}$ increases
monotonically with a sublinear shape. Note that the influence of the aspect
ratio on the relative wave amplitude is almost imperceptible when varying the
rescaled volume of the immersed deposit in figure 5(d), but is non-negligible
when the global Froude number is varied, as in figure 5(c).

Comparison with the experiments. In figure 6(a), the prediction of the wave
amplitude $A_{m}$ using Eq. (5) is compared to the experiments from figure 2.
The numerical prefactors found by the best fit for each initial geometry are
in the range [0.42–0.49], i.e., very close to the value 0.45 of Eq. (5). A
good agreement is observed in shallow water conditions, as the bell-shaped
curves are correctly reproduced. As expected, the present model fails in
predicting the plateau values observed for deep water waves ($A_{m}\lesssim
0.45~{}h_{0}$) since the model is not valid under these conditions. The
relative wave amplitude $A_{m}/h_{0}$ is also reported as a function of
$\mathrm{Fr}_{0}$ in figure 6(b), alongside the model curves obtained for
identical initial granular columns. All experimental data are within the
collection of model curves from Eq. (6), with a sublinear increase of
$A_{m}/h_{0}$ with $\mathrm{Fr}_{0}$ for a given $a$. The influence of the
aspect ratio is also well captured in the investigated range, as shown by the
experiments ($\blacktriangledown$) with the most slender initial columns, that
are significantly below the other data (corresponding to lower values of $a$).
Therefore, at first order, the global Froude number $\mathrm{Fr}_{0}$ governs
the relative wave amplitude, but with a noticeable effect of the initial
aspect ratio. It is also interesting to consider the correlation between the
relative wave amplitude and the relative volume of immersed deposit. To do so,
we first compare the model prediction for the volume of the immersed deposit
$\Delta V_{\mathrm{\infty,th}}$ with the experimental measurements of $\Delta
V_{\infty}$. Figure 6(c) shows that there is a good correlation between them.
A linear trend of slope 0.65 fits well the data, showing that the predicted
volume overestimates the measurements in a systematic way, here again because
the final immersed deposit is not perfectly trapezoidal. Using this result, it
is possible to compare $A_{m}/h_{0}$ to $\Delta V_{\infty}/(W{h_{0}}^{2})$ for
all experiments and to draw the corresponding model curves, as shown in figure
6(d). The agreement between the experiments and the predictions allows one to
explain the strong link existing between $\Delta V_{\infty}/(W{h_{0}}^{2})$
and the generated wave, which was previously reported [43]. This correlation
is a striking result, that may appear surprising, but one needs to keep in
mind that the final deposit is the result of the gravity-driven collapse and
is thus reminiscent of its dynamics. The empirical equation for the wave
amplitude given in [43] is also reported in figure 6(d), to highlight that the
present model better captures the experimental data in the whole range of
shallow water conditions.

Conclusion. The amplitude $A_{m}$ of the impulse wave generated by the
collapse of a granular column in a shallow water of depth $h_{0}$ can be
predicted by combining the spreading dynamics of the grains [36], relating the
initial parameters of the column to the local Froude number $\mathrm{Fr}_{f}$
based on the advancing granular front at the water surface, and the wave
hydrodynamics linking $\mathrm{Fr}_{f}$ to $A_{m}$ [44]. In this situation,
the spreading motion of the grains indeed behaves as a peculiar piston, whose
velocity is coupled with its stroke. The present model explicits the evolution
of $A_{m}$ as a function of the initial parameters $H_{0}$, $L_{0}$, and
$h_{0}$. In addition, it highlights the important role played by the global
Froude number $\mathrm{Fr}_{0}=\sqrt{H_{0}/h_{0}}$ and the relative volume
$\Delta V_{\infty}/(W{h_{0}}^{2})$ of the immersed deposit. It explains why
these different dimensionless numbers have been observed to play a key role in
previous studies [40, 42, 43, 45].

It is worth noting that the numerical prefactors obtained in the model
slightly depend on the considered granular matter. Therefore, one would need
to adjust their values to apply the present model to other materials. Within
these minor adjustments, this predictive model should work in similar
configurations, for instance for the ones explored experimentally [40, 41] or
numerically [45] in previous studies, provided that the shallow water
condition is fulfilled. For deep water conditions, a predictive model remains
to be developed. The present model should be adapted accordingly for the
topographies and bathymetries that are specific in the field, by expressing
the maximum velocity of the granular front in the impact region, to derive the
relevant Froude number in that case. For post-mortem analysis in the field
based on the measurements of the final immersed deposits, the present model
could be used without significant adjustments. Indeed, there is a trace of the
gravity-driven dynamics in the final deposit, as shown by the good correlation
of Eq. (7) with all the data for very different columns with a large range of
investigated aspect ratios. In order to investigate 3D effects that may be
important in real cases, in particular in the collapse of volcanic islands
such as the 2018 Krakatau event [7], the study of waves generated by the
collapse of cylindrical columns instead of rectangular ones would be of great
interest. Indeed, the scaling laws for the final deposit are known to be
different in this case [24, 25]. In addition, a 3D context would make the wave
hydrodynamics more complex.

## Appendix : Critical depth for a maximum wave amplitude

It is interesting here to provide a discussion on the critical water depth at
which the maximum wave amplitude is reached, for the collapse of a granular
column of a given aspect ratio, as observed in the bell-shaped curves of Figs.
5(a)-(b). From Eqs (3)-(4), one can obtain that the critical depth is given by
$h_{0}^{c}=\Delta L_{\infty}H_{\infty}/(2L_{\infty})$, which leads to the
corresponding amplitude $A_{m}^{c}\simeq
0.45\sqrt{\Delta{L_{\infty}}^{2}H_{\infty}/(4L_{\infty})}$. For a column with
a given aspect ratio $a$, leading to a given deposit, we can infer from these
expressions that both $h_{0}^{c}$ and $A_{m}^{c}$ increase linearly with the
initial height $H_{0}$ of the column, as can be observed in Fig. 5(a). To
illustrate this point, fom Eqs. (1) we can write $h_{0}^{c}/H_{0}$ in terms of
the initial aspect ratio $a$ of the granular column in the following manner:

$\displaystyle\frac{h_{0}^{c}}{H_{0}}\simeq\left\\{\begin{array}[]{ll}\displaystyle\frac{0.83a}{1+1.65a}&\mathrm{for\
}a\lesssim 0.93,\\\\[8.0pt]
\displaystyle\frac{0.78a^{1/3}}{1+1.65a}&\mathrm{for\ }0.93\lesssim a\lesssim
3,\\\\[8.0pt] \displaystyle\frac{1.13}{1+2.38a^{2/3}}&\mathrm{for\ }a\gtrsim
3.\end{array}\right.$ (8)

These expressions thus show that $h_{0}^{c}$ is proportional to $H_{0}$ when
the aspect ratio of the column is kept constant.

Figure 7: Critical Froude number $\mathrm{Fr}_{0}^{c}=\sqrt{H_{0}/h_{0}^{c}}$
at which the maximum wave amplitude is reached for the collapse of granular
column, as a function of the initial aspect ratio $a$. The line corresponds to
Eqs. (8) for the three domains separated by the two vertical dotted lines at
$a=0.93$ and $a=3$.

The associated critical Froude number
$\mathrm{Fr}_{0}^{c}=\sqrt{H_{0}/h_{0}^{c}}$ is presented in figure 7 as a
function of the aspect ratio of the column. It shows a nontrivial evolution
with $a$, with a sharp decrease for low aspect ratio ($a\lesssim 0.93$), and
an increasing region for higher values of $a$. It means that, for aspect
ratios below 0.93, the maximum value for the amplitude of the wave is reached
for $h_{0}^{c}$ small when compared to the initial column height $H_{0}$.
Around $a\simeq 0.93$, the curve reaches a global minimum where
$\mathrm{Fr}_{0}^{c}\simeq 1.8$, so that $h_{0}^{c}\simeq H_{0}/3.24$: It is
the highest possible value of $h_{0}^{c}$ when varying the initial aspect
ratio of the column. Then, increasing $a$ beyond this point again implies that
the critical wave amplitude occurs with $h_{0}^{c}$ increasingly small when
compared to $H_{0}$.

###### Acknowledgements.

The authors are grateful to J. Amarni, A. Aubertin, L. Auffray and R. Pidoux
for the elaboration of the experimental setup, and report no conflict of
interest.

## References

  * Tadepalli and Synolakis [1996] S. Tadepalli and C. E. Synolakis, Model for the leading waves of tsunamis, Phys. Rev. Lett. 77, 2141 (1996).
  * Kânoğlu and Synolakis [2006] U. Kânoğlu and C. Synolakis, Initial value problem solution of nonlinear shallow water-wave equations, Phys. Rev. Lett. 97, 148501 (2006).
  * Synolakis and Bernard [2006] C. E. Synolakis and E. N. Bernard, Tsunami science before and beyond Boxing Day 2004, Philos. Trans. Royal Soc. A 364, 2231 (2006).
  * Dutykh  _et al._ [2011] D. Dutykh, C. Labart, and D. Mitsotakis, Long wave run-up on random beaches, Phys. Rev. Lett. 107, 184504 (2011).
  * Ward and Day [2001] S. N. Ward and S. Day, Cumbre Vieja volcano—potential collapse and tsunami at La Palma, Canary Islands, Geophys. Res. Lett. 28, 3397 (2001).
  * Fritz  _et al._ [2009] H. M. Fritz, F. Mohammed, and J. Yoo, Lituya Bay landslide impact generated mega-tsunami 50th anniversary, Pure Appl. Geophys. 166, 153 (2009).
  * Grilli  _et al._ [2019] S. T. Grilli, D. R. Tappin, S. Carey, S. F. L. Watt, S. N. Ward, A. R. Grilli, S. L. Engwell, C. Zhang, J. T. Kirby, L. Schambach, _et al._ , Modelling of the tsunami from the december 22, 2018 lateral collapse of Anak Krakatau volcano in the Sunda Straits, Indonesia, Sci. Rep. 9, 11946 (2019).
  * Waldmann  _et al._ [2021] N. Waldmann, K. Vasskog, G. Simpson, E. Chapron, E. W. N. Støren, L. Hansen, J.-L. Loizeau, A. Nesje, and D. Ariztegui, Anatomy of a catastrophe: Reconstructing the 1936 rock fall and tsunami event in lake Lovatnet, Western Norway, Front. Earth Sci. 9, 671378 (2021).
  * Rauter  _et al._ [2022] M. Rauter, S. Viroulet, S. S. Gylfadóttir, W. Fellin, and F. Løvholt, Granular porous landslide tsunami modelling – the 2014 Lake Askja flank collapse, Nat. Commun. 13, 678 (2022).
  * Dauxois  _et al._ [2021] T. Dauxois, T. Peacock, P. Bauer, C. P. Caulfield, C. Cenedese, C. Gorlé, G. Haller, G. N. Ivey, P. F. Linden, E. Meiburg, N. Pinardi, N. M. Vriend, and A. W. Woods, Confronting grand challenges in environmental fluid mechanics, Phys. Rev. Fluids 6, 020501 (2021).
  * Lajeunesse  _et al._ [2006] E. Lajeunesse, C. Quantin, P. Allemand, and C. Delacourt, New insights on the runout of large landslides in the Valles-Marineris canyons, Mars, Geophys. Res. Lett. 33, L04403 (2006).
  * Roche  _et al._ [2011] O. Roche, M. Attali, A. Mangeney, and A. Lucas, On the run-out distance of geophysical gravitational flows: Insight from fluidized granular collapse experiments, Earth Planet. Sci. Lett. 311, 375 (2011).
  * Viroulet  _et al._ [2014] S. Viroulet, A. Sauret, and O. Kimmoun, Tsunami generated by a granular collapse down a rough inclined plane, EPL 105, 34004 (2014).
  * Langlois  _et al._ [2015] V. J. Langlois, A. Quiquerez, and P. Allemand, Collapse of a two-dimensional brittle granular column: implications for understanding dynamic rock fragmentation in a landslide, J. Geophys. Res. Earth. Surf. 120, 1866 (2015).
  * Lindstrøm [2016] E. K. Lindstrøm, Waves generated by subaerial slides with various porosities, Coast. Eng. 116, 170 (2016).
  * Mulligan and Take [2017] R. P. Mulligan and W. A. Take, On the transfer of momentum from a granular landslide to a water wave, Coast. Eng. 125, 16 (2017).
  * Bullard  _et al._ [2019] G. Bullard, R. Mulligan, and W. Take, An enhanced framework to quantify the shape of impulse waves using asymmetry, J. Geophys. Res. Oceans 124, 652 (2019).
  * Bougouin  _et al._ [2020] A. Bougouin, R. Paris, and O. Roche, Impact of fluidized granular flows into water: Implications for tsunamis generated by pyroclastic flows, J. Geophys. Res. Solid Earth 125, e2019JB018954 (2020).
  * Jop  _et al._ [2006] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, Nature 441, 727 (2006).
  * Bouzid  _et al._ [2013] M. Bouzid, M. Trulsson, P. Claudin, E. Clément, and B. Andreotti, Nonlocal rheology of granular flows across yield conditions, Phys. Rev. Lett. 111, 238301 (2013).
  * Henann and Kamrin [2014] D. L. Henann and K. Kamrin, Continuum modeling of secondary rheology in dense granular materials, Phys. Rev. Lett. 113, 178001 (2014).
  * Kim and Kamrin [2020] S. Kim and K. Kamrin, Power-law scaling in granular rheology across flow geometries, Phys. Rev. Lett. 125, 088002 (2020).
  * Gaume  _et al._ [2020] J. Gaume, G. Chambon, and M. Naaim, Microscopic origin of nonlocal rheology in dense granular materials, Phys. Rev. Lett. 125, 188001 (2020).
  * Lajeunesse  _et al._ [2004] E. Lajeunesse, A. Mangeney-Castelnau, and J. P. Vilotte, Spreading of a granular mass on a horizontal plane, Phys. Fluids 16, 2371 (2004).
  * Lube  _et al._ [2004] G. Lube, H. E. Huppert, R. S. J. Sparks, and M. A. Hallworth, Axisymmetric collapses of granular columns, J. Fluid Mech. 508, 175–199 (2004).
  * Lajeunesse  _et al._ [2005] E. Lajeunesse, J. Monnier, and G. Homsy, Granular slumping on a horizontal surface, Phys. Fluids 17, 103302 (2005).
  * Lube  _et al._ [2005] G. Lube, H. E. Huppert, R. S. J. Sparks, and A. Freundt, Collapses of two-dimensional granular columns, Phys. Rev. E 72, 041301 (2005).
  * Balmforth and Kerswell [2005] N. J. Balmforth and R. R. Kerswell, Granular collapse in two dimensions, J. Fluid Mech. 538, 399–428 (2005).
  * Zenit [2005] R. Zenit, Computer simulations of the collapse of a granular column, Phys. Fluids 17, 031703 (2005).
  * Staron and Hinch [2005] L. Staron and E. Hinch, Study of the collapse of granular columns using two-dimensional discrete-grain simulation, J. Fluid Mech. 545, 1 (2005).
  * Staron and Hinch [2007] L. Staron and E. J. Hinch, The spreading of a granular mass: Role of grain properties and initial conditions, Granul. Matter 9, 205 (2007).
  * Lacaze  _et al._ [2008] L. Lacaze, J. C. Phillips, and R. R. Kerswell, Planar collapse of a granular column: experiments and discrete element simulations, Phys. Fluids 20, 063302 (2008).
  * Lacaze and Kerswell [2009] L. Lacaze and R. R. Kerswell, Axisymmetric granular collapse: A transient 3D flow test of viscoplasticity, Phys. Rev. Lett. 102, 108305 (2009).
  * Lagrée  _et al._ [2011] P.-Y. Lagrée, L. Staron, and S. Popinet, The granular column collapse as a continuum: validity of a two-dimensional Navier–Stokes model with a $\mu$(I)-rheology, J. Fluid Mech. 686, 378 (2011).
  * Man  _et al._ [2021] T. Man, H. E. Huppert, L. Li, and S. A. Galindo-Torres, Deposition morphology of granular column collapses, Granul. Matter 23, 59 (2021).
  * Sarlin  _et al._ [2021a] W. Sarlin, C. Morize, A. Sauret, and P. Gondret, Collapse dynamics of dry granular columns: From free-fall to quasistatic flow, Phys. Rev. E 104, 064904 (2021a).
  * Topin  _et al._ [2012] V. Topin, Y. Monerie, F. Perales, and F. Radjaï, Collapse dynamics and runout of dense granular materials in a fluid, Phys. Rev. Lett. 109, 188001 (2012).
  * Bougouin and Lacaze [2018] A. Bougouin and L. Lacaze, Granular collapse in a fluid: Different flow regimes for an initially dense-packing, Phys. Rev. Fluids 3, 064305 (2018).
  * Jing  _et al._ [2018] L. Jing, G. C. Yang, C. Y. Kwok, and Y. D. Sobral, Dynamics and scaling laws of underwater granular collapse with varying aspect ratios, Phys. Rev. E 98, 042901 (2018).
  * Huang  _et al._ [2020] B. Huang, Q. Zhang, J. Wang, C. Luo, X. Chen, and L. Chen, Experimental study on impulse waves generated by gravitational collapse of rectangular granular piles, Phys. Fluids 32, 033301 (2020).
  * Cabrera  _et al._ [2020] M. A. Cabrera, G. Pinzon, W. A. Take, and R. P. Mulligan, Wave generation across a continuum of landslide conditions from the collapse of partially submerged to fully submerged granular columns, J. Geophys. Res. Oceans 125, e2020JC016465 (2020).
  * Robbe-Saule  _et al._ [2021a] M. Robbe-Saule, C. Morize, R. Henaff, Y. Bertho, A. Sauret, and P. Gondret, Experimental investigation of tsunami waves generated by granular collapse into water, J. Fluid Mech. 907, A11 (2021a).
  * Robbe-Saule  _et al._ [2021b] M. Robbe-Saule, C. Morize, Y. Bertho, A. Sauret, A. Hildenbrand, and P. Gondret, From laboratory experiments to geophysical tsunamis generated by subaerial landslides, Sci. Rep. 11, 18437 (2021b).
  * Sarlin  _et al._ [2021b] W. Sarlin, C. Morize, A. Sauret, and P. Gondret, Nonlinear regimes of tsunami waves generated by a granular collapse, J. Fluid Mech. 919, R6 (2021b).
  * Nguyen [2022] N. H. T. Nguyen, Collapse of partially and fully submerged granular column generating impulse waves: An empirical law of maximum wave amplitude based on coupled multiphase fluid–particle modeling results, Phys. Fluids 34, 013310 (2022).
  * Noda [1970] E. Noda, Water waves generated by landslides, J. Waterw. Harb. Coast. Eng. Div. 96, 835 (1970).
  * Kriaa  _et al._ [2022] Q. Kriaa, S. Viroulet, and L. Lacaze, Modeling of impulse waves generated by a viscous collapse in water, Phys. Rev. Fluids 7, 054801 (2022).
  * [48] See Supplemental Material at [URL will be inserted by publisher] to access movies of the experiments.
  * Robbe-Saule [2019] M. Robbe-Saule,  _Modélisation expérimentale de génération de tsunami par effondrement granulaire_ , Ph.D. thesis, Université Paris Saclay (2019).
  * Li  _et al._ [2021] P. Li, D. Wang, Y. Wu, and Z. Niu, Experimental study on the collapse of wet granular column in the pendular state, Powder Technol. 393, 357 (2021).