# A new closed-form expression for the solution of ODEs in a ring of
distributions and its connection with the matrix algebra

S. Pozzaa CONTACT S. Pozza. Email<EMAIL_ADDRESS>aCharles
University, Sokolovská 83 186, 75 Prague 8, Czech Republic

###### Abstract

A new expression for solving homogeneous linear ODEs based on a generalization
of the Volterra composition was recently introduced. In this work, we extend
such an expression, showing that it corresponds to inverting an infinite
matrix. This is done by studying a particular subring and connecting it with a
subalgebra of infinite matrices.

###### keywords:

Ordinary differential equations; Volterra composition; Legendre polynomials

##  1 Introduction

Let $\tilde{A}(t)$ be an $N\times N$ matrix-valued function analytic over
$t\in\mathcal{I}=[0,1]$ and $I_{N}$ the $N\times N$ identity matrix. Then, the
system of ODEs

$\frac{d}{dt}U_{s}(t)=\tilde{A}(t)U_{s}(t),\quad U_{s}(s)=I_{N},\quad\text{
for }t\geq s,\quad t,s\in\mathcal{I},$ (1)

has a unique solution $U_{s}(t)$. When
$\tilde{A}(\tau_{1})\tilde{A}(\tau_{2})=\tilde{A}(\tau_{2})\tilde{A}(\tau_{1})$
for every $\tau_{1},\tau_{2}\in\mathcal{I}$, $U_{s}(t)$ takes the form

$U_{s}(t)=\exp\left(\int_{s}^{t}\tilde{A}(\tau)\,\text{d}\tau\right).$

In general, however, $U_{s}(t)$ has no known simple expression in terms of
$\tilde{A}(t)$. Indeed, despite systems of non-autonomous linear ODEs are
crucial, common problems that appear in a variety of contexts [1, 2, 3, 4, 5,
6], their solution is surprisingly difficult to formulate by an analytic
expression.

When $\tilde{A}(t)$ is a scalar function, [7] shows that a closed form of the
solution exists in the non-commutative ring $\mathcal{S}$ composed of a
certain distribution set $\mathcal{D}(\mathcal{I})$ [8], the so-called
$\star$-product [9], and the usual addition. The $\star$-product is a
convolution-like operation that generalizes the _Volterra composition_ (e.g.,
[10]). The closed form is given in terms of a $\star$-product inverse in the
ring. Moreover, it is easy to define a $\mathcal{S}$-module of matrices with a
bilinear product that generalizes the results to the case of a matrix-valued
$\tilde{A}(t)$; see [11, 12]. In a few words, in this case, the solution
$U_{s}(t)$ is given by the bilinear product inverse of a matrix in the
$\mathcal{S}$-module. This means that in the framework of the $\mathcal{S}$
ring, it is possible to express $U_{s}(t)$ in a closed form for every matrix-
valued analytic function $\tilde{A}(t,s)$. This new expression has led to
several new symbolic and numerical approaches to the solution of (1) [11, 13,
9, 12, 14, 15, 16, 7].

In the pieces of literature mentioned above, the new expression for the
solution of (1) has not been derived in the ring $\mathcal{S}$, but in
alternative equivalent ways. This paper aims to show the potentiality of
working in the $\mathcal{S}$-ring module. We do that by deriving a new result,
namely, the expression for the solution of the non-homogeneous system of
linear ODEs

$\tilde{A}(t)U_{s}(t)=\frac{d}{dt}U_{s}(t)+\tilde{B}(t),\quad
U_{s}(s)=I_{N},\quad\text{ for }t\geq s,\quad t,s\in\mathcal{I},$ (2)

where $\tilde{B}(t)$ is an $N\times N$ matrix-valued analytic function over
$\mathcal{I}$. Moreover, we will show that there is a subring of $\mathcal{S}$
that corresponds to a subalgebra of infinite matrices, and we will prove the
existence of certain matrix inverses in the subalgebra using the connection
with $\mathcal{S}$.

In Section 2, we define the $\star$-product and the related algebraic
structures, and we derive the new expression of the solution of (2). Section 3
shows the connection between the $\mathcal{S}$ subring and a subalgebra of
infinite matrices. As a consequence, the ODE solution can be obtained by
solving a linear system in the subalgebra. Section 4 concludes the
presentation.

##  2 A $\star$-product solution to non-homogeneous ODEs

Let $\tilde{f}_{1}(t,s),\tilde{f}_{2}(t,s)$ be two bivariate functions and
assume that they are analytic111Note that in the previously appeared works, we
have usually assumed the functions to be smooth. Here we restrict the
assumption to analytic for the sake of simplicity., in both $t$ and $s$, over
$\mathcal{I}=[0,1]$; we denote such a set of functions by
$\mathcal{A}(\mathcal{I})$. The _Volterra composition_ of
$\tilde{f}_{1},\tilde{f}_{2}$, introduced by Vito Volterra (e.g., [10]), is
defined as

$\big{(}\tilde{f}_{2}\star_{v}\tilde{f}_{1}\big{)}(t,s):=\int_{s}^{t}\tilde{f}_{2}(t,\tau)\tilde{f}_{1}(\tau,s)\,\text{d}\tau,\quad
t,s\in\mathcal{I}.$

Note that, from now on, a function marked with a tilde will stand for a
function from $\mathcal{A}(\mathcal{I})$. If we look at it as a product, the
Volterra composition lacks important features. For instance, the identity.
This is why the Volterra composition has been extended to the so-called
$\star$-product [9]. Let $\Theta(t-s)$ be the Heaviside theta function, i.e.,

$\Theta(t-s)=\begin{cases}1,\quad t\geq s\\\ 0,\quad t<s\end{cases}.$

Moreover, let $\delta(\cdot)=\delta^{(0)}(\cdot)$ be the Dirac delta
distribution and $\delta^{(i)}(\cdot)$ be its $i$th derivatives. We denote
with $\mathcal{D}(\mathcal{I})$ the class of the distributions $d$ that can be
expressed as

$d(t,s)=\widetilde{d}(t,s)\Theta(t-s)+\sum_{i=0}^{k}\widetilde{d}_{i}(t,s)\delta^{(i)}(t-s),$

with $\tilde{d},\tilde{d}_{i}\in\mathcal{A}(\mathcal{I})$. The $\star$-product
$\star:\mathcal{D}(\mathcal{I})\times\mathcal{D}(\mathcal{I})\rightarrow\mathcal{D}(\mathcal{I})$
is defined as

$\big{(}f_{2}\star
f_{1}\big{)}(t,s):=\int_{\mathcal{I}}f_{2}(t,\tau)f_{1}(\tau,s)\,\text{d}\tau,\quad
f_{1},f_{2}\in\mathcal{D}(\mathcal{I}),$ (3)

Consider the subclass
$\mathcal{A}_{\Theta}(\mathcal{I})\subset\mathcal{D}(\mathcal{I})$ comprising
those distributions of the form

$f(t,s)=\widetilde{f}(t,s)\Theta(t-s).$

Then, the $\star$-product of $f_{1},f_{2}\in\mathcal{A}_{\Theta}(\mathcal{I})$
is equivalent to the Volterra composition

$\displaystyle\big{(}f_{2}\star f_{1}\big{)}(t,s)$
$\displaystyle=\int_{\mathcal{I}}\widetilde{f}_{2}(t,\tau)\widetilde{f}_{1}(\tau,s)\Theta(t-\tau)\Theta(\tau-s)\,\text{d}\tau,$
$\displaystyle=\Theta(t-s)\int_{s}^{t}\widetilde{f}_{2}(t,\tau)\widetilde{f}_{1}(\tau,s)\,\text{d}\tau=\Theta(t-s)(\tilde{f}_{2}\star_{v}\tilde{f}_{1})(t,s).$

The $\star$-product is well-defined and closed in $\mathcal{D}(\mathcal{I})$;
we refer the reader to [12, 9] for further details. For the goals of this
paper, it will be enough to recall the following properties. Given
$f\in\mathcal{A}_{\Theta}(\mathcal{I})$, then

$\displaystyle\left(\delta^{\prime}(t-s)\ast f\right)(t,s)$
$\displaystyle=\left(\partial_{t}\tilde{f}(t,s)\right)\Theta(t-s)+\tilde{f}(s,s)\delta(t-s);$
(4) $\displaystyle\left(f\ast\delta^{\prime}(t-s)\right)(t,s)$
$\displaystyle=-\left(\partial_{s}\tilde{f}(t,s)\right)\Theta(t-s)+\tilde{f}(t,t)\delta(t-s);$

see [8, 9]. As a consequence,

$\Theta(t-s)\star\delta^{\prime}(t-s)=\delta^{\prime}(t-s)\star\Theta(t-s)=\delta(t-s),$
(5)

i.e., $\delta^{\prime}$ is the $\star$-inverse of $\Theta$. Moreover,

  * •

$\mathcal{D}(\mathcal{I})$ is closed under $\star$-multiplication;

  * •

the $\star$-product is associative over $\mathcal{D}(\mathcal{I})$;

  * •

the Dirac delta distribution $1_{\star}(t,s):=\delta(t-s)$ is the identity of
the $\star$-product.

Therefore,
$\mathcal{S}(\mathcal{I}):=(\mathcal{D}(\mathcal{I}),\star,+,0,1_{\star})$ is
a non-commutative ring. The $\star$-product can also be extended to matrices
and vectors composed of elements from $\mathcal{D}(\mathcal{I})$. This is
easily done by replacing the standard multiplication appearing in the
integrand of (3) with the usual matrix-matrix multiplication [14]. Similarly,
we can define the right (and left) scalar-matrix multiplication. As a result,
we obtain the module of the matrices with elements from
$\mathcal{D}(\mathcal{I})$ with, as bilinear product, the $\star$-product
between matrices, and, as scalar product, the $\star$-product between a scalar
and a matrix.

