*[inlinelist,1]label=(), ETH Zurich, Switzerland and University of Cagliari, Italy<EMAIL_ADDRESS>of Cagliari, Italy and http://blockchain.unica.it bart@unica.ithttps://orcid.org/0000-0003-3796-9774Partially supported by project SERICS (PE00000014) and PRIN 2022 DeLiCE (F53D23009130001) under the MUR National Recovery and Resilience Plan funded by the European Union – NextGenerationEU. Daniele Pusceddu and Massimo Bartoletti [500]Software and its engineering Formal methods [500]Software and its engineering Formal software verification

# Formalizing Automated Market Makers in the Lean 4 Theorem Prover

Daniele Pusceddu  Massimo Bartoletti

###### Abstract

Automated Market Makers (AMMs) are an integral component of the decentralized
finance (DeFi) ecosystem, as they allow users to exchange crypto-assets
without the need for trusted authorities or external price oracles. Although
these protocols are based on relatively simple mechanisms, e.g. to
algorithmically determine the exchange rate between crypto-assets, they give
rise to complex economic behaviours. This complexity is witnessed by the
proliferation of models that study their structural and economic properties.
Currently, most of theoretical results obtained on these models are supported
by pen-and-paper proofs. This work proposes a formalization of constant-
product AMMs in the Lean 4 Theorem Prover. To demonstrate the utility of our
model, we provide mechanized proofs of key economic properties like arbitrage,
that at the best of our knowledge have only been proved by pen-and-paper
before.

###### keywords:

Smart contracts, Ethereum, Verification, Blockchain

###### category:

##  1 Introduction

Automated Market Makers (AMMs) are one of the key applications in the
Decentralized Finance (DeFi) ecosystem, as they allow users to trade crypto-
assets without the need for trusted intermediaries [13]. Unlike traditional
order-book exchanges, where buyers and sellers must find a counterpart, AMMs
enable traders to autonomously swap assets deposited in liquidity pools
contributed by other users, who are incentivized to provide liquidity by a
complex reward mechanism. At the time of writing, there are multiple AMM
protocols controlling several billions of dollars worth of
assets111https://defillama.com/protocols/Dexes. This has made AMMs an
appealing target for attacks, resulting in losses worth billions of dollars
over time222https://chainsec.io/defi-hacks/.

The security of AMMs depends on several factors: besides the absence of
traditional programming bugs, it is crucial that their economic mechanism
gives rise to a rational behaviour of its users that aligns with the AMM ideal
functionality, i.e. providing an algorithmic exchange rate coherent with the
one given by trusted price oracles. Therefore, it is important to obtain
strong guarantees about the economic mechanisms of these protocols. While
formal verification tools for smart contracts based on model-checking are
useful in detecting programming bugs and even in proving some structural
properties of AMMs [8, 10], they are not suitable for verifying, or even
expressing more complex properties regarding the economic mechanism of AMMs.
These economic mechanisms have been studied in several research works, which,
in most cases, provide pen-and-paper proofs of the obtained properties. Given
the complexity of the studied models, it would be desirable to also provide
machine-verified proofs, so that we may rely on the proven properties beyond
any reasonable doubt. To the best of our knowledge, existing mechanized
formalizations [12] focus on verifying relevant structural properties of AMMs
like their state consistency, and not on studying the economic mechanism of
AMMs (see Section 4 for a detailed comparison).

### Contributions

In this paper we formalize Automated Market Makers in the Lean 4 theorem
prover. Our model is based on a slightly simplified version of the Uniswap v2
protocol (one of the leading AMMs), which was studied in [7] with a pen-and-
paper formalization. We provide a Lean specification of blockchain states,
abstracted from any factors that are immaterial to the study of AMMs. Then, we
model the fundamental interactions that users may have with AMMs as well as
the economic notions of price, networth and gain. Finally, we build machine-
checked proofs of economic properties of constant-product AMMs. In particular,
we derive an explicit formula for the economic gain obtained by a user after
an exchange with an AMM. Building upon this formula we prove that, from a
trader’s perspective, aligning a constant-product AMM’s internal exchange rate
with the rate given by the trader’s price oracle implies the optimal gain from
that AMM. This results in the fundamental property of AMMs acting as price
oracles themselves [4]. We then construct the optimal swap transaction that a
rational user can perform to maximize their gain, solving the arbitrage
problem. Our formalization and
proofs333https://github.com/danielepusceddu/lean4-amm are made available in a
public GitHub repository. At the best of our knowledge, this is the first
mechanized formalization of the economic mechanism of AMMs. We finally discuss
some open issues, and alternative design choices for formalizing AMMs.

##  2 Formalization

An Automated Market Maker implements a decentralized exchange between two
different token types. The exchange rate is determined by a smart contract,
which also takes care of performing the exchange itself: namely, the contract
receives from a trader some amount of the input token type, and sends back the
correct amount of the output token type, which is taken from the AMM reserves.
A single smart contract can control many instances of AMMs (also called _AMM
pairs_): we may have a pair for each possible unordered pair of token types.
To create an AMM instance, a user must provide the initial liquidity for the
reserves of that pair of tokens. Liquidity providers are rewarded with a type
of token that specifically represents shares in that AMM’s reserves: we call
these _minted_ token types, while any other token type will be called
_atomic_.

### Blockchain State

We begin by formalizing the blockchain state, abstracting from the details
that are immaterial to the study of AMMs. Then, our model includes the users’
wallets, the AMMs and their reserves (see Section 2). We formalize the
universes of users and atomic token types as the types A and T, respectively,
as structures that encapsulate a natural number. Hereafter, we use
$a,b,\ldots$ to denote users in ${\tt{A}}$, and
$\tau,\tau_{0},\tau_{1},\ldots$ to denote tokens in ${\tt{T}}$. Minted token
types are pairs of T. We represent the funds owned by a user by a _wallet_
that maps token types to non-negative reals. To rule out wallets with infinite
tokens, we use Mathlib’s finitely supported functions444https://leanprover-
community.github.io/mathlib4_docs/Mathlib/Data/Finsupp/Defs.html: in general,
given any type $\alpha$ and any type $M$ with a $0$ element,
$f\in\alpha\to_{0}M$ if $\mathrm{supp}\left(f\right)=\left\\{x\in\alpha\mid
f\left(x\right)\neq 0\right\\}$ is finite.

We define the type ${\tt{W}}_{0}$ of wallets of atomic tokens as a structure
encapsulating ${\tt{T}}\to_{0}\mathbb{R}_{\geq 0}$. This definition induces an
element $0\in{\tt{W}}_{0}$ such that $\mathrm{supp}\left(0\right)=\emptyset$:
this is the _empty wallet_ , which enables us to form the type
${\tt{T}}\to_{0}{\tt{W}}_{0}$. We define the type ${\tt{W}}_{1}$ of wallets of
minted tokens as a structure that encapsulates a function
${\tt{bal}}\in{\tt{T}}\to_{0}{\tt{W_{0}}}$. The intuition is that ${\tt{bal}}\
\tau_{0}\ \tau_{1}$ gives the owned amount of the minted token type created by
the AMM pair with tokens $\tau_{0}$ and $\tau_{1}$. Consistently, the function
bal must satisfy two conditions: ${\tt{bal}}\ \tau_{0}\ \tau_{1}={\tt{bal}}\
\tau_{1}\ \tau_{0}$, meaning that the order of atomic tokens is irrelevant,
and ${\tt{bal}}\ \tau\ \tau=0$, meaning that the two token types in an AMM
must be distinct. Our definition of ${\tt{W}}_{1}$ encapsulates proofs of
these two properties, called unord and distinct, respectively.

We map users to their wallets with the types ${\tt{S}}_{0}$ and
${\tt{S}}_{1}$, which account for the atomic tokens and for the minted tokens,
respectively. Finally, we formalize sets of AMM pairs with the type AMMs. The
definition is strikingly similar to that of ${\tt{W}}_{1}$, but with a changed
constraint. The intuition is that ${\tt{res}}\ \tau_{0}\ \tau_{1}$ gives the
reserves of $\tau_{0}$ in the AMM pair $(\tau_{0}$, $\tau_{1})$, while
${\tt{res}}\ \tau_{1}\ \tau_{0}$ gives the reserves of $\tau_{1}$ in the same
AMM pair. For uninitialized AMM pairs, the obtained reserves must be $0$. The
property posres ensures that either an AMM pair has no reserves of both token
types (i.e., the AMM has not been created yet), or both token types have
strictly positive reserves (i.e., one cannot deplete the reserves of a single
token type in an AMM pair). We combine the previous definitions in the type
$\Gamma$, which represents the state of a blockchain (note that $\Gamma$
abstracts from all the details immaterial for AMMs).

### Token supply

Given a blockchain state $s\in\Gamma$, we define the supply of an atomic token
type $\tau_{0}$ as:

${\tt{atomsupply}}_{s}\left(\tau_{0}\right)\;\;=\;\;\sum_{a\in\mathrm{supp}(s{\tt{.atoms}})}(s{\tt{.atoms}}\;a)\;\tau_{0}\;\;+\;\;\sum_{\tau_{1}\in\mathrm{supp}(s{\tt{.amms}}\,\tau_{0})}(s{\tt{.amms}}\;\tau_{0})\;\tau_{1}$

where the partial application $s{\tt{.amms}}\;\tau_{0}$ gives a map with all
AMM pairs with $\tau_{0}$ as one of their token types. We define the supply of
a minted token type $(\tau_{0},\tau_{1})$ as follows:

${\tt{mintsupply}}_{s}\left(\tau_{0},\tau_{1}\right)\;\;=\;\;\sum_{a\in\mathrm{supp}(s{\tt{.mints}})}(s{\tt{.mints}}\;a)\;\tau_{0}\;\tau_{1}$

The corresponding Lean definitions (in LABEL:list:supplydefs) have been split
in order to facilitate theorem proving and, in particular, the use of Lean’s
simplifier.

Table 1: Fundamental Lean definitions for the state of an AMM system.

##  3 Results

We now present some noteworthy properties of AMMs that we have proven in Lean.

Proposition 3.1 ensures that, in any reachable state, there exists an AMM with
token types $\tau_{0}$ and $\tau_{1}$ if and only if the minted token type
$(\tau_{0},\tau_{1})$ is in circulation, i.e. it has a strictly positive
supply. This result showcases the validity of our model with regards to
reasoning about reachable states. Technically, it also allows us to prove that
${\tt{mintedprice}}_{s}(o,\tau_{0},\tau_{1})$ is strictly positive for any
initialized AMM pair $(\tau_{0},\tau_{1})$ in any reachable state $s$.

######  Proposition 3.1 (Existence of AMMs _vs._ minted token supply).

Let $s^{\prime}\in\Gamma$ be a reachable blockchain state. Then, for any
minted token type $\left(\tau_{0},\tau_{1}\right)$, its supply in $s^{\prime}$
is strictly positive if and only if $s^{\prime}$ has an AMM with token types
$\tau_{0}$ and $\tau_{1}$.

######  Proof 3.2.

By induction on the length of the sequence of transactions leading to
$s^{\prime}$:

  * •

Base case: empty transaction sequence. The proof is trivial for both
directions, since a valid starting state has no initialized AMMs and no minted
tokens in circulation.

  * •

Inductive case: there are several subcases depending on the last transaction
fired in the sequence. Here we consider the creation of the AMM
$(\tau_{0}^{\prime},\tau_{1}^{\prime})$ in reachable state $s^{\prime}$. We
proceed by cases on the truth of the equality
$\left\\{\tau_{0},\tau_{1}\right\\}=\left\\{\tau_{0}^{\prime},\tau_{1}^{\prime}\right\\}$.
If it is a different token pair, then the supply remains unchanged along with
the initialization status, and we can conclude by the induction hypothesis. If
it is the same token pair, then we just incremented its minted token supply
(which is non-negative, so after incrementing it, it must be strictly
positive), and we just initialized the AMM.

  * •

See source code for the other cases. AMMLib/Transaction/Trace.lean:275

Lemma 3.3 allows us to determine the change in the value of a wallet after it
has been updated in some way: the resulting equality is the basis for all the
subsequent proofs. To illustrate it, we briefly introduce two definitions that
have been omitted before: ${\tt{drain_{0}}}(w_{0},\tau_{0})$ is the atomic
token wallet such that ${\tt{drain_{0}}}(w_{0},\tau_{0})(\tau_{0})=0$ and
${\tt{drain_{0}}}(w_{0},\tau_{0})(\tau_{1})=w_{0}(\tau_{1})$ for every other
token $\tau_{1}\neq\tau_{0}$. We define
${\tt{drain_{1}}}(w_{1},(\tau_{0},\tau_{0}))\in{\tt{W_{1}}}$ similarly.

######  Lemma 3.3 (Value expansion).

Let $w_{0}\in{\tt{W}}_{0}$, $o_{0}\in{\tt{T}}\to\mathbb{R}_{>0}$, and
$\tau_{0},\tau_{1}\in{\tt{T}}$. Then, value(w_0, o_0) = w_0(τ_0)⋅o(τ_0) +
value(drain_0(w_0, τ_0), o_0) and, with $w_{1}\in{\tt{W}}_{1}$ and
$o_{1}\in{\tt{T^{2}}}\to\mathbb{R}_{>0}$ such that
$o_{1}(\tau_{0},\tau_{1})=o_{1}(\tau_{1},\tau_{0})$, value(w_1, o_1) =
w_1(τ_0, τ_1)⋅o_1(τ_0, τ_1) + value(drain_1(w_1, τ_0, τ_1), o_1)

######  Proof 3.4.

By definition of value and by properties of the sum over a finite support.
Full Lean proof at AMMLib/State/AtomicWall.lean:116

Lemma 3.5 gives an explicit formula for the gain obtained by a user upon
firing a swap transaction. It is fundamental to all proofs involving gain.

######  Lemma 3.5 (Gain of a swap).

Let $sw\in{\tt{Swap}}\left(sx,s,a,\tau_{0},\tau_{1},x\right)$ and let
$o\in{\tt{T}}\to\mathbb{R}_{>0}$. Then, gain(a, o, s, apply sw ) =
(sw.y⋅o(τ_1) - x ⋅o(τ_0)) ⋅(1 - (s.mints a τ0τ1)mintsupplys(τ0,τ1))

######  Proof 3.6.

By repeated application of Lemma 3.3 in order to isolate the value contributed
by the token types involved in the swap, and by use of Mathlib’s ring_nf
simplifier tactic. AMMLib/Transaction/Swap/Networth.lean:55

Lemma 3.7 establishes a correspondence between the profitability of a swap
transaction (i.e., a positive or negative gain) and the order between the swap
rate and the exchange rate given by the price oracle. In particular, assuming
a trader $a$ who is not a liquidity provider (i.e., $a$ has no minted tokens
for the AMM pair targeted by the swap), Lemma 3.7 states that:

  * •

${\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw\right)<0\;\iff\;sx\left(x,r_{0},r_{1}\right)<\nicefrac{{o\left(\tau_{0}\right)}}{{o\left(\tau_{1}\right)}}$

  * •

${\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw\right)=0\;\iff\;sx\left(x,r_{0},r_{1}\right)=\nicefrac{{o\left(\tau_{0}\right)}}{{o\left(\tau_{1}\right)}}$

  * •

${\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw\right)>0\;\iff\;sx\left(x,r_{0},r_{1}\right)>\nicefrac{{o\left(\tau_{0}\right)}}{{o\left(\tau_{1}\right)}}$

Technically, to formalize this result it is convenient to use Mathlib’s
${\tt{cmp}}$555https://leanprover-
community.github.io/mathlib4_docs/Mathlib/Init/Data/Ordering/Basic.html#cmp,
which gives the order between the two parameters.

######  Lemma 3.7 (Swap rate _vs._ exchange rate).

Let $sw\in{\tt{Swap}}\left(sx,s,a,\tau_{0},\tau_{1},x\right)$ be a swap
transaction, and let $o\in{\tt{T}}\to\mathbb{R}_{>0}$ be a price oracle. For
$i\in\\{0,1\\}$, let
$r_{i}=s.{\tt{amms}}.r_{i}\left(\tau_{0},\tau_{1}\right)$. If
$s.{\tt{mints}}\left(a\right)\left(\tau_{0},\tau_{1}\right)=0$, then
cmp(gain(a, o, s, apply sw), 0) = cmp(sx(x, r_0, r_1), o(τ0)o(τ1))

######  Proof 3.8.

By algebraic manipulation and application of Lemma 3.5.
AMMLib/Transaction/Swap/Networth.lean:86

Lemma 3.9 establishes that, in a constant-product AMM, there exists only one
profitable direction for a swap. Namely, if swapping $\tau_{0}$ for $\tau_{1}$
gives a positive gain, then swapping in the other direction (i.e., $\tau_{1}$
for $\tau_{0}$) will give a negative gain. Note that the inverse does not
hold: a negative gain in a direction does not imply a positive gain in the
other direction.

######  Lemma 3.9 (Unique direction for swap gain).

Let $sw\in{\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{0},\tau_{1},x\right)$
and
$sw^{\prime}\in{\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{1},\tau_{0},x^{\prime}\right)$
be two swap transactions in opposite directions, and let
$o\in{\tt{T}}\to\mathbb{R}_{>0}$. If
${\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw\right)>0$, then
${\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw\right)<0$.

######  Proof 3.10.

By Lemma 3.7. AMMLib/Transaction/Swap/Constprod.lean:160

######  Example 3.11.

Consider a blockchain state $s$ and atomic tokens $\tau_{0}$ and $\tau_{1}$,
and assume that:  1 $a$is a trader with no minted tokens, i.e.
$(s{\tt{.mints}}\;a)\;\tau_{0}\;\tau_{1}=0$;  2 the AMM pair for
$(\tau_{0},\tau_{1})$ has been initialized in $s$ and has reserves
$r_{0}=(s{\tt{.amms}}\;\tau_{0}\;\tau_{1})=18$ and
$r_{1}=(s{\tt{.amms}}\;\tau_{1}\;\tau_{0})=6$;  3 all the AMMs in $s$ use the
constant-product swap rate function;  4 $o$is a price oracle such that
$o\;\tau_{0}=3$ and $o\;\tau_{1}=4$.  Assume that $a$ wants to sell 6 units of
$\tau_{1}$ for some units of $\tau_{0}$ with the swap transaction
$sw\in{\tt{Swap}}({\tt{constprod}},s,a,\tau_{1},\tau_{0},6)$. Then, by Lemma
3.5, $a$’s gain is given by $9\cdot 3-24=3>0$. Coherently with Lemma 3.9, any
swap in the opposite direction would give a negative gain: e.g., if $a$ sells
6 units of $\tau_{0}$, her gain would be $3/2\cdot 4-18=-12$.

We say that a swap transaction
$sw\in{\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{0},\tau_{1},x\right)$ is
_optimal_ for a given price oracle $o$ when, for all
$sw^{\prime}\in{\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{0},\tau_{1},x^{\prime}\right)$
with $x^{\prime}\neq x$ we have

${\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw^{\prime}\right)<{\tt{gain}}\left(a,o,s,{\tt{apply}}\;sw\right)$

Theorem 3.12 gives a sufficient condition for the optimality of swaps in
constant-product AMMs: it suffices that the ratio of the AMM’s reserves after
the swap is equal to the exchange rate given by the price oracle. Intuitively,
the condition in Theorem 3.12 means that the exchange rate between the two
token types induced by the AMM (i.e., the ration between the token reserves)
is aligned with the exchange rate given by the price oracle, and so further
swaps would yield a negative gain. Note that, by definition, if a swap is
optimal then it is also unique, i.e. swapping any other amount would yield a
suboptimal gain.

######  Theorem 3.12 (Sufficient condition for optimal swaps).

Let $sw\in{\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{0},\tau_{1},x\right)$ be
a swap transaction on a constant-product AMM, and let
$o\in{\tt{T}}\to\mathbb{R}_{>0}$ be a price oracle. For $i\in\\{0,1\\}$, let
$r_{i}^{\prime}=\left({\tt{apply}}\;sw\right){\tt{.amms.r}}_{i}\left(\tau_{0},\tau_{1}\right)$
be the AMM reserves after the swap. If
$\nicefrac{{r_{1}^{\prime}}}{{r_{0}^{\prime}}}=\nicefrac{{o(\tau_{0})}}{{o(\tau_{1})}}$,
then $sw$ is optimal.

######  Proof 3.13.

By cases on $x<x^{\prime}$ and by application of Lemma 3.7.
AMMLib/Transaction/Swap/Constprod.lean:184

######  Example 3.14.

Under the assumptions of Example 3.14, the exchange rate between $\tau_{1}$
and $\tau_{0}$ given by the price oracle is $(o\;\tau_{0})/(o\;\tau_{1})=3/4$,
while the exchange rate induced by the AMM (i.e., the ratio between the
reserves) is $r_{1}/r_{0}=1/3$. Hence, to satisfy the equality in Theorem 3.12
we must perform a swap that increases $r_{1}$ and decreases $r_{0}$, i.e. $a$
must sell units of $\tau_{1}$ to buy units of $\tau_{0}$, coherently with the
necessary condition for a positive gain in Example 3.14. Theorem 3.15 below
will establish exactly how many units of $\tau_{1}$ must be traded.

Theorem 3.15 gives an explicit formula for the input amount $x$ that yields an
optimal swap transaction for a given AMM pair, under the assumption that the
user firing the transaction does not hold the AMM’s minted token type. The
other implicit assumption is that the user $a$ firing the swap has the needed
amount $x$ of units of the sold token, i.e.
$x\leq(s{\tt{.atoms}}\;a)\;\tau_{1}$. In practice, this assumption can always
be satisfied with _flash loans_ , which allow $a$ to borrow the amount $x$,
perform the swap, and then return the loan in a single, atomic transaction.

######  Theorem 3.15 (Arbitrage for constant-product AMMs).

Let $sw\in{\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{0},\tau_{1},x\right)$ be
a swap transaction on a constant-product AMM, and let
$o\in{\tt{T}}\to\mathbb{R}_{>0}$ be a price oracle. For $i\in\\{0,1\\}$, let
$r_{i}=s{\tt{.amms.r}}_{i}\left(\tau_{0},\tau_{1}\right)$ be the AMM reserves.
If $a$ has no minted tokens (i.e.,
$s{\tt{.mints}}\left(a\right)\left(\tau_{0},\tau_{1}\right)=0)$ then $sw$ is
optimal if the amount of traded units of $\tau_{0}$ is:

$x\;=\;\sqrt{\frac{o(\tau_{1})\cdot r_{0}\cdot r_{1}}{o(\tau_{0})}}-r_{0}$

######  Proof 3.16.

By algebraic manipulation and Theorem 3.12.
AMMLib/Transaction/Swap/Constprod.lean:316

######  Example 3.17.

Under the assumptions of Example 3.14, we know that to perform a profitable
swap the trader must sell units of $\tau_{1}$, i.e. fire a transaction
${\tt{Swap}}\left({\tt{constprod}},s,a,\tau_{1},\tau_{0},x\right)$. Theorem
3.15 gives the optimal input value $x$, i.e., the number of sold units of
$\tau_{1}$:

$x\;=\;\sqrt{\frac{3\cdot 6\cdot 18}{4}}\,-\,6\;=\;3$

Then, the output amount is given by $sw{\tt{.y}}=3\cdot 18/(6+3)=6$ and, by
Lemma 3.5, the gain of $a$ is $6\cdot 3-3\cdot 4=6$, which maximizes it. Note
that, in the new state, the exchange rate given by the AMM coincides with that
given by the price oracle:
$(r_{1}+x)(r_{0}-sw{\tt{.y)}}=(o\;\tau_{0})/(o\;\tau_{1})$.

##  4 Related work

The closest work to ours is [12], which proposes a methodology for developing
and verifying AMMs in the Coq proof assistant. In this approach, an AMM is
decomposed into multiple interacting smart contract: e.g., each minted token
is modelled as a single smart contract, following way fungible tokens are
encoded in blockchains platforms that do not provide custom tokens natively,
as e.g. in Ethereum and Tezos. These smart contracts are then implemented as
Coq functions on top of ConCert [6], a generic model of blockchain platforms
and smart contracts mechanized in Coq. Concert is used to verify behavioural
properties of smart contracts, either in isolation or composed with other
smart contracts: this is fundamental for [12], where AMMs are specified as
compositions of multiple contracts.