The system of ODEs in (2) can be rewritten in the form

$\partial_{t}U(t,s)=\tilde{A}(t)\Theta(t-s)U(t,s)+\tilde{B}(t,s)\Theta(t-s),\quad
U(s,s)=I_{N},\quad t,s\in\mathcal{I},$ (6)

with $U_{s}(t)=U(t,s)$. Note that the matrices
$U(t,s)=\tilde{U}(t,s)\Theta(t-s)$, $A(t,s):=\tilde{A}(t)\Theta(t-s)$, and
$B(t,s):=\tilde{B}(t,s)\Theta(t-s)$ are all composed of elements from
$\mathcal{A}_{\Theta}(\mathcal{I})$. Therefore, by exploiting formula (4) and
(5), equation (6) becomes:

$\delta^{\prime}(t-s)\star U(t,s)=\tilde{A}(t)U(t,s)+I_{\star}(t,s)+B(t,s),$
(7)

where $I_{\star}(t,s)$ is the identity matrix $I_{N}$ multiplied by
$\delta(t-s)=1_{\star}(t,s)$. Once the problem has been rewritten into the
$\star$-framework, we can derive a formula for its solution by working in the
$\mathcal{S}(\mathcal{I})$-module. If we $\star$-multiplying (7) from the left
by $\Theta(t-s)$ we obtain

$U(t,s)=\Theta(t-s)\star\left(\tilde{A}(t)U(t,s)+I_{\star}(t,s)+B(t,s)\right).$
(8)

Now, by replacing $U(t,s)$ in the right-hand side of (8) with the right-hand
side of (8) itself, we get the following iterations (we drop the dependency
from $t,s$ for the sake of readability)

$\displaystyle U$
$\displaystyle=\Theta\star\left(\tilde{A}U+I_{\star}+B\right),$
$\displaystyle=\Theta\star\left(\tilde{A}\left(\Theta\star\left(\tilde{A}U+I_{\star}+B\right)\right)+I_{\star}+B\right),$
$\displaystyle=\Theta\star\left(A\star\left(\tilde{A}U+I_{\star}+B\right)+I_{\star}+B\right),$
$\displaystyle=\Theta\star\left(A\star\tilde{A}U+\left(A+I_{\star}\right)\star(I_{\star}+B)\right).$

Note that the equality

$\tilde{A}\left(\Theta\star\left(\tilde{A}U+I_{\star}+B\right)\right)=A\star\left(\tilde{A}U+I_{\star}+B\right)$

holds since $\tilde{A}(t)$ does not depend on $s$. Repeating the iterations
$k$ times222In fact, such iterations are Picard iterations, see [16, Section
2]., we obtain

$\displaystyle U$ $\displaystyle=\Theta\star\left(A\star
A\star\tilde{A}U+\left(A\star A+A+I_{\star}\right)\star(I_{\star}+B)\right),$
$\displaystyle\,\,\,\vdots$
$\displaystyle=\Theta\star\left(A^{k\star}\star\tilde{A}U+\left(A^{k\star}+\dots+A+I_{\star}\right)\star(I_{\star}+B)\right),$
(9)

with $A^{k\star}$ the $k$th $\star$-power of $A$. As shown in [14],

$\max_{t,s\in\mathcal{I}}\left\|\left(A^{k\star}\right)(t,s)\right\|\leq\left(\max_{t,s\in\mathcal{I}}\|A(t,s)\|\right)^{k}\frac{(t-s)^{k-1}}{(k-1)!},\quad
k\geq 1,$

for any induced matrix norm. Therefore (9) uniformly converges to the
expression

$U(t,s)=\Theta(t-s)\star
R_{\star}(A)(t,s)\star\left(I_{\star}(t,s)+B(t,s)\right),$ (10)

where $R_{\star}(A)$ is the $\star$-resolvent of $A$, i.e.,

$R_{\star}(A)=I_{\star}+\sum_{k=1}^{\infty}\left(A^{\star k}\right)(t,s).$

Noticing that

$R_{\star}(A)\star(I_{\star}-A)=\left(I_{\star}+\sum_{k\geq 1}A^{\star
k}\right)\star(I_{\star}-A)=I_{\star},$ (11)

that is, $R_{\star}(A)=(I_{\star}-A)^{-\star}$ (the $\star$-inverse of
$(I_{\star}-A)$), we get

$U(t,s)=\Theta(t-s)\star(I_{\star}-A)^{-\star}(t,s)\star\left(I_{\star}(t,s)+B(t,s)\right),$
(12)

that is a closed-form expression in the $\mathcal{S}$-module.

Finally, since the matrix $A$ is composed of elements from the subset

$\mathcal{A}_{\Theta}^{t}(\mathcal{I}):=\left\\{f\in\mathcal{A}_{\Theta}(\mathcal{I}):f(t,s)=\tilde{f}(t)\Theta(t-s)\right\\}\subset\mathcal{A}_{\Theta}(\mathcal{I}),$

it is useful to define the set

$\mathcal{D}^{t}_{0}(\mathcal{I}):=\left\\{f(t,s)=\alpha
1_{\star}+\sum_{i=1}^{n}\big{(}g_{i,1}\star\dots\star
g_{i,m_{n}}\big{)},\;g_{i,j}\in\mathcal{A}_{\Theta}^{t}(\mathcal{I}),\;\alpha\in\mathbb{C}\right\\}.$

and the related subring
$(\mathcal{D}^{t}_{0}(\mathcal{I}),\star,+,0,1_{\star})$.

##  3 The $\star$-product and the matrix algebra

Let $\\{p_{k}\\}_{k}$ be a sequence of orthonormal shifted Legendre
polynomials over the bounded interval $\mathcal{I}=[0,1]$, i.e.,

$\displaystyle\int_{\mathcal{I}}p_{k}(\tau)p_{\ell}(\tau)d\tau=\delta_{k,\ell}=\begin{cases}0,\quad\text{if
}k\neq\ell\\\ 1,\quad\text{if }k=\ell\end{cases}.$

Despite the fact that the functions $p_{k}$ are not in
$\mathcal{D}(\mathcal{I})$, with a small abuse of notation, we can still
define the product

$p_{k}(s)\star
p_{\ell}(t)=\int_{\mathcal{I}}p_{k}(\tau)p_{\ell}(\tau)\;\textrm{d}\tau=\delta_{k,\ell}.$

Given a function
$f(t,s)=\tilde{f}(t,s)\Theta(t-s)\in\mathcal{A}_{\Theta}(\mathcal{I})$, we can
expand it into the series

$f(t,s)=\sum_{k=0}^{\infty}\sum_{\ell=0}^{\infty}f_{k,\ell}\,p_{k}(t)p_{\ell}(s),\quad
t\neq s,\quad t,s\in\mathcal{I},$ (13)

with coefficients

$f_{k,\ell}=\int_{\mathcal{I}}\int_{\mathcal{I}}f(\tau,\rho)p_{k}(\tau)p_{\ell}(\rho)\;\textrm{d}\rho\;\textrm{d}\tau;$

see, e.g., [17, p. 55]. The expansion can be rewritten in the matrix form

$\displaystyle
f(t,s)=\sum_{k=0}^{\infty}\sum_{\ell=0}^{\infty}f_{k,\ell}\,p_{k}(t)p_{\ell}(s)=\phi(t)^{T}F\,\phi(s)\quad
t\neq s,\quad t,s\in\mathcal{I}.$

where the _coefficient matrix_ $F$ and the vector $\phi(\tau)$ are defined as
follows

$F:=\begin{bmatrix}f_{0,0}&f_{0,1}&f_{0,2}&\dots\\\
f_{1,0}&f_{1,1}&f_{1,2}&\dots\\\ f_{2,0}&f_{2,1}&f_{2,2}&\dots\\\
\vdots&\vdots&\vdots&\ddots\end{bmatrix},\quad\phi_{M}(\tau):=\begin{bmatrix}p_{0}(\tau)\\\
p_{1}(\tau)\\\ p_{2}(\tau)\\\ \vdots\end{bmatrix}.$ (14)

Note that each element of $F$ can be bounded by

$|f_{k,\ell}|\leq\max_{t,s\in[0,1]}|\tilde{f}(t,s)|\sqrt{2k+1}\sqrt{2\ell+1}$
(15)

In particular, an univariate function $\tilde{f}(t)$ can be expanded as

$\tilde{f}(t):=\sum_{d=0}^{\infty}\alpha_{d}p_{d}(t),\quad\text{with
}\alpha_{d}=\int_{-1}^{1}\tilde{f}(t)p_{d}(t)dt.$

Let $B^{(k)}$ be the coefficient matrix of $p_{k}(t)\Theta(t-s)$. Then, the
coefficient matrix of $f(t,s)=\tilde{f}(t)\Theta(t-s)$ can be also expanded
into the series

$F=\sum_{k=0}^{\infty}\alpha_{k}B^{(k)}.$

Note that each $B^{(k)}$ is a banded matrix with bandwidth $k+1$, [7].
Moreover, the Fourier coefficients $\\{\alpha_{k}\\}_{k\geq 0}$ decay
geometrically [18]. Indeed, there exist $0<\rho<1,C>0$ such that

$|a_{k}|\leq C\rho^{k}.$ (16)

As a consequence, each element of $F$ can be bounded as follows

$\displaystyle|f_{k,\ell}|$
$\displaystyle=\left|\sum_{j=|k-\ell|+2}^{\infty}\alpha_{j}B^{(j)}_{k,\ell}\right|\leq\sum_{j=|k-\ell|+2}^{\infty}|\alpha_{j}||B^{(j)}_{k,\ell}|$
$\displaystyle\leq
C\max_{t,s\in[0,1]}|\tilde{f}(t,s)|\sum_{j=|k-\ell|+2}^{\infty}\rho^{j}\sqrt{2k+1}\sqrt{2\ell+1}$
$\displaystyle\leq K\rho^{|k-\ell|+2},$ (17)

for some $K>0$. This means that $F$ is characterized by a geometric decay of
the element magnitude as we move away from the diagonal.

Consider
$f(t,s)=\tilde{f}(t)\Theta(t-s),g(t,s)=\tilde{g}(t)\Theta(t-s)\in\mathcal{A}_{\Theta}^{t}(\mathcal{I})$
and the related coefficient matrices (14) $F,G$, respectively. By
orthogonality, $\phi(s)\star\phi(t)^{T}=I$, with $I$ the identity matrix.
Therefore, for every $t\neq s$,

$\displaystyle(f\star g)(t,s)$
$\displaystyle=\left(\phi(t)^{T}F\,\phi(s)\right)\star\left(\phi(t)^{T}G\,\phi(s)\right),$
$\displaystyle=\phi(t)^{T}F\left(\phi(s)\star\phi(t)^{T}\right)G\,\phi(s),$
$\displaystyle=\phi(t)^{T}FG\,\phi(s).$

Thus, the coefficient matrix $H$ of the function $h=f\star
g\in\mathcal{D}_{0}^{t}(\mathcal{I})$ is given by the matrix-matrix product
$FG$. It is important to note that the product $FG$ is well-defined. Indeed,
the series is convergent since by (17) there exist $K>0$ and $0<\rho<1$ so
that

$\displaystyle\left|(FG)_{k,\ell}\right|$
$\displaystyle=\left|\sum_{j=1}^{\infty}F_{k,j}G_{j,\ell}\right|\leq\sum_{j=1}^{\infty}\left|F_{k,j}\right|\left|G_{j,\ell}\right|$
$\displaystyle\leq\sum_{j=1}^{\infty}K_{f}K_{g}\rho_{f}^{|k-j|+2}\rho_{g}^{|\ell-j|+2}\leq\sum_{j=1}^{\infty}K\rho^{|k-j|+|\ell-j|+4}.$

Moreover, since $\min_{j=1,2\dots}|k-j|+|\ell-j|=|k-\ell|$, there exist
$K_{fg}>0$ and $0<\rho_{fg}<1$ so that

$\displaystyle\left|(FG)_{k,\ell}\right|\leq K_{fg}\rho_{fg}^{|k-\ell|}.$ (18)

This latter bound show that $FG$ is also characterized by a geometric decay of
the element magnitude as we move away from the diagonal. Therefore, given
$F,G,H$ coefficient matrices of
$f(t,s)=\tilde{f}(t)\Theta(t-s),g(t,s)=\tilde{g}(t)\Theta(t-s),h(t,s)=\tilde{h}(t)\Theta(t-s)$,
the matrix product $FGH$ is also well-defined and characterized by a off-
diagonal geometric decay. As a consequence, the set $\mathcal{F}$ of all the
coefficient matrices of functions from $\mathcal{D}_{0}^{t}(\mathcal{I})$ is a
subalgebra (with the usual sum, product, and matrix product) and it
corresponds to the subring
$(\mathcal{D}_{0}^{t}(\mathcal{I}),\star,+,0,1_{\star})$.

Consider now the $N\times N$ matrix-valued functions
$A(t,s)=[f_{i,j}(t,s)]_{i,j=1}^{N},B(t,s)=[g_{i,j}(t,s)]_{i,j=1}^{N}\in\mathbb{C}$
composed of elements from $\mathcal{D}_{0}^{t}(\mathcal{I})$. The functions
$f_{i,j}$ and $g_{i,j}$ are associated with their coefficient matrices
$F^{(i,j)},G^{(i,j)}$, respectively. By extending the arguments presented
above, we get the following expression for the (matrix) $\star$-product
$C=A\star B=[h_{i,j}]_{i,j=1}^{N}$,

$\left(A\star B\right)_{k,\ell}(t,s)=\sum_{j=1}^{N}(f_{k,j}\star
g_{j,\ell})(t,s)=\sum_{j=1}^{N}\phi(t)^{T}F^{(k,j)}G^{(j,\ell)}\,\phi(s),\quad
t\neq s.$

Therefore, the coefficient matrices $H^{(k,\ell)}$ of the functions
$h_{k,\ell}(t,s)$ are given by

$H^{(k,\ell)}=\sum_{j=1}^{N}F^{(k,j)}G^{(j,\ell)}.$ (19)

Defining the block matrices $\mathbf{A}=[F_{i,j}]_{i,j=1}^{N}$,
$\mathbf{B}=[G_{i,j}]_{i,j=1}^{N}$, $\mathbf{C}=[H_{i,j}]_{i,j=1}^{N}$, we
obtain the relation:

$\mathbf{C}=\mathbf{A}\mathbf{B}.$ (20)

Note that, despite the blocks having an infinite size, the product in (20) is
well-defined since the matrix products in (19) are well-defined.

In Section 2, we have seen the crucial role played by the $\star$-resolvent
$R_{\star}(A)$. Let $C:=\sum_{k\geq 1}A^{\star k}=[h_{i,j}]_{i,j=1}^{N}$, with
$h_{i,j}\in\mathcal{A}_{\Theta}$. Therefore, using the notation above,

$\displaystyle\left(R_{\star}(A)\right)_{k,\ell}(t,s)$
$\displaystyle=\phi(t)^{T}\left(I+H^{(k,\ell)}\right)\phi(s).$

Hence, for $t\neq s$, we obtain

$\displaystyle\left(R_{\star}(A)\star(I_{\star}-A)\right)_{k,\ell}$
$\displaystyle=\sum_{j=1}^{N}\phi(t)^{T}\left(I+H^{(k,j)}\right)\phi(s)\star\phi(t)^{T}\left(I-F^{(j,\ell)}\right)\phi(s),$
(21)
$\displaystyle=\phi(t)^{T}\sum_{j=1}^{N}\left(I+H^{(k,j)}\right)\left(I-F^{(j,\ell)}\right)\phi(s).$
(22)

Then, relation (11) implies

$\displaystyle\phi(t)^{T}\sum_{j=1}^{N}\left(I+H^{(k,j)}\right)\left(I-F^{(j,\ell)}\right)\phi(s)=1_{\star}\delta_{k,\ell}=\phi(t)^{T}I\,\phi(s)\delta_{k,\ell}.$
(23)

In the following, we prove that
$\sum_{j=1}^{N}\left(I+H^{(k,j)}\right)\left(I-F^{(j,\ell)}\right)=I$.

######  Lemma 3.1.

Let $D$ be an infinite matrix so that $\phi(t)^{T}D\phi(s)=\delta(t-s)$, then
$D=I$.

###### Proof.

First of all, since $\delta(t-s)$ is a generalized function, the convergence
of the series $\phi(t)^{T}D\phi(s)$ is intended in a weak sense. This means
that, for every $\tilde{f}(t)$ analytic on $\mathcal{I}$,

$\lim_{N\rightarrow\infty}\int_{\mathcal{I}}\sum_{k,\ell=1}^{N}d_{k,\ell}p_{k}(\tau)p_{\ell}(s)\tilde{f}(\tau)d\tau=\tilde{f}(s)=\int_{\mathcal{I}}\delta(\tau-s)\tilde{f}(\tau)d\tau.$

Setting $\tilde{f}(t)=p_{j}(t)$ gives

$\displaystyle p_{j}(s)$
$\displaystyle=\lim_{N\rightarrow\infty}\sum_{k,\ell=1}^{N}d_{k,\ell}\int_{\mathcal{I}}p_{k}(\tau)p_{\ell}(s)p_{j}(\tau)d\tau$
$\displaystyle=\lim_{N\rightarrow\infty}\sum_{k,\ell=1}^{N}d_{k,\ell}p_{\ell}(s)\int_{\mathcal{I}}p_{k}(\tau)p_{j}(\tau)d\tau,=\sum_{\ell=1}^{\infty}d_{j,\ell}p_{\ell}(s).$

As the Legendre expansion of $p_{j}(s)$ is unique,
$d_{j,\ell}=\delta_{j,\ell}$, for $j,\ell=1,2,\dots$ . ∎

######  Theorem 3.2.

Let $A(t,s)$ be an $N\times N$ matrix-valued function composed of elements
from $\mathcal{D}_{0}^{t}(\mathcal{I})$ and let
$\mathbf{A}=[F^{(k,\ell)}]_{k,\ell=1}^{N}$ be the related block matrix, with
$F^{(k,\ell)}$ the coefficient matrix of $A_{k,\ell}(t,s)$. Moreover, consider
the matrix-valued function $C(t,s)=\sum_{k\geq 1}A^{\star k}(t,s)$, and let
$\mathbf{C}=[H^{(k,\ell)}]_{k,\ell=1}^{N}$ be the related block matrix, with
$H^{(k,\ell)}$ the coefficient matrix of $C_{k,\ell}(t,s)$. Then

$(I+\mathbf{C})(I-\mathbf{A})=I,$

that is, $(I-\mathbf{A})$ is invertible.

###### Proof.

Equations (21), (22) and (23) show that

$\displaystyle\phi(t)^{T}\sum_{j=1}^{N}\left(I+H^{(k,j)}\right)\left(I-F^{(j,\ell)}\right)\phi(s)=\phi(t)^{T}D^{(k,\ell)}\phi(s)=\delta(t-s)\delta_{k,\ell},$