More specifically, [12] applies the proposed methodology to the Dexter2
protocol, which implements a constant-product AMM based on Uniswap v1 on the
Tezos blockchain666https://gitlab.com/dexter2tz/dexter2tz/-/tree/master/. The
main properties of AMMs proved in [12] are correspondences between the state
of AMMs and the sequences of transactions executed on the blockchain. In
particular, they prove that:

  * •

the balance recorded in the main AMM contract is coherent with the actual
balance resulting from the execution of the sequence of transactions;

  * •

the supply of the minted token recorded in the main AMM contract is equal to
the actual supply resulting from the execution;

  * •

the state of the minted token contract is coherent with the execution.

Furthermore, starting from the Coq specification of the Dexter2 AMM, [12]
extracts verified CamlLigo code, which is directly deployable on the Tezos
blockchain.

Although both our work and [12] involve AMM formalizations within a proof
assistant, the ultimate goals are quite different. The formalization in [12]
closely follows the concrete implementation of a particular AMM instance
(Dexter2) and produces a deployable implementation that is provably coherent
with the proposed Coq specification. By contrast, we start from a more
abstract specification of AMM, with the goal of studying the properties that
must be satisfied by any implementation coherent with the specification. An
advantage of our approach over [12] is that it provides a suitable level of
abstraction where proving properties about the economic mechanisms of AMMs,
i.e. properties about the gain of users and of equilibria among users’
strategies. In particular, Theorem 3.15 establishes a paradigmatic property of
AMMs, which explains the economic mechanism underlying their design.

Besides these main differences, our AMM model and that in [12] have several
differences. A notable difference is that our model is based on Uniswap v2,
while [12] is based on Uniswap v1. In particular, this means that our AMMs can
handle arbitrary token pairs, while in Dexter2 any AMM pairs a token with the
blockchain native crypto-currency.

At the best of our knowledge, [12] and ours are currently the only mechanized
formalizations of AMMs in a proof assistant. Many other works study economic
properties of AMMs that go beyond those proved in this paper. The works [5]
and [4] study, respectively, an AMM model based on Uniswap similar to ours,
and a generalization of the model where the AMM is parameterised over a
_trading function_ of the AMM reserves, which must remain constant before and
after any swap transactions (in Uniswap v1 and v2, the trading function is
just the product between the AMM reserves). The work [7] studies another
generalization of Uniswap v2, where the relation between input and output
tokens of swap transactions is determined by an arbitrary _swap rate function_
, studying the properties of this function that give rise to a sound economic
mechanism of AMMs. While both [4] and [7] share the common goal of providing
general models of AMMs wherein to study their economic behaviour, they largely
diverge on the formalization: [4] is based on concepts related to convex
optimization problems, while [7] borrows formalization and reasoning
techniques from concurrency theory. The Lean 4 model proposed in this paper
follows the formalization in [7], which has the advantage of requiring far
less mathematical dependencies: although Mathlib is equipped to reason about
convex sets and functions777https://leanprover-
community.github.io/mathlib4_docs/Mathlib/Analysis/Convex/Function, it is
currently lacking in advanced convex optimization definitions and results as
those used in [4].

Our AMM model is based on Uniswap v2, one of the most successful AMMs so far.
We briefly discuss some alternative AMM constructions. Balancer [2]
generalizes the constant-product function used by Uniswap to a constant
(weighted geometric) mean
$f(r_{1},\cdots,r_{n})=\prod_{i=1}^{n}r_{i}^{w_{i}}$, where the weight $w_{i}$
reflects the relevance of a token $\tau_{i}$ in a tuple
$(\tau_{1},\cdots,\tau_{n})$. Curve [9] mixes constant-sum and constant-
product functions, aiming at a swap rate with small fluctuations for large
amounts of swapped tokens. The work [11] studies a variant of the constant-
product swap rate invariant, where the rate adjusts dynamically based on
oracle prices, with the goal of reducing the need for arbitrage transactions.
Other approaches aiming at the same goal are studies in [1, 11], and
implemented in [3]. Extending our Lean formalization and results to these
alternative AMM designs would require a substantial reworking of our model and
proofs.

##  5 Conclusions

In this work we have provided a formalization of AMMs in Lean. Blockchain
states are represented as structures containing wallets and AMMs, and
transactions as dependent types equipped with a function that defines the
state resulting from firing the transaction. Based on this, we have modeled
the key economic notions of price, networth and gain. We have then focused on
the economic properties of AMMs, constructing machine-checkable proofs. In
Lemma 3.5 we have given an explicit formula for the economic gain of a user
after firing a swap transaction. In Theorem 3.12 we have proved that the
rational strategy for traders leads to the alignment between the AMM internal
exchange rate and that given by price oracles. Finally, in Theorem 3.15 we
have derived the amount of tokens that a trader should sell to maximize the
gain from a constant-product AMM.

### Design choices

Before coming up with our Lean formalization, we have experimented with a few
alternative definitions. Currently, we use the two pairs of atomic token types
$\left(\tau_{0},\tau_{1}\right)$ and $\left(\tau_{1},\tau_{0}\right)$ to
represent the same minted token type. Initially, we used the type $M$ of sets
of atomic tokens of cardinality 2, and modeled wallets of minted tokens by the
type $M\to_{0}\mathbb{R}_{\geq 0}$. However, using $M$ in the definition of
AMMs turned out not to be as easy. The type $M\to_{0}\mathbb{R}_{\geq 0}^{2}$
would obviously not work since we would not know which value corresponds to
the reserve of which token. On the other hand, the dependent type
$\left(m:M\right)\to_{0}{\tt{Option}}\left(m\to\mathbb{R}_{>0}\right)$ would
work after defining $0:={\tt{None}}$, at the cost of losing the
straightforward definition for supply of atomic tokens. We have opted for the
custom subtype $\mathbb{R}_{>0}$ to represent the positive reals, since they
simplify writing certain definitions and proof passages (e.g., avoiding the
use of garbage outputs in the definitions where negative inputs would not make
sense). This choice however turned out to have some cons, since it makes using
Mathlib more complex, and in some cases we have to coerce back to the reals
anyway (e.g. when reasoning about the gain). Using Mathlib’s reals would
perhaps lead to a smoother treatment.

### Limitations

Compared to real-world AMM implementations, our Lean formalization introduces
a few simplifications, that overall contribute to keeping our proofs
manageable. Bridging the gap with real AMMs would require several extensions,
which we discuss below as directions for future work. AMMs typically implement
a trading fee $\phi\in[0,1]$ that represents the portion of the swap amount
kept by the AMM. While modeling the fee would be easy ($\phi$ would be an
additional parameter to ${\tt{SX}}$ functions, to ${\tt{Swap}}$ types, and to
transactions ${\tt{Tx}}$), it would require a major reworking of all the
results that deal with swaps. Our swaps have _zero-slippage_ , in that either
a swap gives exactly the amount of tokens required by a user, or they are
aborted. While on the one hand this is desirable (e.g., it rules out sandwich
attacks), on the other hand it has drawbacks related to liveness, since a user
may need to repeatedly send swap transactions until one is accepted. Real-
world AMM implementations allow users to specify a slippage tolerance in the
form of the minimum amount of tokens they expect from a swap. Extending our
model to encompass slippage tolerance would require to add a parameter to each
transaction type, and a minor reworking of the results. Some AMM
implementations allow users to create AMMs pairs involving _minted_ token
types. Consequently, the tokens minted by these AMMs in general are “nestings”
of token types. Extending our model in this direction would require to replace
price oracles in our results with a suitable price function.

## References

  * [1] Improving frontrunning resistance of x*y=k market makers, 2018.  https://ethresear.ch/t/improving-front-running-resistance-of-x-y-k-market-makers/1281.
  * [2] Balancer whitepaper, 2019.  https://balancer.finance/whitepaper/.
  * [3] Mooniswap whitepaper, 2020.  https://mooniswap.exchange/docs/MooniswapWhitePaper-v1.0.pdf.
  * [4] Guillermo Angeris and Tarun Chitra.  Improved price oracles: Constant function market makers.  In ACM Conference on Advances in Financial Technologies (AFT), pages 80–91. ACM, 2020.  https://arxiv.org/abs/2003.10001.  doi:10.1145/3419614.3423251.
  * [5] Guillermo Angeris, Hsien-Tang Kao, Rei Chiang, Charlie Noyes, and Tarun Chitra.  An analysis of Uniswap markets.  Cryptoeconomic Systems, 1(1), 2021.  doi:10.21428/58320208.c9738e64.
  * [6] Danil Annenkov, Jakob Botsch Nielsen, and Bas Spitters.  Concert: a smart contract certification framework in Coq.  In ACM SIGPLAN International Conference on Certified Programs and Proofs (CPP), pages 215–228. ACM, 2020.  doi:10.1145/3372885.3373829.
  * [7] Massimo Bartoletti, James Hsin yu Chiang, and Alberto Lluch-Lafuente.  A theory of automated market makers in DeFi.  Logical Methods in Computer Science, Volume 18, Issue 4, dec 2022.  URL: https://doi.org/10.46298%2Flmcs-18%284%3A12%292022, doi:10.46298/lmcs-18(4:12)2022.
  * [8] Certora.  Formal verification of Compound’s open-oracle with Uniswap anchor.  https://files.safe.de.fi/safe/files/audit/pdf/CompoundUniswapAnchoredOpenOracleAug2020.pdf, 2020.
  * [9] Michael Egorov.  Stableswap - efficient mechanism for stablecoin, 2019.  https://curve.fi/files/stableswap-paper.pdf.
  * [10] Uri Kirstein.  Detecting corner cases in Compound V3 with formal specifications.  https://medium.com/certora/detecting-corner-cases-in-compound-v3-with-formal-specifications-b7abf137fb15, 2022.
  * [11] Bhaskar Krishnamachari, Qi Feng, and Eugenio Grippo.  Dynamic curves for decentralized autonomous cryptocurrency exchanges.  In International Symposium on Foundations and Applications of Blockchain (FAB), volume 92 of OASIcs, pages 5:1–5:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.  doi:10.4230/OASIcs.FAB.2021.5.
  * [12] Eske Hoy Nielsen, Danil Annenkov, and Bas Spitters.  Formalising decentralised exchanges in Coq.  In ACM SIGPLAN International Conference on Certified Programs and Proofs (CPP), pages 290–302. ACM, 2023.  doi:10.1145/3573105.3575685.
  * [13] Jiahua Xu, Krzysztof Paruch, Simon Cousaert, and Yebo Feng.  Sok: Decentralized exchanges (DEX) with automated market maker (AMM) protocols.  ACM Comput. Surv., 55(11):238:1–238:50, 2023.  doi:10.1145/3570639.


# Shack-Hartmann wavefront sensing: A new approach to time-resolved
measurement of stress intensity during dynamic fracture of small brittle
specimens

Liuchi Li<EMAIL_ADDRESS>Velat Kilic  Milad Alemohammad  K.T. Ramesh  Mark A.
Foster  Todd C. Hufnagel<EMAIL_ADDRESS>Hopkins Extreme Materials
Institute, Johns Hopkins University, Baltimore, MD 21218, USA  Department of
Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD
21218, USA  Department of Materials Science and Engineering, Johns Hopkins
University, Baltimore, MD 21218, USA  Department of Mechanical Engineering,
Johns Hopkins University, Baltimore, MD 21218, USA

###### Abstract

The stress intensity factor describes the stress state around a crack tip in a
solid material and is important for understanding crack initiation and
propagation. Because stresses cannot be measured directly, the
characterization of the stress intensity factor relies on the measurement of
deformation around a crack tip. Such measurements are challenging for dynamic
fracture of brittle materials where the deformation is small and the crack tip
velocity can be high ($>$1\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$).
Digital gradient sensing (DGS) is capable of full-field measurement of surface
deformation with sub-microsecond temporal resolution, but it is limited to
centimeter-scale specimens and has a spatial resolution of only
$\sim$1\text{\,}\mathrm{m}\mathrm{m}$$. This limits its ability to measure
deformations close to the crack tip. Here, we demonstrate the potential of
Shack-Hartmann wavefront sensing (SHWFS), as an alternative to DGS, for
measuring surface deformation during dynamic brittle fracture of millimeter-
scale specimens. Using an opaque commercial glass ceramic as an example
material, we demonstrate the capability of SHWFS to measure the surface slope
evolution induced by a propagating crack on millimeter-scale specimens with a
micrometer-scale spatial resolution and a sub-microsecond temporal resolution.
The SHWFS apparatus has the additional advantage of being physically more
compact than a typical DGS apparatus. We verify our SHWFS measurements by
comparing them with analytical predictions and phase-field simulations of the
surface slope around a crack tip. Then, fitting the surface slope measurements
to the asymptotic crack-tip field solution from linear elastic fracture
mechanics, we extract the evolution of the apparent stress intensity factor
associated with the propagating crack tip. We conclude by discussing potential
future enhancements of this technique and how its compactness could enable the
integration of SHWFS with other characterization techniques including x-ray
phase-contrast imaging (XPCI) for multi-modal characterization of dynamic
fracture.

###### keywords:

Dynamic fracture, Brittle material, Shack-Hartmann wavefront sensor, Surface
slope measurement, Stress intensity factor, Phase-field simulation

##  1 Introduction

The stress intensity factor (SIF or $K_{\mathrm{I}}$ in Mode-I condition) is a
critical parameter for understanding the fracture of materials. Although it
cannot be measured directly, the stress intensity factor can be obtained
indirectly from measurements of displacements around a crack tip. An ideal
measurement technique would provide full-field displacement measurement with
good spatial resolution, be fast enough to allow measurement of the stress
intensity around the tip of a propagating crack, and be suitable for use on an
arbitrary material. Particularly challenging are studies of brittle materials,
where the in-plane displacements around the crack tip can be small (on the
order of tens of nanometers) and the crack-tip speed is high
($>$1\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$$).

Over the past fifty years, extensive research has resulted in the development
of a series of increasingly powerful optical techniques for characterizing the
stress intensity of propagating cracks in a variety of materials. Notable
examples include measurements relying on the photoelastic effect [1], the
method of caustics [2, 3, 4], coherent gradient sensing (CGS) [5, 6, 7, 8]),
digital image correlation (DIC) [9, 10]), and, most recently, digital gradient
sensing (DGS) [11, 12, 13]. Any of these techniques may be effective,
depending on the material under examination, the specimen geometry, and the
loading conditions, but all have limitations [5, 11].

Measuring the stress intensity during dynamic fracture of brittle materials
(defined here as materials showing little or no plastic deformation under the
loading conditions of interest) is especially challenging because the
displacements around the crack tip are small, the time available for the
measurement is limited by the high crack tip speeds (often in excess of
$1\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$), and there is the potential
for crack branching [14]. One of the most powerful techniques is digital
gradient sensing (DGS), which can provide a full-field measurement around a
dynamically propagating crack tip for both transparent and opaque brittle
materials [11, 14]. DGS has been applied to the study of dynamic fracture of a
variety of materials including polymers [15], oxide glasses [16], and polymer-
matrix composites [17, 18].

When performed in reflection from an opaque material (or transparent material
with a reflective coating), DGS measures the deflection of light rays due to
the surface profile and, in particular, the slope associated with out-of-plane
displacements due to stress concentration around the crack tip. From these
displacements, the stress intensity factor can be extracted using fracture
mechanics theory and the instant crack velocity [19]. DGS has an advantage
over DIC for measurements on brittle materials, in that a long optical path
(several meters) amplifies the small displacements enabling higher
sensitivity. The field of view of DGS measurements is typically on the order
of a few centimeters, with a spatial resolution of around a millimeter [14,
16, 20]. These characteristics make it challenging to use DGS on small
samples, and require a significant amount of laboratory space for
implementation.

Here we propose a new approach (as an alternative to DGS) to measure surface
profile gradients during dynamic fracture of brittle materials, based on
Shack-Hartman wavefront sensing (SHWFS) [21, 22]. Our technique retains the
advantages of gradient-based methods but provides significantly improved
spatial resolution
($\sim$10\text{\,}\mathrm{\SIUnitSymbolMicro}\,\mathrm{m}$$) while making the
experimental setup more compact (with a total optical path of
$<$1\text{\,}\mathrm{m}$$). This combination makes SHWFS ideal for multi-modal
studies, particularly in combination with x-ray phase contrast imaging (XPCI)
performed at synchrotron facilities [23, 24, 25]. In this paper, we describe
the SHWFS technique, apparatus, and data analysis; illustrate its application
to the study of dynamic fracture of a commercial glass ceramic (including
extraction of the apparent stress intensity factor); and discuss potential
extensions and enhancements of the technique.

##  2 Shack-Hartmann wavefront sensor (SHWFS)

Figure 1: (a) A 2D schematic along the $x$ direction demonstrating the working
principle of SHWFS. The distorted wavefront goes through a collection of
microlenses (with a uniform focal length $f$), each of which performs a lens-
area-averaged estimation of the wavefront’s gradient, denoted as
$\bar{w}_{x}$. (b) A 3D schematic demonstrating the working principle of
SHWFS.

###  2.1 Working principle

The Shack-Hartman wavefront sensor (SHWFS) was originally developed to measure
distortions in telescope images due to atmospheric turbulence and enabled the
field of adaptive optics [22]. Here we use SHWFS in a manner similar to
reflection-mode DGS to measure surface slopes and extract fracture parameters
during dynamic fracture. The major difference from DGS is in how the surface
slope is measured. In DGS, an image of a reference pattern is recorded in
reflection off of the sample; when the sample is deformed, the pattern becomes
distorted relative to the original image, and from the distortion, the surface
gradients are determined. In SHWFS, instead of imaging a reference pattern, a
microlens array (MLA) is used to measure distortions of a planar wavefront
reflected from the sample.

A basic SHWFS consists of a micro-lens array and an image sensor, as shown in
Fig. 1. With the sensor positioned at the focal plane of the MLA and
illuminated by a plane wave, each micro-lens creates a focused spot on its
respective optical axis. A distorted wavefront on the other hand produces a
focused spot displaced from the optical axis. The displacement of the centroid
of each spot is then proportional to the orientation of the wavefront averaged
across the aperture of the corresponding microlens. The microlenses thus map
angular deformations of the wavefront to centroid shifts at the focal plane as
depicted in Fig. 1(a):

$\displaystyle\bar{w}_{x}\approx\Delta x/f,$ (2.1)

where $\bar{w}_{x}$ is the averaged gradient of the wavefront across a
particular microlens, $\Delta x$ is the displacement of the centroid on the
detector plane, and $f$ is the focal length of the microlens. Fig. 1(b)
extends Fig. 1(a) to a 3D schematic, where we can also compute similarly the
gradient in the transverse direction,

$\displaystyle\bar{w}_{y}\approx\Delta y/f.$ (2.2)

In essence, by sampling the wavefront with an MLA all of these shifts can be
measured, and the whole wavefront can be reconstructed provided the scene is
properly sampled both spatially and temporally [26], see section. 2.2 for a
more detailed discussion.

###  2.2 Application in the context of dynamic fracture

We use one surface of a sample to generate a wavefront by reflecting off
incoming collimated light rays. For this work, we consider only the Mode-I
loading condition and use the single-notched three-point bend configuration
for our experiment. As shown in Fig. 2(a), initially, when the sample is
stress-free, the surface is perfectly flat, giving a planar wavefront shown as
a set of equally spaced spots on the detector plane. Next, when the sample is
loaded in-plane, due to Poisson’s effect and heterogeneous in-plane stress
distribution, the sample surface will deform out-of-plane heterogeneously,
leading to non-zero surface slope and shifts of those spots on the detector
plane, as shown in Fig. 2(b). During a dynamic fracture event, a propagating
crack at every time instant will cause a spatial variation of the surface
slope (especially near the crack tip due to stress concentration), which we
calculate by first measuring the shifts $\Delta x$ and $\Delta y$ and next
feeding into Eqns. 2.1 and 2.2. For successful detection of $\Delta x$ and
$\Delta y$, the following requirements should be met:

Figure 2: Application of SHWFS in the context of dynamic fracture using a
single-notch three-point bending configuration. The planar (reference)
wavefront is generated by the sample’s smooth surface reflecting incoming
collimated light rays, shown in (a). In contrast, the distorted (deformed)
wavefront is caused by the surface deforming out-of-plane due to Poisson’s
effect, as shown in (b). In (a), the polar coordinate associated with the
crack tip is also shown.

  * 1.

_Spatial sampling_ — The MLA pitch (the distance from the center of one
lenslet to the other) should be smaller than the spatial correlation length of
the distorted wavefront so that the wavefront can be safely assumed to be
locally planar. Indeed, any function can be approximated locally linearly,
provided it is sampled sufficiently densely. Although a smaller pitch length
is preferred, we note that diffraction puts a lower bound on the pitch length
for an MLA to function properly. Thus, we can think of the pitch length as
similar to the correlation window size in DIC: a large window size reduces
measurement noise but is prone to average out locally high deformation
gradients, the opposite is true for a small window size. There is also a lower
bound for the correlation window size with regard to accurate measurements,
which is related to the average speckle feature size and contrast [27].

  * 2.

_Temporal sampling_ — In dynamic fracture experiments the sample surface can
deform significantly and rapidly. To measure surface distortions at a given
time instant, the incident wavefront must be essentially constant, which
requires that it be sampled at a rate faster than its correlation time. If the
incident wavefront varies significantly over a single exposure time of the
recording camera, the focused spot will move as a result, and the integrated
image will be blurred. The same principle also applies to DIC analysis (and
any other imaging techniques), where the in-plane deformation can be viewed as
the incident wavefront, whose variation should not be significant within a
single exposure time of the recording camera. Therefore, a high-speed camera
is needed. For imaging dynamic fracture in glass and ceramics, a frame rate on
the order of one million per second is sufficient [14, 16].

  * 3.

_Spot size and shape_ — A tightly focused spot array on the detector plane is
essential for accurate detection, and it ideally requires collimated
illumination from a spatially confined source (fiber-coupled laser in our
case), a smooth sample surface, and a large enough lenslet diameter. A
collimated source is necessary to achieve a tightly focused spot array since
the source plane is mapped onto the detector plane. However perfect
collimation is not a strict requirement since the spot positions are measured
differentially (i.e., referenced with respect to the positions before crack
initiation). A smooth sample surface ensures that the in-coming planar (i.e.,
collimated) illumination remains largely unchanged after reflection. In
practice, a mirror-like finish (with a roughness average
$R_{a}<$0.3\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$) is sufficient for
our purpose. Lastly, a small lenslet will lead to a large (as opposed to
tightly focused) spot shape in the Fourier plane due to diffraction, which is
undesired.

  * 4.

_Dynamic range_ — The MLA pitch and focal length place a limit on the surface
tilt that can be measured. In particular, the shift of the focused spots must
be large enough to be measured, but not so large that spots from adjacent
lenses impinge on one another. In that case, it becomes challenging to
disambiguate whether a change in the centroid position is caused by a small
deformation, or a deformation so large that the centroid is coupled to the
neighboring lenslet area. This effectively puts a limit on the range of spot
shifts that can be measured, with a lower bound set by the camera pixel size
and an upper bound set by the MLA pitch. We note, though, that computational
methods could be used to extend this range. Indeed, the spot centroids may be
located with a sub-pixel resolution if the spot is spatially sampled
sufficiently densely by the camera pixels and with a sufficient signal-to-
noise ratio (SNR) [28].

##  3 Sample preparation and experimental setup

###  3.1 Sample preparation

We use Macor™, a commercially available glass ceramic, as an example material
to test the utility of SHWFS because its mechanical properties (such as the
fracture toughness $K_{\text{IC}}$ and Young’s modulus $E$) are well
characterized and because it is highly machinable, making sample fabrication
convenient. We purchased Macor™™ plates (Master-Carr 8489K231) with an initial
thickness of $0.0625\text{\,}\mathrm{i}\mathrm{n}$
($1.58\text{\,}\mathrm{m}\mathrm{m}$). We used a diamond wire saw machine to
cut each plate into rectangular bars $12\text{\,}\mathrm{m}\mathrm{m}$ long
(along the $y$ direction), $3\text{\,}\mathrm{m}\mathrm{m}$ wide (along the
$x$ direction), and $1.58\text{\,}\mathrm{m}\mathrm{m}$ thick (denoted as $h$
in Fig. 2(a), along the $z$ direction). We next used the same saw to create a
$1\text{\,}\mathrm{m}\mathrm{m}$ deep notch in each bar to produce the single-
notched three-point bending configuration mentioned in the previous section.
The wire used has a diameter of
$250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, creating a notch width of
$250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and a semicircular notch
tip. We next polished one face (the $x$-$y$ plane) of each bar to a mirror-
like finish with diamond lapping films down to a
$0.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ grade). This forms the
reflective surface necessary for the SHWFS approach. After polishing,
specimens are found to have a reduced thickness of $1.5\pm 0.2$ mm. Note that
it is not necessary to deposit a reflective coating, at least not for these
specimens. For transparent materials such as glass, a reflective coating is
desired.