where the products $(I+H^{(k,j)})(I-F^{(j,\ell)})$ are well-defined since
$(I+H^{(k,j)})_{m,n}$ is bounded by $1+K_{h}\sqrt{2m+1}\sqrt{2n+1}$ and
$(I-F^{(j,\ell)})_{m,n}$ by $K_{f}\rho^{|k-\ell|}$, see (15) and (18). Thus,
for $k=\ell$,
$D^{(k,k)}=\sum_{j=1}^{N}\left(I+H^{(k,j)}\right)\left(I-F^{(j,\ell)}\right)=I$
by Lemma 3.1. For $k\neq\ell$, $D^{(k,\ell)}=0$. ∎

Consider the solution $U(t,s)$ of (6) and let us define $\mathbf{T}$ as the
block diagonal matrix with blocks all equal to $T$, the coefficient matrix of
$\Theta(t-s)$. Using the notation of Theorem 3.2 and Expression (10), we can
transform the ODE (6) into the matrix problem:

$\mathbf{U}=\mathbf{T}(I+\mathbf{C})(I+\mathbf{B})=\mathbf{T}(I-\mathbf{A})^{-1}(I+\mathbf{B}),$
(24)

where the block matrix $\mathbf{B}=[G^{(k,\ell)}]_{k,\ell=1}^{N}$ is composed
of the coefficient matrices $G^{(k,\ell)}$ of the functions $B_{k,\ell}(t,s)$.
Hence, the solution of (6) can be expressed by

$U_{k,\ell}(t,s)=\phi(t)^{T}Y^{(k,\ell)}\,\phi(s),\quad k,\ell=1,\dots,N,$

with $\mathbf{U}=[Y^{(k,\ell)}]_{k,\ell=1}^{N}$.

To conclude the presentation, we need to discuss the convergence of the
expansion (13). Indeed, since $f$ is discontinuous for $t=s$, the expansion
does not converge to $f(t,t)$ and, moreover, it converges only linearly for
$t\neq s$; see, e.g., [17, 18].

######  Lemma 3.3.

Consider $f(t,s)\in\mathcal{A}_{\Theta}$ and the related expansion in
orthonormal shifted Legendre polynomials (13). Then,