###  3.2 Experimental setup

We use a custom-designed loading apparatus [25] to induce dynamic fracture
under three-point bending. In this apparatus, an indenter is connected to a
piezoelectric actuator (Cedrat Technologies PPA40M), with the actuator being
driven by a voltage signal that is generated by a function generator
(Tektronix AFG3252) and then amplified by a high-speed voltage amplifier
(PiezoDrive PD200). We place each sample over a loading plate with rolling
support that sits above a vertical translation stage (OptoSigma TSD-653DMUU).
We use the function generator to output a single linear voltage ramp from
$0\text{\,}\mathrm{V}$ to $7\text{\,}\mathrm{V}$ over a time window of
$256\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$, resulting in an
indentation speed of approximately $0.13\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$.
(Details of the speed measurement are provided in the appendix.) This
apparatus is compact enough for operations in synchrotron facilities, and it
has been used to characterize dynamic fracture of different brittle materials
via XPCI [25, 29]. With this loading apparatus inducing crack propagation
within a sample, we use a high-speed camera (Shimadzu HPV-X2) to capture the
distorted wavefront (subsequently $\bar{w}_{x}$ and $\bar{w}_{y}$) evolution
around the propagating crack tip.

Figure 3 shows a schematic and a photograph of the optical setup. Illumination
of the sample was provided by a pulsed laser beam (SI-LUX 640) with a fiber-
coupled collimator (M92L01 and RC08FC-F01 from Thorlabs). The reflection from
the sample was imaged with an infinity-corrected microscopy system (Mitutoyo
$5\times$ objective with a ThorLabs TTL200-A tube lens), a microlens array
(ThorLabs MLA150-5C-M with a pitch of
$150\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$), and a telescopic relay
lens pair (ThorLabs MP105075-A) with $1.5\times$ magnification. Note the
compactness of the experimental setup, which fits in considerably less space
than typically required for DGS. A typical spot pattern from an undeformed
specimen is shown in Fig. 3, where the field of view is approximately
$$1.6\text{\,}\mathrm{m}\mathrm{m}$\times$1\text{\,}\mathrm{m}\mathrm{m}$$,
with the spacing between spots about
$45\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$. The spot spacing can be
calculated using the pitch of the microlens array and the magnification ratio
of the microscope and the relay:
$$150\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$/5\times
1.5=$45\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$.

Figure 3: (a) Schematic of the SHWFS design. (b) Top-down photograph of the
SHWFS setup. (c) Image formed on the camera detector showing the planar
wavefront generated from the polished surface of the undeformed sample.

##  4 Calculation of shift and surface slope from SHWFS recording

###  4.1 Shift and surface slope

As a crack initiates and propagates through a sample, it will cause spots
focused on the image plane to shift laterally. (Refer to the Supplemental
Material for a representative recording.) We measure these shifts and later
use them to extract the apparent stress intensity factor (SIF). For a
particular recording, we start from the initial reference image and perform a
spatial partition such that each focused spot belongs to one and only one
cell, as shown in Fig. 4(a). For images captured after crack initiation, we
use this partition to calculate the shift $\Delta x$ and $\Delta y$ of each
spot by finding its centroid in the reference (denoting as $[x_{0},y_{0}]$)
and the deformed configuration at a particular time instance $t$ (denoting as
$[x_{t},y_{t}]$):

$\displaystyle[\Delta x,\Delta y]|_{t}=[x_{t}-x_{0},y_{t}-y_{0}].$ (4.1)

The procedure for finding the centroid of each spot is illustrated in Fig. 4.
For a given spot we calculate the intensity-weighted average of all pixels
that constitute that spot, in a way similar to Ref. [30]. Taking one spot in
the reference image as an example, its centroid is given by

$\displaystyle x_{0}=\frac{\sum_{i,j\in S}I_{ij}\times i}{\sum_{i,j\in
S}I_{ij}},y_{0}=\frac{\sum_{i,j\in S}I_{ij}\times j}{\sum_{i,j\in S}I_{ij}},$
(4.2)

where $(i,j)$ denotes the location of a pixel, $S$ denotes the cell captured
by an individual lenslet, and $I_{ij}$ denotes the intensity of the pixel at
location $(i,j)$. We then apply Eqn. 4.2 to find $(x_{t},y_{t})$ for the same
spot in each subsequent image (Fig. 4(b)). Lastly, we calculate $\Delta x$ and
$\Delta y$ using Eqn. 4.1, as shown in Fig. 4(c).

Figure 4: (a) Spatial partitioning of light spots, each assigned to a unique
square cell. (b) Examples of calculating the centroid of one spot in the
reference (top) and deformed (bottom) configurations. (c) An example of
calculating the shifts $\delta x$ and $\delta y$ by overlaying the two images
from (b). Areas colored in either green or purple are mismatches between the
two images.

The next step in the data reduction is converting the spot shifts ($\Delta x$
and $\Delta y$) to surface slopes ($\bar{w}_{x}$ and $\bar{w}_{y})$. Care must
be taken in this process because the microscope provides transverse
magnification for increased spatial resolution, but this leads to angular de-
magnification which in turn reduces the sensitivity to out-of-plane
deformations. In Fourier optics terms, magnifying a 2D signal must necessarily
reduce its spatial bandwidth. On the other hand, the relay used after the MLA
images the focal plane of the MLA and increases system sensitivity to out-of-
plane deformations. Putting these effects together,

$\displaystyle\frac{s_{r}}{s_{m}}\bar{w}_{x}=\frac{d\,\Delta
x}{f}\qquad\text{and}\qquad\frac{s_{r}}{s_{m}}\bar{w}_{y}=\frac{d\,\Delta
y}{f},$ (4.3)

where $s_{m}$ and $s_{r}$ are the magnification factors of the microscope
objective ($s_{m}=5$ here) and the relay lens pair ($s_{r}=1.5$ here), $d$ is
the pixel size of the camera sensor ($d=30\,\,\mu$m), $\Delta x$ (or $\Delta
y$) is the shift in pixels, and $f$ is the focal length of the MLA
($f=$5.6\text{\,}\mathrm{m}\mathrm{m}$$ here). In contrast, if we want to
calculate in-plane information such as the distance $L$ between any two pixels
$(i_{1},j_{1})$ and $(i_{2},j_{2})$, we will have:

$\displaystyle s_{r}s_{m}L=d\sqrt{(i_{1}-i_{2})^{2}+(j_{1}-j_{2})^{2}}.$ (4.4)

###  4.2 Extraction procedure for $K_{\text{I}}$ using slope data

Once we have measured the shifts $\Delta x$ and $\Delta y$ of each spot and
computed the corresponding surface slope $\bar{w}_{x}$ and $\bar{w}_{y}$ using
Eqn. 4.3, we can extract the apparent stress intensity factor $K_{\text{I}}$
by fitting the slope data to the asymptotic crack-tip stress field [31, 32]
for a Mode-I crack opening under plane stress condition:

$\displaystyle\begin{split}-\frac{2E}{\nu
h}\bar{w}_{x}=\frac{\partial\left(\hat{\sigma}_{xx}+\hat{\sigma}_{yy}\right)}{\partial
x}=-\frac{1}{2}r^{-\frac{3}{2}}g(V)A_{1}\text{cos}\left(\frac{3}{2}\theta\right)+\sum_{N=2}^{\infty}A_{N}\left(\frac{N}{2}-1\right)r^{\left(\frac{N}{2}-2\right)}\text{cos}\left[\left(\frac{N}{2}-2\right)\theta\right],\\\
-\frac{2E}{\nu
h}\bar{w}_{y}=\frac{\partial\left(\hat{\sigma}_{xx}+\hat{\sigma}_{yy}\right)}{\partial
y}=-\frac{1}{2}r^{-\frac{3}{2}}g(V)A_{1}\text{sin}\left(\frac{3}{2}\theta\right)+\sum_{N=2}^{\infty}A_{N}\left(\frac{N}{2}-1\right)r^{\left(\frac{N}{2}-2\right)}\text{sin}\left[\left(\frac{N}{2}-2\right)\theta\right],\end{split}$
(4.5)

where $\nu,E,h$ are the Poisson’s ratio, the Young’s modulus, and the
thickness of the sample, respectively; $\hat{\sigma}_{xx}$ and
$\hat{\sigma}_{yy}$ are the thickness averages of stress components along the
$x$ and $y$ directions; $(r,\theta)$ denotes the position relative to the
crack-tip in polar coordinates (see Fig. 2(a)), where $r$ at each data point
can be computed following Eqn. 4.4; $A_{1}=K_{\text{I}}\sqrt{\frac{2}{\pi}}$
where $K_{\text{I}}$ is the Mode-I stress intensity factor (SIF); $V$ is the
crack tip speed with $g(V)$ being a function of the instantaneous crack tip
speed to account for the velocity dependency of the K-field [32]. It takes the
following expression:

$\displaystyle
g(V)=\frac{1+\nu}{1-\nu}\frac{\left(1-\alpha_{L}^{2}\right)\left(1+\alpha_{S}^{2}\right)}{4\alpha_{L}\alpha_{S}-\left(1+\alpha_{S}^{2}\right)^{2}},\,\,\text{with}\,\,\alpha_{L}^{2}=1-\frac{V^{2}}{C_{L}^{2}}\,\,\text{and}\,\,\alpha_{S}^{2}=1-\frac{V^{2}}{C_{S}^{2}},$
(4.6)

where $C_{L}$ and $C_{S}$ are longitudinal and shear wave speeds in the
material. For relatively slow cracks ($V<0.4C_{R}$ where $C_{R}$ is the
Rayleigh wave speed [33]), $g(V)\simeq 1$ and the velocity dependence is
negligible [32]. Because our experiments are conducted within this slow crack
velocity regime
($V\sim$1\text{\times}{10}^{2}\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$$ while
$C_{R}>$1\text{\times}{10}^{3}\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$$), we take
$g(V)=1$ and neglect higher-order terms in Eqn. 4.5, resulting in the
following expressions:

$\displaystyle\begin{split}-\frac{2E}{\nu
h}\bar{w}_{x}=\frac{\partial\left(\hat{\sigma}_{xx}+\hat{\sigma}_{yy}\right)}{\partial
x}\simeq-\frac{K_{\text{I}}}{\sqrt{2\pi}}r^{-\frac{3}{2}}\text{cos}\left(\frac{3}{2}\theta\right),\\\
-\frac{2E}{\nu
h}\bar{w}_{y}=\frac{\partial\left(\hat{\sigma}_{xx}+\hat{\sigma}_{yy}\right)}{\partial
y}\simeq-\frac{K_{\text{I}}}{\sqrt{2\pi}}r^{-\frac{3}{2}}\text{sin}\left(\frac{3}{2}\theta\right).\end{split}$
(4.7)

With these expressions in hand, we can determine the stress intensity
($K_{\text{I}}$) during crack propagation from the surface slope data using an
iterative procedure. Our approach to data selection and fitting is similar to
that presented in Ref. [16], and has three basic steps:

  * 1.

Identify the crack tip location P based on the spatial distribution of
$\bar{w}_{x}$ and $\bar{w}_{y}$. We use a parameter $r_{c}$ to assess the
uncertainty associated with P, as shown in Figs. 5(a) and (b);

  * 2.

Choose a subset of the $\bar{w}_{x}$ and $\bar{w}_{y}$ data to use in
determining $K_{\text{I}}$. The choice is described by four parameters:
$r_{\text{min}}$, $r_{\text{max}}$, $\theta_{x}$, and $\theta_{y}$. The points
chosen for $\bar{w}_{y}$ (Fig. 5(a)) correspond to $r_{\text{min}}\leq r\leq
r_{\text{max}}$ and either $\theta_{y}\leq\theta\leq\pi-\theta_{y}$ or
$\pi+\theta_{y}\leq\theta\leq 2\pi-\theta_{y}$. The points chosen for
$\bar{w}_{x}$ (Fig. 5(b)) are similar but with a different range of angles
(specified by $\theta_{x}$ rather than $\theta_{y}$). Note that the number of
points chosen for $\bar{w}_{y}$ and $\bar{w}_{x}$ may be different in general.

  * 3.

Fit the data from the selected points using Eqn. 4.5 to determine
$K_{\text{I}}$. To do so, we first calculate the coefficients (denoted as
$b_{1}$ and $b_{3}$) associated with the first two nonzero terms in Eqn. 4.5
(_i.e._ those with $A_{1}$ with and $A_{3}$) for each data point selected from
the previous step. $b_{1}$ and $b_{3}$ are determined from the geometrical
location ($r$ and $\theta$) of each selected point. Once we have these
coefficients, we perform a least-squares minimization to find the vector
$\mathbf{A}=[A_{1},A_{3}]^{T}$ (subsequently extracting
$K_{\text{I}}=\sqrt{\pi/2}A_{1}$), whose solution is given by the pseudo
inverse of $\mathbf{B}$ shown in Fig. 4.5(c):

$\displaystyle\mathbf{A}^{*}=\min_{\mathbf{A}}\left\lVert\mathbf{B}^{(y)}\mathbf{A}-\mathbf{W}^{(y)}\right\lVert\quad\text{or}\quad\mathbf{A}^{*}=\min_{\mathbf{A}}\left\lVert\mathbf{B}^{(x)}\mathbf{A}-\mathbf{W}^{(x)}\right\lVert,$
(4.8)

where $\mathbf{B}^{(y)}$ (or $\mathbf{B}^{(x)}$) is a matrix containing
$b_{1}$ and $b_{3}$, and $\mathbf{W}^{(y)}$ (or $\mathbf{W}^{(x)}$) is a
vector containing the surface slope data for the $y$ direction (or the $x$
direction), as shown in Fig. 5(c). Alternatively, we can merge the two
minimization problems into one:
$\mathbf{A}^{*}=\min_{\mathbf{A}}\left\lVert[\mathbf{B}^{(y)};\mathbf{B}^{(x)}]\mathbf{A}-[\mathbf{W}^{(y)};\mathbf{W}^{(x)}]\right\lVert$.
This fitting procedure using the out-of-plane information is similar to that
which uses in-plane information (such as from DIC measurements) to extract
$K_{\text{I}}$ for polymers [10].

Figure 5: Selection of points for determining $K_{\mathrm{I}}$ from surface
slope data. Parts (a) and (b) show contour plots of $\bar{w}_{y}$ and
$\bar{w}_{x}$, respectively (see 5.1 for a description of these data). The
superimposed circles represent positions from which the surface slope data
were used to determine $K_{\mathrm{I}}$. The region of data used is specified
by four parameters: the position of the crack tip, $\mathbf{P}$; inner and
outer radii, $r_{\text{min}}$ and $r_{\text{max}}$; and the half-angle of an
excluded region, $\theta_{y}$ and $\theta_{x}$. The parameter $r_{c}$ reflects
uncertainty in determining $\mathbf{P}$, as described in the text. (c) A flow
chart showing steps for calculating the apparent stress intensity factor
$K_{\text{I}}$.

As we can see from Eqn. 4.5 the terms dominating the solution of Eqn. 4.8 are
those from points close to the crack tip, which have large $r^{-3/2}$ and
$r^{-1/2}$ as well as large $\bar{w}_{x}$ and $\bar{w}_{y}$. This makes
precise determination of the position of the crack tip ($\mathbf{P}$)
critical. On the other hand, we want to exclude data that are too close to the
crack tip and near the crack path, because the SHWFS measurements become
unreliable there due to the highly localized deformation and discontinuity.
Thus, selecting an appropriate value of $r_{\text{min}}$ balancing these
considerations is also important. The other parameters ($r_{\text{max}}$,
$\theta_{x}$, and $\theta_{y}$) are much less important than $\mathbf{P}$ and
$r_{\text{min}}$.

In practice, we have found that choosing $r_{\text{min}}=20$ pixels
($0.08\text{\,}\mathrm{m}\mathrm{m}$ on the specimen), which excludes data
from the closest two or three spots to the crack tip, works well with our
setup. Other parameters chosen are $r_{\text{max}}$= 110 pixels (equivalent to
$0.4\text{\,}\mathrm{m}\mathrm{m}$ on the sample, which covers essentially the
whole horizontal field of view), $\theta_{y}=\pi/6$, and $\theta_{x}=\pi/4$.
Lastly, we locate $\mathbf{P}$ for each snapshot based on the symmetric
spatial pattern of $w_{y}$ and $w_{x}$ shown in Fig. 6(a) and Fig. 7(a), in a
way similar to Ref. [13]. We denote this location of $\mathbf{P}$ as the base
position, whose location is indicated by the black star in Figs. 5(a) and (b).

Of course, there is uncertainty associated with both $\mathbf{P}$ and
$r_{\text{min}}$, especially the former, and it is important to quantify how
this affects the value $K_{\text{I}}$ determined from the measurements. We
assess the uncertainty in $\mathbf{P}$ by sampling multiple positions of
$\mathbf{P}$ in a circle of a radius of $r_{c}$ centered at the base position
of $\mathbf{P}$. Similarly, we apply a variation $r_{d}$ to the value of
$r_{\text{min}}$. We calculate the value of $K_{\text{I}}$ over each of these
ranges and report the average value. In practice, we have found that setting
$r_{c}=r_{d}=5$ pixels ($0.02\text{\,}\mathrm{m}\mathrm{m}$) works well.
Finally, we do not analyze data from the initial ten images starting from
crack initiation (i.e., crack lengths less than
$0.3\text{\,}\mathrm{m}\mathrm{m}$) due to the difficulty in identifying
$\mathbf{P}$ for such short cracks. We note that this could be alleviated by
using fatigue-precracked specimens, where the more pronounced stress
concentration would aid in the determination of $\mathbf{P}$.

##  5 Results

###  5.1 Measurement of surface slope $\bar{w}_{x}$ and $\bar{w}_{y}$

Fig. 6(a) shows the distribution of $\bar{w}_{y}$ across the field of view at
three instants during a representative dynamic fracture experiment from a
notched three-point bend sample of Macor™. We use the crack initiation time as
a reference ($t=0$). For this experiment with images collected at
$5\text{\,}\mathrm{M}\mathrm{f}\mathrm{p}\mathrm{s}$, the time between frames
is $0.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$. The spatial pattern
(shape and magnitude) of $\bar{w}_{y}$ is consistent with a calculation (Fig.
6(b)) using Eqn. 4.7 and the mechanical properties of Macor™
($E=$66.9\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$, $\nu=0.29$, and
$K_{\text{IC}}=$1.53\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}\,\mathrm{m}^{1/2}$$).
As an additional comparison we compute $\bar{w}_{y}$ (Fig. 6(c)) from dynamic
phase-field simulations [34, 35] of fracture using a custom code [36].
(Details of the phase-filed simulations are provided in the Appendix). This
computation is achieved by plugging in the in-plane stresses $\sigma_{xx}$ and
$\sigma_{yy}$ obtained from the simulation into Eqn. 4.7. We believe that the
highly localized surface slope profile in the simulation results in the wake
of the crack is an artifact of the phase-field algorithm, which smears the
sharp crack over a certain region, in which the stress does not completely
drop to zero.

Figure 6: (a) Spatial distribution of surface slope $\bar{w}_{y}$ determined
from SHWFS measurements at three different times after crack initiation
($t=0$). (b) Spatial distribution of $w_{y}$ calculated from the near-tip
stress field (analytical solution), given a crack tip at the center of the
image. (c) A schematic showing the dynamic phase-field simulation and three
plots showing the computed distribution of $\bar{w}_{y}$ at times comparable
to the experimental data in (a). The unit for all color bars is radians, and
that for all axes is millimeters.

Corresponding SHWFS measurements of $\bar{w}_{x}$ are shown in Fig. 7. Again
the measurements are qualitatively consistent with both analytical
calculations and phase-field simulations. However, the quantitative agreement
is not as good as that for $\bar{w}_{y}$ (Fig. 6). In particular, the SHWFS
measurement seems to underestimate the distribution of negative $\bar{w}_{x}$
behind the crack tip. We are still investigating the cause, but believe that
it may result from tilting of the specimen about the $y$ axis during loading.
Such a tilt will cause the front surface of the sample to rotate around the
$y$ axis, leading to a change in $\Delta x$ but leaving $\Delta y$ unchanged.
In the particular case shown in Fig. 7(b), the sample likely tilted along the
positive direction of $y$ (based on the right-hand rule), adding a positive
value to $\Delta x$. We note here, and show in the next section, that this
also affects the apparent stress intensity factor $K_{\text{I}}$ extracted
using $\bar{w}_{x}$.

Figure 7: Similar to Fig. 6 but for surface slopes $w_{x}$. (a) Spatial
distribution of $\bar{w}_{x}$ determined from SHWFS measurements. (b)
Analytical solution for spatial distribution of $w_{x}$. (c) A schematic
showing the dynamic phase-field simulation and three plots showing the
computed distribution of $\bar{w}_{x}$. The unit for all color bars is
radians, and that for all axes is millimeters.

###  5.2 Extraction of $K_{\text{I}}$

Fig. 8(a) shows how $K_{\text{I}}$ varies as a function of the apparent crack
length. When the crack length is less than $0.9\text{\,}\mathrm{m}\mathrm{m}$,
$K_{\text{I}}$ approaches from the below the critical stress intensity for
Macor™,
$K_{\text{IC}}=$1.53\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}\,\mathrm{m}^{1/2}$$.
This is possibly due to unloading [16] and the time required after crack
initiation to establish the singular stress field [37]. After the crack length
reaches $0.9\text{\,}\mathrm{m}\mathrm{m}$, $K_{\text{I}}$ is essentially
constant. We note that the value of
$K_{\text{IC}}=$1.53\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}\,\mathrm{m}^{1/2}$$
for Macor™ is that for a quasi-static crack growth, which is a reasonable
approximation for these conditions since the crack tip velocity is a small
fraction of the Rayleigh wave speed. The uncertainty in $K_{\text{I}}$ (shown
as error bars) in Fig. 8(a) is computed as the root mean square of the scatter
stemming from varying $\mathbf{P}$ and $r_{\text{min}}$, whose uncertainties
are quantified by $r_{c}$ and $r_{d}$ as discussed in Section. 4.2.

Although the $K_{\text{I}}$ trends determined using $\bar{w}_{y}$ and
$\bar{w}_{x}$ are similar for crack lengths smaller than about
$1\text{\,}\mathrm{m}\mathrm{m}$, they diverge beyond that point. In
particular, the value determined from $\bar{w}_{y}$ plateaus at the static
critical stress intensity as noted above, but the value determined from
$\bar{w}_{x}$ drops. We are still investigating the cause of this drop, but we
believe that it may result from the sample tilting around the $y$ axis during
loading (see Fig. 2(a)). This tilt may be due to a combination of imperfect
sample machining (along the sample thickness direction) and sample
disintegration due to crack propagation. It only adds changes to $\Delta x$
but not $\Delta_{y}$, and this can help explain the lack of negative
$\bar{w}_{x}$ values measured behind the crack tip, as shown in Fig. 8(c), and
in comparison with Fig. 8(b).

Figure 8: (a) Evolution of the apparent stress intensity factor $K_{\text{I}}$
as a function of the apparent crack length, using both the $\bar{w}_{y}$ data
(blue circles) and the $\bar{w}_{x}$ data (red circles). The error bar
associated with each data point indicates the scatter coming from varying
$\mathbf{P}$ and $r_{\text{min}}$. (b) A visualization of a particular
snapshot showing the spatial contour plot of $\bar{w}_{y}$, the slope data
selected for extracting $K_{\text{I}}$ (hollow circles), the base location of
the crack tip (black star), and the crack trajectory (grey stars). (c) A
similar visualization to (b) but showing the spatial contour plot of
$\bar{w}_{x}$.

##  6 Summary and outlook