$\lim_{N\rightarrow\infty}\sum_{k=0}^{N}\sum_{\ell=0}^{N}f_{k,\ell}\,p_{k}(t)p_{\ell}(s)=\left\\{\begin{array}[]{lc}f(t,s),&t\neq
s\\\ f(t,t)/2,&t=s\end{array}\right.,\quad t,s\in(0,1).$

###### Proof.

The proof is direct consequence of Theorem 1 and Remark 1 in Section 4.7 of
[17]. ∎

However, for the fixed $s=0$, the univariate function
$f(t,0)=\tilde{f}(t,0)\Theta(t-0)=\tilde{f}(t,0)$ is analytic over $[0,1]$.
Therefore, defining $a_{k}=\sum_{\ell=0}^{\infty}(f_{k,\ell}\,p_{\ell}(0))$,
we get the Legendre expansion

$f(t,0)=\sum_{k=0}^{\infty}p_{k}(t)\sum_{\ell=0}^{\infty}f_{k,\ell}\,p_{\ell}(0)=\sum_{k=0}^{\infty}a_{k}p_{k}(t),\quad
t\in[0,1].$

Therefore, the truncated series $\sum_{k=0}^{M}a_{k}p_{k}(t)$ converges
geometrically to $f(t,0)$. As a consequence, in a numerical setting, we can
approximate the function $f(t,0)$ by using $F_{M}$, the principal leading
submatrix of $F$, obtaining the approximation

$f(t,0)\approx\phi_{M}(t)^{T}F_{M}\,\phi_{M}(0),$

with $\phi_{M}(t)$ the first $M$ elements of $\phi(t)$. In this case, we
expect to reach a good enough accuracy for a (relatively) small $M$. Note that
for $s>0$ this is not possible, as we expect the emergence of the Gibbs
phenomenon; see, e.g., [18].

By considering the principal leading submatrix of each of the blocks in
formula (24), for $s=0$, we get the following approximated solution to (2)

$U_{0}(t)\approx(I_{N}\otimes\phi_{M}(t)^{T}T_{M})(I_{M}-\mathbf{A}_{M})^{-1}(I_{M}+\mathbf{B}_{M})(I_{N}\otimes\phi_{M}(0));$
(25)

see also [15]. The numerical approach for the solution of a non-autonomous
linear ODE system derived from (25) can be found in [15, 7, 16] where several
numerical examples show its efficacy.

##  4 Conclusion

In this paper, we have presented a new expression for the solution of a (non-
homogeneous, non-autonomous) system of linear ODEs using the so-called
$\star$-product. The $\star$-product, the usual sum, and a specific set of
distributions constitute a ring $\mathcal{S}$. We have also shown that a
certain subring of $\mathcal{S}$ corresponds to a subalgebra of infinite
matrices. Thanks to this correspondence, we have expressed the solution of the
linear ODE system in the infinite matrix algebra. Such a solution is obtained
by inverting a determined infinite matrix. The connection between the
$\star$-product ring and the matrix subalgebra helped us to show that such an
inverse always exists. By truncating the infinite matrices, it is possible to
derive numerical methods for the solution of ODEs. This paper complements the
results we are developing in the truncated case by placing it in the general
framework of the infinite matrix algebra.

## Acknowledgements

This work was supported by Charles University Research programs UNCE/SCI/023
and PRIMUS/21/SCI/009 and by the Magica project ANR-20-CE29-0007 funded by the
French National Research Agency.

## References

  * [1] Autler SH, Townes CH. Stark effect in rapidly varying fields. Phys Rev. 1955 Oct;100:703–722. Available from: https://doi.org/10.1103/PhysRev.100.703.
  * [2] Benner P, Cohen A, Ohlberger M, et al. Model reduction and approximation: Theory and algorithms. Philadelphia: SIAM; 2017. Computational Science and Engineering.
  * [3] Blanes S. High order structure preserving explicit methods for solving linear-quadratic optimal control problems. Numer Algorithms. 2015 Jun;69(2):271–290. Available from: https://doi.org/10.1007/s11075-014-9894-0.
  * [4] Kwakernaak H, Sivan R. Linear optimal control systems. Vol. 1. New York: Wiley-interscience; 1972.
  * [5] Lauder M, Knight P, Greenland P. Pulse-shape effects in intense-field laser excitation of atoms. Opt Acta. 1986;33(10):1231–1252. Available from: https://doi.org/10.1080/713821874.
  * [6] Shirley JH. Solution of the Schrödinger equation with a Hamiltonian periodic in time. Phys Rev. 1965 May;138:B979–B987. Available from: https://doi.org/10.1103/PhysRev.138.B979.
  * [7] Pozza S, Van Buggenhout N. A $\star$-product solver with spectral accuracy for non-autonomous ordinary differential equations ; 2022. Available from: https://arxiv.org/abs/2209.15533.
  * [8] Schwartz L. Théorie des distributions. Paris: Hermann; 1978.
  * [9] Giscard PL, Pozza S. Lanczos-like algorithm for the time-ordered exponential: the $*$-inverse problem. Appl Math. 2020;65(6):807–827. Available from: https://doi.org/10.21136/AM.2020.0342-19.
  * [10] Volterra V, Pérès J. Leçons sur la composition et les fonctions permutables. Paris: Éditions Jacques Gabay; 1928.
  * [11] Giscard PL, Lui K, Thwaite SJ, et al. An exact formulation of the time-ordered exponential using path-sums. J Math Phys. 2015;56(5):053503. Available from: https://doi.org/10.1063/1.4920925.
  * [12] Giscard PL, Pozza S. Tridiagonalization of systems of coupled linear differential equations with variable coefficients by a Lanczos-like method. Linear Algebra Appl. 2021;624:153–173. Available from: https://doi.org/10.1016/j.laa.2021.04.011.
  * [13] Giscard PL, Bonhomme C. Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR. Phys Rev Research. 2020 Apr;2:023081. Available from: https://link.aps.org/doi/10.1103/PhysRevResearch.2.023081.
  * [14] Giscard PL, Pozza S. A Lanczos-like method for non-autonomous linear ordinary differential equations. Bol Unione Mat Ital. 2022;Available from: https://doi.org/10.1007/s40574-022-00328-6.
  * [15] Pozza S, Van Buggenhout N. A new matrix equation expression for the solution of non-autonomous linear systems of ODEs ; 2022. Available from: https://arxiv.org/abs/2210.07052.
  * [16] Pozza S, Van Buggenhout N. The *-product approach for linear ODEs: a numerical study of the scalar case ; 2022. Available from: https://arxiv.org/abs/2209.13322.
  * [17] Lebedev N, Silverman R. Special functions and their applications. Dover Publications; 1972. Dover Books on Mathematics.
  * [18] Trefethen LN. Approximation theory and approximation practice. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 2013.


# Softer than soft: diving into squishy granular matter

Jonathan Barés  LMGC, Université de Montpellier, CNRS, Montpellier, France.
Manuel Cárdenas-Barrantes  LMGC, Université de Montpellier, CNRS, Montpellier,
France.  David Cantor  Department of Civil, Geological and Mining Engineering,
Polytechnique Montréal, Montréal, QC, Canada.  Mathieu Renouf  LMGC,
Université de Montpellier, CNRS, Montpellier, France.  Émilien Azéma  LMGC,
Université de Montpellier, CNRS, Montpellier, France.  Institut Universitaire
de France (IUF), Paris, France.

###### Abstract

Softer than soft, squishy granular matter is composed of grains capable of
significantly changing their shape (typically larger than $10$% of
deformation) without tearing or breaking. Because of the difficulty to test
these materials experimentally and numerically, such a family of discrete
systems remains largely ignored in the granular matter physics field despite
being commonly found in nature and industry. Either from a numerical,
experimental, or analytical point of view, the study of highly deformable
granular matter involves several challenges covering, for instance: ($i$) the
need to include a large diversity of grain rheology, ($ii$) the need to
consider large material deformations, and ($iii$) the analysis upon the
effects the large body distortion has on the global scale. In this article, we
propose a thorough definition of these squishy granular systems, and we
summarize the upcoming challenges in their study.

##  1 Introduction

The importance of granular matter for human activities is quietly striking. It
not only covers $90$% of the Earth’s crust [1], it adds up to about $70$% of
the materials used in the industry, represents the second-largest volume of
matter carried by human beings – just behind water [2] –, and more energy is
used to produce it than that consumed in human transportation [3]. For all
these reasons, a thorough understanding of the behavior of granular matter
made of all sorts of materials is of tremendous importance for humankind, from
environmental preservation and control of geohazards to space exploration.
Yet, this is an extremely complex task since these materials span over several
scales, from colloids to screes or asteroids. These materials can also be
quite different – made of virtually any sort of matter – they can behave like
a solid, a liquid, or a gas depending on loading conditions and local
interactions. Over the last three decades, a framework and state equations to
describe the flow of granular systems made of rigid, virtually undeformable,
particles have been set up [4, 5, 6]. Similarly, the jamming transition of
these rigid grains has been fruitfully explored and the so obtained results
have been extended to the case of soft grains [7, 8, 9, 10]. Nevertheless, it
is unsettling that in much of the literature involving deformable granular
media the material these particles are made of is assumed to stay in the
linear elastic regime so that the Hertz contact assumption remains valid [11,
12]. This means that the level of this material deformation does not exceed
$10$%. However, we can naturally wonder about the consequences of such an
assumption in granular systems undergoing more intense loading, well above
$10$% of deformation [13]. The existing framework mostly neglects the grains’
change in shape, limiting the number of real-life scenarios that can be
described: This approach is good at describing the loading of corn grains but
not of pop-corn.As shown in Fig. 1, many systems, ranging from biology to
industrial applications, are made of particles brought to large deformation
regimes. To differentiate these systems from those composed of soft grains, we
propose to call squishy granular materials all discrete matter whose
individual bodies undergo large nonlinear deformations. This definition
implies that the global stresses undergone by the system are substantially
higher than the Young modulus of the particles, i.e.,
$E_{\rm{material}}\ll\sigma_{\rm{loading}}$.

Figure 1: Squishy granular systems. (a) Highly deformed 2D soap bubbles from
[14]. (b) Sintered ceramic from [15]. (c) Jammed layer of biologic cells from
[16]. (d) Plastic particles made of agar gel deep in the jammed state from
[17].

In addition, stating that a granular system is squishy also means that the
particles’ change of shape must be considered to describe their local
interactions, displacement, and, more generally, to understand the global
behavior of the granular medium. In other words, the exact nature of a squishy
material cannot be hidden behind a linear elastic law. The intrinsic rheology
of the particles becomes paramount and must be considered through the
nonlinear theory of elasticity [18]. One step further, the exact geometry of
the particles must be thoroughly considered since it can, for example, induce
discontinuous and sudden phenomena such as buckling. This is the case for
shell particles [19, 20] or, more generally, for particles made of any
metamaterial [21, 22]. Naturally, a large diversity of geometries and material
behaviors have to be considered to properly handle squishy media. This can
encompass thermal and anisotropic effects, as well as magnetic interactions
[23, 24]. As a limitation in the broad category of squishy particles, we
assume that the number of entities in the system stays stable all along the
loading process. The squeezable nature of the particles allows them to support
any loading with no merging or fragmentation. Thus, the edge of each particle
is continuously described by a closed function.

Fortunately, the granular matter physics field is not starting from zero as
far as it concerns squishy particles. The most advanced studies in this domain
come from the rheology of liquid foams [25], materials composed of gas bubbles
within liquids. Numerous experimental and numerical works have significantly
contributed to the understanding of these materials under different loadings
[26, 27, 28, 14] and direct links have been drawn to granular matter physics
[29]. However, emulsions and colloids wildly differ from other discrete
materials when considering their material rheology and inter-particle
frictional properties [30, 31, 32, 33]. Consequently, it is difficult to
generalize the results found for foams, colloids, and emulsions to the case of
squishy granular matter. Motivated by technological applications, very soft
grains have also been studied under compression as in the case of sintered
material manufacturing. The compaction of ceramic, metal, gels, and drug
powders has been central in the research of compacted materials [34, 35, 36]
or self-formed colloids [37]. In all these cases, however, studies are still
limited to the plastic behavior of particles under compression. Other examples
of squishy grains can be found in living organisms. Cells forming a tissue or
clogging in a vein can certainly be described using squishy particle
assemblies [31, 38, 39, 40, 16].

From a numerical point of view, only recently a number of models have emerged
for squishy granular matter. Some of them, based on the discrete element
method (DEM), have allowed researchers to catch a glimpse of the behavior of
multiparticle, highly deformable systems [41, 20, 42], but remain
unsatisfactory in terms of material mimicking. Other methods based on finite-
elements (FEM) are better at mimicking particles’ materials, but have
difficulties to properly take contact interactions into account [43, 44, 45,
46]. In both cases, major improvements are necessary to simulate any squishy
granular system.

The rest of this article is dedicated to summarizing the different loading
conditions that squishy granular matter can undergo and the different
approaches to deal with their characterization, simulation, and analysis. More
particularly, we focus on the experimental and numerical challenges to
understand the behavior of squishy grains under compression and shearing.

##  2 Browsing deep in the jammed state

###  2.1 Compressing squishy particle assemblies

The mechanical behavior of a granular system is ruled by the individual
behavior of each particle within it. For example, when hard particles – i.e.,
made of a material with high stiffness relatively to the applied stress – are
compressed inside a box, they tend to rearrange compactly. In this case, the
rearrangement ends when a stable force network is reached. This state, called
the jammed state, depends mainly on the morphology of the particles, their
size and shape [47], and their interaction laws [48, 47, 49]. If one considers
assemblies of soft particles with low stiffness relatively to the applied
stress, they are able to undergo large deformations, and the compression can
continue beyond the jammed state. Exploring the compaction of granular matter
and, especially, systems very deep in the jammed state is challenging, mainly
due to the nature of the particle rearrangements and the evolution of the
force-chain network. For assemblies of squishy particles, the compactness –
i.e., the ratio between the volume of the particles and the volume of the box
– surpasses the compactness of the jammed state, the random close packing, and
tends to a maximal value close to $1$. Naturally, this high compactness comes
with an increasing number of contacts and contact surfaces. Therefore, it is
seen as an evolution of the probability distribution of the forces and an
evident reorganization of the force-chain network that drastically change the
mechanical properties of the assembly.

Innovative experiments and advanced numerical methods have allowed one to take
a small step forward in understanding the microstructural evolution beyond the
jamming point (see Section 3). However, theoretical modeling of the compaction
process remains a complicated task. Although many equations have been proposed
in the last decades to try to link the confining pressure to the compactness
[50, 51, 52, 53], most of them are based on empirical strategies with no
apparent physical bases, requiring several fitting parameters. Only recently,
a systematic approach, based entirely on micromechanical parameters and free
of any ad hoc parameters, has been proposed and numerically tested in the case
of hyperelastic circular, polygonal, and spherical shapes under different
frictionnal conditions [45, 46]. This method was derived from the
micromechanical expression of the granular stress tensor together with the
approximation of the Hertz contact law between two particles at small strain
and the theoretical solution of the hydro-static compaction of a single
particle.

Although the above micromechanical approach constitutes a substantial progress
in the challenge of understanding and modeling the compaction of deformable
particles, it only describes convex particles under isotropic compaction with
continuous laws of elasticity. Many unknowns remain unresolved concerning
particles with a more complex shape or composition characteristics, such as
non-convex particles or particles that do not follow elastic or continuous
deformation behavior (for example, particles suffering buckling). Also, we
still need to better understand the behavior of granular systems under
different loading conditions, such as uniaxial compaction or shear loading.

###  2.2 Shearing squishy particle assemblies

The second and more complex loading that a squishy granular system can undergo
is shearing. This loading type induces not only large particle deformations
but also large and often erratic particle rearrangements. When this latter
regime dominates, the granular system yields and flows. Below the jammed
state, for rigid grains, the flowing regime is described in many
configurations by well-known state equations [4, 5, 6]. Although squishy
grains could present state equations reminiscent of those for rigid grains,
this has not been shown yet. So far, for a dense flow of soft grains, it has
been shown that grains rearrange through a succession of local mechanisms
called T1-events [14, 54] inducing specific glassy rheology. These T1-events,
presented in Fig.2(e), mainly happen in a relatively narrow area of the system
called the shear band. This has been observed for liquid foams [14] and living
cells [54], but one can still wonder how the shear-band evolves depending on
the material rheology. Do these T1-events still exist for other material types
such as plastic ones?

The shearing of granular systems often exhibits erratic dynamic mechanisms at
the global scale – also known as avalanches – even under constant loading [55,
56]. For soft grains, depending on the particle stiffness, strain rate, and
confining pressure, the global dynamics can display different statistical
signatures [57]. In densely packed squishy systems, the space for particles to
rearrange can be extremely reduced and the material rheology has a strong
influence on the possible rearrangement mechanisms. As such, it is expected
that the avalanche dynamics are strongly modified by the exact particle nature
and frictional properties. But we still do not know how or up to which extent
this is the case.

In the case of rigid or soft particles, plasticity, seen as an irreversible
deformation of the material, happens at the meso-scale: grains rearrange so
that they cannot recover their previous position when the loading is released.
In the case of squishy granular systems not only this plastic mechanism is at
play, but also the intrinsic plasticity of the particles. Hence, the
plasticity of these amorphous materials is by definition of a different nature
as the one of soft grains.

Concerning the shearing of squishy particles, the main challenge is to
evidence the local rearrangement mechanisms, their dynamics and spatial
distribution, and to understand how they induce the flowing of the sheared
system. The second challenge consists in understanding the coupling between
meso and micro-scale plasticity processes in squishy amorphous materials

##  3 Experimental & Numerical challenges

###  3.1 Physical experiments

Studying experimentally granular matter is challenging because its multi-scale
nature, so tests have to be able to catch main features across scales. One
needs to follow the behavior of individual grains that end up inducing the
global material behavior from specific particle properties. In the case of
squishy granular matter, one also needs to go one scale smaller since it is
necessary to consider the material behavior inside each grain at the sub-
granular scale [17]. Although several studies are possible without analyzing
this sub-granular level or measuring the local stresses – as it has been done
in the early stages of the granular matter physics – a comprehensive study of
squishy granular matter requires such a multiscale approach.

So far, whether it is in a 2D or 3D configuration, the existing ways to
measure grain displacements, contacts, and stresses are all based on an
inverse problem method. Usually, the Hertz contact law [11, 12], or a certain
surface tension in the case of liquid foams [14], is assumed and some local
field measurements permit to detect contacts and guess interaction forces.
These local measurements can be either the grain boundary deformation [30, 14,
58], a photoelastic [59, 60] or thermoelastic [61] signal, or a digital image
correlation (DIC) rough measurements [62]. By definition, these measurements
are limited by the method sensitivity and the domain of validity of the
assumed contact law. For example in the case of photoelasticity, when
deformation is too large, interference fringes overlap and the measured signal
saturates. Assuming that there is no saturation of the measured signal for
large deformations, it would be possible to extend these methods by finding
contact laws that include complex material rheologies. For example, if the
Hertz contact law could be extended beyond the case of a linear elastic
behavior with a finite friction coefficient [13], this would permit extending
these methods to squishy grains. Still, one can imagine that the number of
parameters to guess in the problem inversion would be much larger, making it
complicated to implement.

Recently it has been proposed not to make any assumption and to directly
measure displacement fields at the grain scale, even in the large deformation
regime. This process is based on a large field and an accurate imaging method
[13, 17, 63], which allows allows one to measure all the tensors of the non-
linear elasticity theory and to guess the stress tensor knowing the loading
history and the material rheology via large deformation DIC. However, this
only works for small systems ($\sim 100$ grains) under quasistatic conditions.
This method – and the fact that it provides strain fields inside the grains –
comes with a new outlook on this topic upon the notion of force chain and its
relevance or validity in the study of densely packed configuration.

Although recent experimental approaches have permitted to follow the evolution
of squashy granular systems, the development of an experimental method to
characterize squishy particle assemblies under dynamic loading remains a
challenging task.

###  3.2 Numerical simulations

Figure 2: (a to d) Modeling squishy particle assembly: (a) using the
“Deformable-Particle” (DP) method, a class of BPM-like method (with permission
of C. O’hern et al. [20]), (b) using the meshless approach MPM (with
permission of S. Nezamabadi et al. [41]), (c) using the meshless “soft-DEM”
approach (with permission of G. Mollon [64]), (d) using the 3D Non-Smooth
Contact Dynamics approach by M. Cardenas-Barrantes et al. [46]. (e) Example of
a T1-event running from time $t_{1}$ to $t_{3}$. Particles react to shear
loading by this peculiar rearrangement: particles in the expansion direction
split apart and permit the ones in the compression direction to get closer up
to the point where they touch each other.

The modeling of granular media is historically associated with discrete-
element approaches (DEM) to describe the dynamical evolution of a collection
of hard grains. DEM can be split into two main approaches: the ($i$) smooth
and ($ii$) non-smooth methods. This classification is related to the way a
contact between two particles is considered. In the case ($i$), regularised
contact laws link the contact force to a slight overlap between particles in
contact. In the case ($ii$), the unilaterality of contacts (non-regularized
laws) traduces the non-penetration of the particles. Both approaches have
proven to produce equivalent results as far as the basic hypotheses of each
method are respected – i.e., binary contact and small grain deformations
compared to the grain size.

Extending the numerical approaches, first to soft, and further to squishy
grains, implies new assumptions and new numerical strategies. In the framework
of the smooth approach, one alternative is to increase the allowed overlapping
between bodies that can be later interpreted as a virtual grain deformation.
Although this approach is debatable in terms of realistic representation of
the grain shapes, it has shown to be valid in the small deformation domain, as
in the case of the compaction of hard grains assemblies just after reaching
the jammed state [65, 66, 67]. Nevertheless, going deeper into the jammed
state leads to the breaking of some of the model assumptions and the classical
methods need to be extended to meet this new range of deformation.

As a first approximation, a squishy grain can be considered as an ‘aggregate’
of elements (in the numerical sense) bonded together by bilateral or cohesive
bonds. The bonding law can be adjusted to mimic the macroscopic behaviour of a
given grain material. This aggregate representation can be performed
throughout the volume of the grain. This approach is known as the Bonded-
Particle Method (BPM) [68, 69, 70, 71]. Another alternative is to consider
sub-elements only on the particle’s surface, as proposed in the Deformable
Particle Method (DPM) [72, 20] (see Fig. 2-a). The main advantages of these
methods are their low computational cost and ease of implementation in the DEM
framework. However, in the case of ’aggregate’ methods, very large
deformations can lead to non-physical situations. Another option consists in
using constitutive models to describe the bulk behavior of the particles.
Then, it is possible to use standard numerical methods such as a finite
elements (FEM) or meshless approaches.

In a meshed approach, a volumetric rheological model associated with a state
equation is considered for each grain. In this case, FEM-like methods coupled
with DEM are used to describe the particle deformations, keeping – or not –
its integrity [73, 74, 75, 76]. The advantage of these strategies is to
closely mimic the geometry and bulk properties of the particles. However, they
end up being computationally expensive. As with the classical DEM approaches,
it is possible to distinguish two classes of methods: the Multi-Particle
Finite Element Method (MPFEM) [77, 78, 73, 79], based on regularized contact
interactions, and the Non-Smooth Contact Dynamic Method (NSCD) [80, 81], based
on non-regularized contact laws. MPFEM simulations have been performed by
Mesarovic _et al._ [82] to study the deformation mechanisms in a periodic
configuration of dissimilar elastoplastic spheres. For squishy granular
packings, the first MPFEM studies were conducted by Gethin et al. in 2002 [77]
and Chen et al. in 2006 [83] in 2D and 3D, in the context of metal compaction.
The NSCD was used for the first time in 2000 by Acary et al. [84] with a work
on masonry. They were followed by Cao et al. in 2011 [85], who tested the
numerical approach for the compaction and shearing of soft packings. The
latest works related to compaction of squishy granular packings using NSCD
corresponds to Vu et al. [44, 86] in 2D, followed by Cardenas et al. in 2D and
3D [46] (see Fig. 2-d).

In the framework of mesh-free models, squishy grains have been simulated via
the Material Point Method (MPM) [41] (see Fig. 2(b)), or the Soft Discrete
Element Method (SDEM) [64] (see Fig.2(c)), where simple kinematics are
postulated in order to represent the ovalisation of the grains and their local
deformations around interparticle contacts. In the last case, a continuous
deformation field is interpolated between the nodes belonging to a grain by
means of a moving least square algorithm. These approaches require the
calibration of some fundamental computational parameters from experimental
data, such as the stiffness and viscosity [87, 64]. Despite the development of
these alternative methods, the different results obtained via FEM-DEM compare
better to experiments and show that the behavior of squishy grain collections
can only be understood using these advanced numerical approaches despite being
extremely time-consuming.

Most of today’s simulations involving squishy grains are performed in the
quasi-static regime and under uniaxial deformation. These conditions remain
relatively simple and involve few dynamical effects. Other loading conditions
such as shearing would rapidly increase the complexity of the simulations,
making it difficult to conduct numerical experiments at higher shear rates or
compression levels. Moreover, while several thousands of squishy grains can be
modeled in 2D, only a few hundred can be simulated in 3D in a reasonable
computing time. Therefore, additional challenges concern ($i$) the development
of scalable and highly parallel algorithms, ($ii$) the implementation of fully
periodic conditions, allowing for the volume changes to account for the
expansion/contraction of the particles, and ($iii$) considering realistic
mixtures that include grains with different size, shape, and bulk behavior
within the same system.

##  4 Summary & Perspectives

Deep in the jammed state, when grain deformation is well above $10$%, the soft
granular matter framework is not relevant to describe the evolution of
granular systems anymore. As packing fraction grows way beyond the jamming
point, it is necessary to consider the large particle change in shape and it
effects: ($i$) creation of new contacts, ($ii$) interaction stresses depending
on the particle material, and ($iii$) possible large particle rearrangements.
To make up for the lack of a correct framework for these systems, we propose
the concept of squishy granular materials in which the linear elastic behavior
is overpassed to consider hyperelastic, plastic or other complex material
rheologies.

Considering squishy materials rapidly increases the complexity of the
experimental, numerical and theoretical approaches needed for their study.
This involves several technical and scientific challenges including the
understanding of the local deformation mechanisms that lead to a global
compaction law. Such a law is expected to include the effect of variability in
the particles’ rheology. Concerning the shear process of squishy systems, the
main challenge consists in getting insights into the local rearrangement
mechanisms, generalizing the T1-event concept. This should work as a building
block for the description of the deformation of amorphous squishy materials.
Studying the dynamics of this phenomenon and the spatial distribution of these
events will also constitute a major step forward in understanding the plastic
processes in these systems.

These challenges cannot be overcome without tackling serious numerical and
experimental difficulties. The first one consists in finding an experimental
method to measure the evolution of dynamically loaded 2D and 3D systems down
to the sub-granular scale. The numerical counterpart of this task is equally
challenging since simulation methods can currently mainly deal with the
quasistatic compression of a limited number of grains. The complexity of the
simulations can rapidly increase as soon as shearing tests are considered, or
realistic and strongly deformed particles geometries ought to be implemented.

Setting up a new framework to analyse the behavior of squishy particle systems
deep in the jammed state is an open and complex task yet to be developed in
granular matter research. Some aspects concerning the experimental and
numerical study of squishy matter are currently maturing, but major challenges
are still ahead. Fortunately, these systems are so ubiquitous that scientists
brave enough to take up these challenges can count on the support of a broad
and diverse community gathering physicists, geophysicists, biologists and
manufacturers to name a few.

## References

  * [1] H. C. Nataf and J. Sommeria.  La physique de la terre.  Belin, 2000.
  * [2] F. Lechenault.  La transition de” Jamming” dans un milieu granulaire bidimensionnel: Statique et dynamique d’un système athermique modèle.  PhD thesis, Université Paris Sud-Paris XI, 2007.
  * [3] É. Guyon, J.-Y. Delenne, and F. Radjaï.  Built on Sand: The Science of Granular Materials.  Mit Press, 2020.
  * [4] GDR MiDi.  On dense granular flows.  The European Physical Journal E, 14:341–365, 2004.
  * [5] P. Jop, Y. Forterre, and O. Pouliquen.  A constitutive law for dense granular flows.  Nature, 441(7094):727–730, 2006.
  * [6] O. Coquand and M. Sperl.  Rheology of granular liquids in extensional flows: Beyond the $\mu$ (i)-law.  Physical Review E, 104(1):014604, 2021.
  * [7] A. J. Liu and S. R. Nagel.  Jamming is not just cool any more.  Nature, 396(6706):21–22, 1998.
  * [8] C. S. O’hern, L. E. Silbert, A. J. Liu, and S. R. Nagel.  Jamming at zero temperature and zero applied stress: The epitome of disorder.  Physical Review E, 68(1):011306, 2003.
  * [9] D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer.  Jamming by shear.  Nature, 480(7377):355–358, 2011.
  * [10] Y. Zhao, J. Barés, H. Zheng, J. E.S. Socolar, and R. P. Behringer.  Shear-jammed, fragile, and steady states in homogeneously strained granular materials.  Physical review letters, 123(15):158001, 2019.
  * [11] J. F. Archard.  Elastic deformation and the laws of friction.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 243(1233):190–205, 1957.
  * [12] K. L. Johnson.  Contact mechanics.  1985\.
  * [13] T.-L. Vu, J. Barés, S. Mora, and S. Nezamabadi.  Deformation field in diametrically loaded soft cylinders.  Experimental Mechanics, 59(4):453–467, 2019.
  * [14] A. Kabla, J. Scheibert, and G. Debregeas.  Quasi-static rheology of foams. part 2. continuous shear flow.  Journal of Fluid Mechanics, 587:45–72, 2007.
  * [15] W. H. Rhodes.  Phase chemistry in the development of transparent polycrystalline oxides.  Phase Diagrams for Advanced Ceramics, pages 1–42, 1995.
  * [16] M. J. O’Sullivan, T.-K. N. Phung, and J.-A. Park.  Bronchoconstriction: a potential missing link in airway remodelling.  Open Biology, 10(12):200254, 2020.
  * [17] T.-L. Vu and J. Barés.  Soft-grain compression: Beyond the jamming point.  Physical Review E, 100(4):042907, 2019.
  * [18] L. A. Taber.  Nonlinear theory of elasticity: applications in biomechanics.  World Scientific, 2004.
  * [19] J. Jose, G. A. Blab, A. Van Blaaderen, and A. Imhof.  Jammed elastic shells–a 3d experimental soft frictionless granular system.  Soft Matter, 11(9):1800–1813, 2015.
  * [20] A. Boromand, A. Signoriello, F. Ye, C. S. O’Hern, and M. D. Shattuck.  Jamming of deformable polygons.  Physical review letters, 121(24):248003, 2018.
  * [21] B. Florijn, C. Coulais, and M. van Hecke.  Programmable mechanical metamaterials.  Physical review letters, 113(17):175503, 2014.
  * [22] Pedro M Reis, Heinrich M Jaeger, and Martin Van Hecke.  Designer matter: A perspective.  Extreme Mechanics Letters, 5:25–29, 2015.
  * [23] D. Hernández-Enríquez, G. Lumay, and F. Pacheco-Vázquez.  Discharge of repulsive grains from a silo: experiments and simulations.  In EPJ Web of Conferences, volume 140, page 03089. EDP Sciences, 2017.
  * [24] M. Cox, D. Wang, J. Barés, and R. P. Behringer.  Self-organized magnetic particles to tune the mechanical behavior of a granular system.  EPL (Europhysics Letters), 115(6):64003, 2016.
  * [25] R. Höhler and S. Cohen-Addad.  Rheology of liquid foam.  Journal of Physics: Condensed Matter, 17(41):R1041, 2005.
  * [26] G. Debregeas, H. Tabuteau, and J.-M. Di Meglio.  Deformation and flow of a two-dimensional foam under continuous shear.  Physical Review Letters, 87(17):178305, 2001.
  * [27] M. Asipauskas, M. Aubouy, J. A. Glazier, F. Graner, and Y. Jiang.  A texture tensor to quantify deformations: the example of two-dimensional flowing foams.  Granular Matter, 5(2):71–74, 2003.
  * [28] R. J. Clancy, E. Janiaud, D. Weaire, and S. Hutzler.  The response of 2d foams to continuous applied shear in a couette rheometer.  The European Physical Journal E, 21(2):123–132, 2006.
  * [29] D. Weaire, V. Langlois, M. Saadatfar, and S. Hutzler.  Foam as granular matter.  In Granular and Complex Materials, pages 1–26. World Scientific, 2007.
  * [30] J. Brujić, S. F. Edwards, D. V. Grinev, I. Hopkinson, D. Brujić, and H. A. Makse.  3d bulk measurements of the force distribution in a compressed emulsion system.  Faraday discussions, 123:207–220, 2003.
  * [31] K. A. Newhall, L. L. Pontani, I. Jorjadze, S. Hilgenfeldt, and J. Brujic.  Size-topology relations in packings of grains, emulsions, foams, and biological cells.  Physical review letters, 108(26):268001, 2012.
  * [32] T. Krebs, D. Ershov, C. GPH. Schroen, and R. M. Boom.  Coalescence and compression in centrifuged emulsions studied with in situ optical microscopy.  Soft Matter, 9(15):4026–4035, 2013.
  * [33] R. Höhler and D. Weaire.  Can liquid foams and emulsions be modeled as packings of soft elastic particles?  Advances in colloid and interface science, 263:19–37, 2019.
  * [34] A. R. Cooper Jr and L. E. EATON.  Compaction behavior of several ceramic powders.  Journal of the American Ceramic Society, 45(3):97–101, 1962.
  * [35] K. Kawakita and K.-H. Lüdde.  Some considerations on powder compression equations.  Powder technology, 4(2):61–68, 1971.
  * [36] C.-Y. Wu, O. M. Ruddy, A. C. Bentham, B. C. Hancock, S. M. Best, and J. A. Elliott.  Modelling the mechanical behaviour of pharmaceutical powders during compaction.  Powder technology, 152(1-3):107–117, 2005.
  * [37] J. Fan, S.-H. Kim, Z. Chen, S. Zhou, E. Amstad, T. Lin, and D. A. Weitz.  Creation of faceted polyhedral microgels from compressed emulsions.  small, 13(31):1701256, 2017.
  * [38] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, et al.  Unjamming and cell shape in the asthmatic airway epithelium.  Nature materials, 14(10):1040–1048, 2015.
  * [39] T. PJ. Wyatt, A. R. Harris, M. Lam, Q. Cheng, J. Bellis, A. Dimitracopoulos, A. J. Kabla, G. T. Charras, and B. Baum.  Emergence of homeostatic epithelial packing and stress dissipation through divisions oriented along the long cell axis.  Proceedings of the National Academy of Sciences, 112(18):5726–5731, 2015.
  * [40] J. Mauer, S. Mendez, L. Lanotte, F. Nicoud, M. Abkarian, G. Gompper, and D. A. Fedosov.  Flow-induced transitions of red blood cell shapes under shear.  Physical review letters, 121(11):118103, 2018.
  * [41] S. Nezamabadi, T. H. Nguyen, J.-Y. Delenne, and F. Radjai.  Modeling soft granular materials.  Granular Matter, 19(1):1–12, 2017.
  * [42] D. Wang, J.D. Treado, A. Boromand, B. Norwick, M. P. Murrell, M. D. Shattuck, and C. S. O’Hern.  The structural, vibrational, and mechanical properties of jammed packings of deformable particles in three dimensions.  Soft Matter, 17(43):9901–9915, 2021.
  * [43] G. Mollon.  Mixtures of hard and soft grains: micromechanical behavior at large strains.  Granular Matter, 20(3):1–16, 2018.
  * [44] T.-L. Vu, J. Barés, S. Mora, and S. Nezamabadi.  Numerical simulations of the compaction of assemblies of rubberlike particles: A quantitative comparison with experiments.  Physical Review E, 99(6):062903, 2019.
  * [45] D. Cantor, M. Cárdenas-Barrantes, I. Preechawuttipong, M. Renouf, and E. Azéma.  Compaction model for highly deformable particle assemblies.  Physical Review Letters, 124(20):208003, 2020.
  * [46] Manuel Cárdenas-Barrantes, David Cantor, Jonathan Barés, Mathieu Renouf, and Emilien Azéma.  Three-dimensional compaction of soft granular packings.  Soft Matter, 18:312–321, 2022.
  * [47] M. van Hecke.  Jamming of soft particles: geometry, mechanics, scaling and isostaticity.  Journal of Physics: Condensed Matter, 22(3):033101, 2009.
  * [48] S. Torquato and F. H. Stillinger.  Jammed hard-particle packings: From kepler to bernal and beyond.  Rev. Mod. Phys., 82:2633–2672, Sep 2010.
  * [49] Leonardo E. Silbert.  Jamming of frictional spheres and random loose packing.  Soft Matter, 6:2918–2924, 2010.
  * [50] R W Heckel.  Density-Pressure Relations in Powder Compaction.  Trans Metal Soc Aime, 221(4):671–675, 1961.
  * [51] K.T. Kim and M.M. Carroll.  Compaction equations for strain hardening porous materials.  International Journal of Plasticity, 3(1):63–73, 1987.
  * [52] J. Secondi.  Modelling powder compaction: From a pressure-density law to continuum mechanics.  Powder Metallurgy, 45(3):213–217, 2002.
  * [53] Ramon Cabiscol, Hao Shi, Isabell Wünsch, Vanessa Magnanimo, Jan Henrik Finke, Stefan Luding, and Arno Kwade.  Effect of particle size on powder compaction and tablet strength using limestone.  Advanced Powder Technology, 2020.
  * [54] D. Bi, J. H. Lopez, J. M. Schwarz, and L. M. Manning.  Energy barriers and cell migration in densely packed tissues.  Soft matter, 10(12):1885–1890, 2014.
  * [55] K. E. Daniels and N. W. Hayman.  Force chains in seismogenic faults visualized with photoelastic granular shear experiments.  Journal of Geophysical Research: Solid Earth, 113(B11), 2008.
  * [56] J. Barés, D. Wang, D. Wang, T. Bertrand, C. S. O’Hern, and R. P. Behringer.  Local and global avalanches in a two-dimensional sheared granular medium.  Physical Review E, 96(5):052902, 2017.
  * [57] A. Abed-Zadeh, J. Barés, and R. P. Behringer.  Crackling to periodic dynamics in granular media.  Physical Review E, 99(4):040901, 2019.
  * [58] N. Brodu, J. A. Dijksman, and R. P. Behringer.  Spanning the scales of granular materials through microscopic force imaging.  Nature communications, 6(1):1–6, 2015.
  * [59] K. E. Daniels, J. E. Kollmer, and J. G. Puckett.  Photoelastic force measurements in granular materials.  Review of Scientific Instruments, 88(5):051808, 2017.
  * [60] A. Abed-Zadeh, J. Barés, T. A. Brzinski, K. E. Daniels, J. Dijksman, N. Docquier, H. O. Everitt, J. E. Kollmer, O. Lantsoght, D. Wang, et al.  Enlightening force chains: a review of photoelasticimetry in granular matter.  Granular Matter, 21(4):1–12, 2019.
  * [61] P. Jongchansitto, X. Balandraud, M. Grédiac, C. Beitone, and I. Preechawuttipong.  Using infrared thermography to study hydrostatic stress networks in granular materials.  Soft matter, 10(43):8603–8607, 2014.
  * [62] R. Hurley, E. Marteau, G. Ravichandran, and J. E. Andrade.  Extracting inter-particle forces in opaque granular materials: beyond photoelasticity.  Journal of the Mechanics and Physics of Solids, 63:154–166, 2014\.
  * [63] M. Cárdenas-Barrantes, J. Barés, M. Renouf, and E. Azéma.  Experimental validation of a micromechanically-based compaction law for soft/hard grain mixtures.  arXiv preprint arXiv:2111.01568, 2021.
  * [64] G. Mollon.  The soft discrete element method.  Granular Matter, 24(1):1–20, 2022.
  * [65] Ivana Agnolin and Jean Noël Roux.  On the elastic moduli of three-dimensional assemblies of spheres: Characterization and modeling of fluctuations in the particle displacement and rotation.  International Journal of Solids and Structures, 45:1101 – 1123, 2008.
  * [66] Roger T Bonnecaze and Michel Cloitre.  Micromechanics of soft particle glasses.  High Solid Dispersions, pages 117–161, 2010.
  * [67] J. Lopera-Perez, C. Kwok, and K. Senetakis.  Micromechanical analyses of the effect of rubber size and content on sand-rubber mixtures at the critical state.  Geotextiles and Geomembranes, 45(2):81–97, 2017.
  * [68] D.O. Potyondy and P.A. Cundall.  A bonded-particle model for rock.  International Journal of Rock Mechanics and Mining Sciences, 41(8):1329–1364, 2004.
  * [69] S. Utili and R. Nova.  Dem analysis of bonded granular geomaterials.  International Journal for Numerical and Analytical Methods in Geomechanics, 32(17):1997–2031, 2008.
  * [70] N. Cho, C.D. Martin, and D.C. Sego.  A clumped particle model for rock.  International Journal of Rock Mechanics and Mining Sciences, 44(7):997–1010, 2007.
  * [71] Mohsen Asadi, Ahmad Mahboubi, and Klaus Thoeni.  Discrete modeling of sand-tire mixture considering grain-scale deformability.  Granular Matter, 20(2), 05 2018.
  * [72] Y. Chélin, K. Azzag, P. Cañadas, J. Averseng, S. Baghdiguian, and B. Maurin.  Simulation of cellular packing in non-proliferative epithelia.  Journal of Biomechanics, 46(6):1075 – 1080, 2013.
  * [73] Adam T. Procopio and Antonios Zavaliangos.  Simulation of multi-axial compaction of granular media from loose to high relative densities.  Journal of the Mechanics and Physics of Solids, 53(7):1523–1551, 2005.
  * [74] B. Harthong, J.-F. Jérier, P. Dorémus, D. Imbault, and F.-V. Donzé.  Modeling of high-density compaction of granular materials by the discrete element method.  International Journal of Solids and Structures, 46(18):3357–3364, 2009.
  * [75] Fen Huang, Xizhong An, Yuxi Zhang, and A. B. Yu.  Multi-particle FEM simulation of 2D compaction on binary Al/SiC composite powders.  Powder Technology, 314:39–48, 2017.
  * [76] Defeng Wang, Xizhong An, Peng Han, Haitao Fu, Xiaohong Yang, and Qingchuan Zou.  Particulate scale numerical investigation on the compaction of tic-316l composite powders.  Mathematical Problems in Engineering, 2020:5468076, Feb 2020.
  * [77] David T Gethin, Roland W Lewis, and Rajesh S Ransing.  A discrete deformable element approach for the compaction of powder systems.  Modelling and Simulation in Materials Science and Engineering, 11(1):101–114, dec 2002.
  * [78] X.J. Xin, P. Jayaraman, G.S. Daehn, and R.H. Wagoner.  Investigation of yield surface of monolithic and composite powders by explicit finite element simulation.  International Journal of Mechanical Sciences, 45(4):707–723, 2003\.
  * [79] Göran Frenning.  Towards a mechanistic contact model for elastoplastic particles at high relative densities.  Finite Elements in Analysis and Design, 104:56–60, 2015.
  * [80] J. Moreau.  Some numerical methods in multibody dynamics: application to granular materials.  European Journal of Mechanics A Solids, 19:93–114, 1994.
  * [81] Michel Jean.  The non-smooth contact dynamics method.  Computer Methods in Applied Mechanics and Engineering, 177(3-4):235–257, 1999.
  * [82] Sinisa Dj Mesarovic and Norman A Fleck.  Frictionless indentation of dissimilar elastic-plastic spheres.  International Journal of Solids and Structures, 37(46):7071–7091, 2000.
  * [83] Yuan Chen, Didier Imbault, and Pierre Dorémus.  Numerical simulation of cold compaction of 3d granular packings.  In Progress in Powder Metallurgy, volume 534 of Materials Science Forum, pages 301–304. Trans Tech Publications Ltd, 1 2007.
  * [84] Vincent Acary and Michel Jean.  Numerical simulation of monuments by the contact dynamics method.  In DGEMN-LNEC-JRC, editor, Monument-98, Workshop on seismic perfomance of monuments, pages 69–78, Lisbon, Portugal, November 1998. Laboratório Nacional de engenharia Civil (LNEC), Lisboa, Portugal.
  * [85] Hong-Phong Cao.  Modélisation par éléments discrets rigides et/ou déformables des milieux granulaires et des troisiémes corps solides : Influence du comportement local sur le comportement global.  PhD thesis, Mécanique, Energétique, Génie Civil, Acoustique, INSA de Lyon, Lyon, 6 2011.
  * [86] Thi-Lo Vu, Saeid Nezamabadi, and Serge Mora.  Effects of particle compressibility on structural and mechanical properties of compressed soft granular materials.  Journal of the Mechanics and Physics of Solids, 146:104201, 2021\.
  * [87] Saeid Nezamabadi, Xavier Frank, Jean Yves Delenne, Julien Averseng, and Farhang Radjai.  Parallel implicit contact algorithm for soft particle systems.  Computer Physics Communications, 237:17–25, 2019.