In this work, we demonstrate a new approach for measuring stress intensity
during dynamic fracture of brittle solids using a Shack-Hartmann wavefront
sensor (SHWFS). With the help of a high-speed camera and high-intensity pulsed
laser illumination, we have successfully detected the surface slope
($\bar{w}_{x}$ and $\bar{w}_{y}$) evolution during the dynamic fracture of a
millimeter-scale Macor™ sample with a spatial resolution of
$45\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and a temporal resolution
of $0.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$. The value of the
stress intensity factor ($K_{\text{I}}$) we extract from SHWFS measurement of
$\bar{w}_{y}$ also agrees with the static fracture toughness of Macor™,
provided the crack propagates over a distance larger than about
$1\text{\,}\mathrm{m}\mathrm{m}$ after initiation. Compared to a typical DGS
apparatus, the SHWFS apparatus has improved spatial resolution while requiring
less physical space to implement. As such, SHWFS can be especially useful for
small-scale fracture experiments and is ideal for multi-modal studies,
particularly in combination with x-ray phase contrast imaging (XPCI) performed
at synchrotron facilities. One potential application in this regard is
studying the dynamic fracture of materials under high-rate loadings using
miniaturized desktop Kolsky bar setups [38] at synchrotron facilities [25]. In
addition, unlike DIC or DGS, SHWFS does not require the preparation of speckle
patterns (which is particularly challenging for small-scale experiments), and
it can be used on any material so long as the surface can be polished to a
mirror-like finish. The setup of SHWFS can also be modified to study
transparent materials without the need for coating a specimen’s sample to be
reflective, in a way similar to transmission-mode DGS [11]. As such, the
flexibility of SHWFS with regard to sample preparation also makes it adaptable
to studying a broader range of brittle materials. We discuss two directions
for improvement:

  * 1.

Increasing the resolution of the measurement. There are two aspects in this
direction: increasing the smallest detectable deformation (i.e., sensitivity)
and increasing the spatial resolution of detection. The smallest detectable
deformation is related to the shift of the reflected spots on the camera
sensor. Given a fixed field of view determined by the microscope objective,
the shift is tied to the focal length (denoted as $f$) of the multi-lens array
(MLA), with a larger focal length producing a larger shift (Eqn. 2.1 and 2.2).
The spatial resolution is tied to the pitch of the MLA (denoted as $a$,
essentially how large each microlens is). Since each microlens performs an
average measurement of surface slope across its area, we will want an MLA with
a small pitch (i.e., a smaller microlens) to increase the spatial resolution.
However, diffraction effects limit the smallest pitch (denoted as $a$) and the
largest focal length (denoted as $f$) of an MLA. Therefore, a tradeoff exists
between a small $a$ and a large $f$. Specifically, the Fresnel number
$\frac{a^{2}}{f\lambda}$ must be larger than unity where $\lambda$ is the
wavelength. Our current setup
($a=$150\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$, $\lambda=640$ nm,
and $L=5.6$ mm) leads to a Fresnel number of about 6.2. Choosing a different
MLA for our setup may lead to a better detection sensitivity with minimal
sacrifice to the spatial resolution or vice versa. Nevertheless, it is
important to note that the choice of MLA is application-specific: depending on
the specific experiment details (such as material type and size), a different
MLA may be needed for the desired detection sensitivity and spatial
resolution. Lastly, note that we also have the option of enlarging the field
of view (by using a microscope objective with a smaller magnification, see
Eqn. 4.3), but this can require a larger sample and lead to a change of the
loading apparatus.

  * 2.

Increasing the accuracy of the crack tip position determination. We have shown
that an accurate estimation of $K_{\text{I}}$ depends heavily on accurately
identifying the crack tip location (as also highlighted in [39, 40]). This is
quite challenging to do at the millimeter scale based solely on the
distribution pattern of the surface slope, especially given the presence of
noise in experimental measurements. In this regard, one possibility is
exploiting the distribution property of the out-of-plane asymptotic field, in
a way the extends the novel approach proposed by [40] which takes advantage of
the separability of the in-plane asymptotic field (e.g., the in-plane
displacement field measured by DIC). Another possibility is integrating SHWFS
and x-ray phase contrast imaging (XPCI), with the latter providing more
accurate time-resolved identification of the crack tip location. In XPCI,
phase perturbation introduced by the sample is exploited to modulate the
intensity recorded at the image detector plane [41]. Since the surface of a
crack induces a steep phase gradient, a phase-contrast image can have a
significantly enhanced contrast compared to conventional radiography. XPCI has
already been used to characterize the crack dynamics within different brittle
materials from ceramics [25] to glass [24, 29] taking advantage of spatially-
coherent high-energy X-rays provided by synchrotron facilities.

##  7 CRediT Authorship Contribution Statement

L Li. Conceptualization, Methodology, Experiment, Numerical Simulation, Formal
Analysis, Investigation, Visualization, Writing – Original & Draft. V. Kilic.
Methodology, Experiment, Formal Analysis, Investigation, Visualization,
Writing – Original & Draft. M. Alemohammad. Methodology, Experiment, Writing –
Review & Editing. KT Ramesh. Conceptualization, Resources, Supervision,
Funding Acquisition, Project Administration, Writing – Review & Editing. MA
Foster. Conceptualization, Methodology, Resources, Supervision, Project
Administration, Writing – Review & Editing. TC Hufnagel. Conceptualization,
Resources, Supervision, Funding acquisition, Project Administration, Writing –
Review & Editing.

##  8 Acknowledgements

We thank Dr. Jason Harris, Dr. Charlene Smith, and Dr. Xinyi Xu from Corning
Inc. for stimulating discussions on glass ceramics, as well as Nand Vommi for
his help with sample preparation, Justin Moreno and Matthew Shaeffer for their
assistance to our experiments, and Dr. Chengyun Miao for insightful
discussions on DGS. We gratefully acknowledge financial support provided by
the Corning Research and Development Corporation. The project was sponsored in
part by the Department of Defense, Defense Threat Reduction Agency under the
Materials Science in Extreme Environments University Research Alliance,
HDTRA1-20-2-0001. The content of the information does not necessarily reflect
the position or the policy of the federal government, and no official
endorsement should be inferred.

##  Appendix A Measuring the indentation speed

We recourse to image analysis to measure the indentation speed of the piezo
actuator (which drives the indenter), following the same procedure the lead
author has used to measure the angle of repose of granular materials [42]. We
record the motion of the bottom edge of the piezo actuator using a high-speed
camera, shown in Fig. 9(a), and we compute the position variation, shown in
Fig. 9(b) as a function of time by binarizing each image and finding the
boundary of the actuator in each image. From the position data, we can
calculate the instantaneous velocity, shown in Fig. 9(c), using a finite-
difference scheme. We see that the velocity is fairly constant. We calculate
the average velocity of the actuator following the above procedure, both with
and without indenting a sample. Fig. 9(d) shows the result. In short, the
indentation velocity stays largely unchanged, equaling around 0.13 m/s
regardless of whether indenting a sample or not.

Figure 9: (a)A schematic showing the steps of calculating the indentation
speed. (b) The position evolution of the actuator’s bottom edge as a function
of time without indenting a sample. (c) The corresponding speed evolution of
the actuator’s bottom edge as a function of time. (d) The time-averaged
indentation speed without and with indenting a sample (which is repeated three
times).

##  Appendix B Variational phase-field approach to fracture

The variational phase-field method to fracture [43] uses a scalar field
$\phi\in[0,1]$ to describe the spatial distribution of damage ($\phi=0$ being
intact and $\phi=1$ being completely damaged) of a simulation domain through a
length scale parameter $\ell$. Consequently, a crack (a sharp interface) is
smeared through a smooth variation of $\phi$ over $\ell$. The parameter $\ell$
represents active mechanisms in the process zone, determining the threshold
for crack nucleation [43]. Essentially, to model a dynamic fracture problem,
we want to minimize the following incremental Lagrange energy functional
$\mathbf{I}_{\ell}$ by the principle of least action [35]:

$\displaystyle\mathbf{I}_{\ell}(u,\dot{u},\phi)=\int_{t_{1}}^{t_{2}}\left\\{\int_{\Omega}\left[\frac{\rho}{2}|\dot{u}|^{2}-\mathcal{W}^{e}(u,\phi)-G_{\text{C}}\gamma_{\ell}(\phi,\nabla\phi)+\rho
b\cdot u\right]\text{d$\Omega$}+\int_{\partial\Omega}t\cdot
u\text{dS}\right\\}\text{dt},$ (B.1)

under the constraint $\dot{\phi}>0$ to account for the irreversibility of a
fracture process. Above, $u$ is the displacement field, with
$\dot{u}=\frac{\partial u}{\partial t}$ the velocity field, $\rho$ the
material density, $\phi$ the phase field parameter indicating the degree of
material damage, $\mathcal{W}^{e}$ the elastic strain energy density,
$G_{\text{C}}$ the critical energy release rate (or fracture toughness) [44],
$\gamma_{\ell}$ the (regularized) fracture energy density, $b$ the
gravitational constant, and $t$ the surface traction. We follow [34] for
modeling brittle solids, decomposing the elastic strain energy shown in Eqn.
B.1 to a tensile part ($``+$”) and a compressive part ($``-"$), with the
phase-field acting only on the former:

$\displaystyle\mathcal{W}^{e}(\epsilon_{ij},\phi)=[(1-k)(1-\phi)^{2}+k]\mathcal{W}^{e,+}(\epsilon_{ij})+\mathcal{W}^{e,-}(\epsilon_{ij}),$
(B.2)

where $\epsilon^{d}$ is the $d$-th eigenvalue of $\epsilon$, $n^{d}$ is the
corresponding eigenvector, $\langle x\rangle_{+}$ stands for $(x+|x|)/2$, and
$\langle x\rangle_{-}$ stands for $(x-|x|)/2$ with $|x|$ being the absolute
value of $x$. We can then express $\mathcal{W}^{e,+}(\epsilon_{ij})$ and
$\mathcal{W}^{e,-}(\epsilon_{ij})$ as the following:

$\displaystyle\mathcal{W}^{e,+}(\epsilon_{ij})$
$\displaystyle=\frac{1}{2}\lambda\langle\epsilon_{kk}\rangle_{+}^{2}+\mu\epsilon^{+}_{kj}\epsilon^{+}_{jl}\delta_{kl},$
(B.3) $\displaystyle\mathcal{W}^{e,-}(\epsilon_{ij})$
$\displaystyle=\frac{1}{2}\lambda\langle\epsilon_{kk}\rangle_{-}^{2}+\mu\epsilon^{-}_{kj}\epsilon^{-}_{jl}\delta_{kl},$
(B.4)

where $\lambda$ and $\mu$ are the Lamé constants that can be determined from
the Young’s modulus $E$ and the Poisson’s ratio $\nu$. The (regularized)
fracture energy density $\gamma_{\ell}$ takes the following form:

$\displaystyle\gamma_{\ell}=\frac{1}{4c_{w}\ell}\left(w(\phi)+\ell^{2}|\nabla\phi|^{2}\right),\text{with}\,c_{w}=\frac{1}{2}\,\,\text{which
implies}\,\,w(\phi)=\phi^{2}.$ (B.5)

Applying the principle of least action to Eqn. B.1 with $\mathcal{W}^{e}$ and
$\gamma_{\ell}$ expressed by Eqns. B.2, B.3, B.3, and B.5, we arrive at the
following two governing equations:

$\displaystyle\sigma_{ij,j}+b_{i}$ $\displaystyle=\rho\ddot{u}_{i},$ (B.6)
$\displaystyle\left[1+\frac{4c_{w}\ell(1-k)}{G_{\text{C}}}\mathcal{W}^{e,+}\right]\phi-\ell^{2}\phi_{,ii}$
$\displaystyle=\frac{4c_{w}\ell(1-k)}{G_{\text{C}}}\mathcal{W}^{e,+},$ (B.7)

where $\sigma_{ij}=\partial\mathcal{W}^{e}/\partial{\epsilon_{ij}}$. We
enforce the irreversible growth condition $\dot{\phi}>0$ using a strain-
history field [34] over the simulation domain:

$\displaystyle\mathcal{H}(x,t)=\max_{s\in[0,t]}\mathcal{W}^{e,+}\left(\epsilon(x,s)\right)\,\,\forall\,x\in\Omega.$
(B.8)

Replacing $\mathcal{W}^{e,+}$ with $\mathcal{H}(x,t)$ in Eqn. B.7 we then want
to solve:

$\displaystyle\sigma_{ij,j}+b_{i}$ $\displaystyle=\rho\ddot{u}_{i},$ (B.9)
$\displaystyle\left[1+\frac{4c_{w}\ell(1-k)}{G_{\text{C}}}\mathcal{H}\right]\phi-\ell^{2}\phi_{,ii}$
$\displaystyle=\frac{4c_{w}\ell(1-k)}{G_{\text{C}}}\mathcal{H},$ (B.10)

together with the following Neumann boundary conditions (plus any existing
Dirichlet boundary conditions) :

$\displaystyle\sigma_{ij}n_{j}$
$\displaystyle=t_{i}\,\,\text{on}\,\,\partial\Omega,$ (B.11)
$\displaystyle\phi_{,i}n_{i}$
$\displaystyle=0\,\,\text{on}\,\,\partial\Omega.$ (B.12)

We solve Eqns. B.9 and B.10 weakly follow a standard finite element
discretization and calculation procedure, using the alternating minimization
(or staggered) scheme that runs parallel on high-performance computer
clusters. We use the mechanical properties
$E=66.9\,\,\text{GPa},K_{\text{IC}}=$1.53\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}\,\mathrm{m}^{1/2}$$,
and $\nu=0.29$ in our PF simulations. The notched beam in the simulation
shares the same geometry as the experiment specimen, having a notch of
$250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ in width and 1 mm in
length. Since $\ell$ affects crack nucleation, it will affect values of
$w_{x}$ and $w_{y}$. Essentially, a smaller value of $\ell$ suggests more
delayed crack nucleation, leading to larger values of $w_{x}$ and $w_{y}$
because of a larger stress build-up before crack nucleation. We find
$\ell=0.035$ mm to be a reasonable choice, and the corresponding results are
presented in the main text. We pick an element size of $\delta h=0.015$ mm
near the crack propagation region, which is small enough ($\delta h\leq\ell/2$
[34]) to resolve the crack geometry. However, we emphasize that $\ell$ does
not change the spatial pattern of $w_{x}$ and $w_{y}$, which is a direct
consequence of the balance of linear momentum. Figs. 10(c) and 11(c) shows
results obtained from $\ell=0.15$ mm. Compared to Figs. 6(c) and 7(c), the
spatial pattern of $\bar{w}_{x}$ and $\bar{w}_{y}$ stays the same, but the
magnitude is smaller.

Figure 10: A similar figure to Fig. 6 but with plots shown in (c) obtained
from $\ell=0.15$ mm. Figure 11: A similar figure to Fig. 7 but with plots
shown in (c) obtained from $\ell=0.15$ mm.

## References

  * [1] James W Dally.  Dynamic photoelastic studies of fracture: Dynamic photoelasticity is applied to characterize the behavior of cracks propagating in polymers and a metal with an -k relation.  Experimental Mechanics, 19:349–361, 1979.
  * [2] K Ravi-Chandar.  An experimental investigation into the mechanics of dynamic fracture.  PhD thesis, California Institute of Technology, 1982.
  * [3] K Ravi Chandar and WG Knauss.  Dynamic crack-tip stresses under stress wave loading—a comparison of theory and experiment.  International Journal of Fracture, 20:209–222, 1982.
  * [4] AJ Rosakis, J Duffy, and LB Freund.  The determination of dynamic fracture toughness of aisi 4340 steel by the shadow spot method.  Journal of the Mechanics and Physics of Solids, 32(6):443–460, 1984\.
  * [5] Hareesh V Tippur, Sridhar Krishnaswamy, and Ares J Rosakis.  A coherent gradient sensor for crack tip deformation measurements: analysis and experimental results.  International journal of fracture, 48:193–204, 1991.
  * [6] Hareesh V Tippur, Sridhar Krishnaswamy, and Ares J Rosakis.  Optical mapping of crack tip deformations using the methods of transmission and reflection coherent gradient sensing: a study of crack tip k-dominance.  International Journal of Fracture, 52:91–117, 1991.
  * [7] Sridhar Krishnaswamy, Hareesh V Tippur, and Ares J Rosakis.  Measurement of transient crack-tip deformation fields using the method of coherent gradient sensing.  Journal of the Mechanics and Physics of Solids, 40(2):339–372, 1992\.
  * [8] RJ Butcher, C-E Rousseau, and HV Tippur.  A functionally graded particulate composite: preparation, measurements and failure analysis.  Acta Materialia, 47(1):259–268, 1998.
  * [9] Madhu S Kirugulige, Hareesh V Tippur, and Thomas S Denney.  Measurement of transient deformations using digital image correlation method and high-speed photography: application to dynamic fracture.  Applied optics, 46(22):5083–5096, 2007.
  * [10] Gabriele Albertini, Mathias Lebihain, François Hild, Laurent Ponson, and David S Kammer.  Effective toughness of heterogeneous materials with rate-dependent fracture energy.  Physical Review Letters, 127(3):035501, 2021.
  * [11] Chandru Periasamy and Hareesh V Tippur.  Full-field digital gradient sensing method for evaluating stress gradients in transparent solids.  Applied Optics, 51(12):2088–2097, 2012.
  * [12] C Miao and HV Tippur.  Higher sensitivity digital gradient sensing configurations for quantitative visualization of stress gradients in transparent solids.  Optics and Lasers in Engineering, 108:54–67, 2018.
  * [13] Balamurugan M Sundaram and Hareesh V Tippur.  Dynamic fracture of soda-lime glass: A full-field optical investigation of crack initiation, propagation and branching.  Journal of the Mechanics and Physics of Solids, 120:132–153, 2018\.
  * [14] S Dondeti and HV Tippur.  A comparative study of dynamic fracture of soda-lime glass using photoelasticity, digital image correlation and digital gradient sensing techniques.  Experimental Mechanics, 60:217–233, 2020.
  * [15] Balamurugan M Sundaram and Hareesh V Tippur.  Dynamic mixed-mode fracture behaviors of pmma and polycarbonate.  Engineering Fracture Mechanics, 176:186–212, 2017.
  * [16] S Dondeti and HV Tippur.  Cascading crack bifurcations in soda-lime glass: Quantification of fracture mechanics-based precursors using digital gradient sensing.  International Journal of Solids and Structures, 234:111252, 2022\.
  * [17] Wenfeng Hao, Can Tang, Yanan Yuan, Xuefeng Yao, and Yinji Ma.  Experimental study on the fiber pull-out of composites using digital gradient sensing technique.  Polymer Testing, 41:239–244, 2015.
  * [18] Chengyun Miao and Hareesh V Tippur.  Fracture behavior of carbon fiber reinforced polymer composites: an optical study of loading rate effects.  Engineering Fracture Mechanics, 207:203–221, 2019.
  * [19] Chandru Periasamy and Hareesh V Tippur.  A full-field reflection-mode digital gradient sensing method for measuring orthogonal slopes and curvatures of thin structures.  Measurement Science and Technology, 24(2):025202, 2013.
  * [20] Chengyun Miao and Hareesh V Tippur.  A simplified reflection-mode digital gradient sensing technique for measuring surface slopes, curvatures and topography.  Optics and Lasers in Engineering, 124:105843, 2020.
  * [21] Roland V Shack.  Production and use of a lenticular hartmann screen.  In Spring Meeting of Optical Society of America, 1971, volume 656, 1971.
  * [22] Ben C Platt and Roland Shack.  History and principles of shack-hartmann wavefront sensing, 2001.
  * [23] SN Luo, BJ Jensen, DE Hooks, K Fezzaa, KJ Ramos, JD Yeager, K Kwiatkowski, and T Shimada.  Gas gun shock experiments with single-pulse x-ray phase contrast imaging and diffraction at the advanced photon source.  Review of Scientific Instruments, 83(7), 2012.
  * [24] Niranjan D Parab, John T Black, Benjamin Claus, Matthew Hudspeth, Jianzhuo Sun, Kamel Fezzaa, and Weinong W Chen.  Observation of crack propagation in glass using x-ray phase contrast imaging.  International Journal of Applied Glass Science, 5(4):363–373, 2014\.
  * [25] Andrew FT Leong, Andrew K Robinson, K Fezzaa, T Sun, N Sinclair, DT Casem, PK Lambert, CJ Hustedt, Nitin P Daphalapurkar, KT Ramesh, et al.  Quantitative in situ studies of dynamic fracture in brittle solids using dynamic x-ray phase contrast imaging.  Experimental Mechanics, 58:1423–1437, 2018.
  * [26] Lawrence E Schmutz.  Hartmann sensing at adaptive optics associates.  In Electromechanical system interaction with optical design, volume 779, pages 13–17. SPIE, 1987.
  * [27] Bing Pan, Huimin Xie, Zhaoyang Wang, Kemao Qian, and Zhiyong Wang.  Study on subset size selection in digital image correlation for speckle patterns.  Optics express, 16(10):7037–7048, 2008.
  * [28] Kim I Mortensen, L Stirling Churchman, James A Spudich, and Henrik Flyvbjerg.  Optimized localization analysis for single-molecule tracking and super-resolution microscopy.  Nature methods, 7(5):377–381, 2010.
  * [29] M Kang, W Li, AFT Leong, M Guan, K Fezzaa, JT Harris, KT Ramesh, and TC Hufnagel.  Crack nucleation and growth during dynamic indentation of chemically-strengthened glass.  Extreme Mechanics Letters, 38:100754, 2020.
  * [30] Sung-Hoon Baik, Seung-Kyu Park, Cheol-Jung Kim, and Byungheon Cha.  A center detection algorithm for shack–hartmann wavefront sensor.  Optics & Laser Technology, 39(2):262–267, 2007.
  * [31] Max L Williams.  On the stress distribution at the base of a stationary crack.  1957\.
  * [32] Lambert Ben Freund.  Dynamic fracture mechanics.  Cambridge university press, 1998.
  * [33] Ilya Svetlizky and Jay Fineberg.  Classical shear cracks drive the onset of dry frictional motion.  Nature, 509(7499):205–208, 2014.
  * [34] Christian Miehe, Martina Hofacker, and Fabian Welschinger.  A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits.  Computer Methods in Applied Mechanics and Engineering, 199(45-48):2765–2778, 2010.
  * [35] Michael J Borden, Clemens V Verhoosel, Michael A Scott, Thomas JR Hughes, and Chad M Landis.  A phase-field description of dynamic brittle fracture.  Computer Methods in Applied Mechanics and Engineering, 217:77–95, 2012.
  * [36] Liuchi Li, Jack Rao, Todd Hufnagel, and KT Ramesh.  Meso-scale size effects of material heterogeneities on crack propagation in brittle solids.  arXiv preprint arXiv:2309.12916, 2023.
  * [37] Krishnaswa Ravi-Chandar and Wolfgang G Knauss.  An experimental investigation into dynamic fracture: Iii. on steady-state crack propagation and crack branching.  International Journal of fracture, 26:141–154, 1984.
  * [38] Sh Akhondzadeh, Minju Kang, Ryan B Sills, KT Ramesh, and Wei Cai.  Direct comparison between experiments and dislocation dynamics simulations of high rate deformation of single crystal copper.  Acta Materialia, 250:118851, 2023.
  * [39] A Taylor Owens and HV Tippur.  An image processing technique to identify crack tip position and automate fracture parameter extraction using dic: Application to dynamic fracture.  Experimental Mechanics, 63(3):445–466, 2023.
  * [40] Swati Gupta, Grant West, Mark A Wilson, Scott J Grutzik, and Derek H Warner.  Identifying crack tip position and stress intensity factors from displacement data.  International Journal of Fracture, pages 1–17, 2023.
  * [41] Marco Endrizzi.  X-ray phase-contrast imaging.  Nuclear instruments and methods in physics research section A: Accelerators, spectrometers, detectors and associated equipment, 878:88–98, 2018\.
  * [42] Liuchi Li and José E Andrade.  Identifying spatial transitions in heterogenous granular flow.  Granular Matter, 22(2):52, 2020.
  * [43] Blaise Bourdin, Gilles A Francfort, and Jean-Jacques Marigo.  The variational approach to fracture.  Journal of elasticity, 91:5–148, 2008.
  * [44] Alan Arnold Griffith.  Vi. the phenomena of rupture and flow in solids.  Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character, 221(582-593):163–198, 1921.