# Mechanics of micropillar confined thin film plasticity

Abhishek Arora Rajat Arora Amit Acharya Department of Civil & Environmental
Engineering, Carnegie Mellon University, Pittsburgh, PA 15213Siemens Corporate
Technology, Princeton, NJDepartment of Civil & Environmental Engineering, and
Center for Nonlinear Analysis, Carnegie Mellon University, Pittsburgh, PA
15213, email<EMAIL_ADDRESS>
###### Abstract

Micropillar compression experiments probing size effects in confined
plasticity of metal thin films, including the indirect imposition of
‘canonical’ simple shearing boundary conditions, show dramatically different
responses in compression and shear of the film. The Mesoscale Field
Dislocation Mechanics (MFDM) model is confronted with this set of experimental
observations and shown to be capable of modeling such behavior, without any
ad-hoc modification to the basic structure of the theory (including boundary
conditions), or the use of extra fitting parameters. This is a required
theoretical advance in the current state-of-the art of strain gradient
plasticity models. It is also shown that significantly different inhomogeneous
fields can display qualitatively similar size effect trends in overall
agreement with the experimental results. The (plastic) Swift and (elastic)
Poynting finite deformation effects are also demonstrated.

##  1 Introduction

Size-effects in confined plasticity of metal thin films sandwiched in ceramic
micropillars have been demonstrated in the work of Meng and co-workers [17,
19, 18]. Shear failure testing of the interfacial regions of CrN/Cu/Si and
CrN/Ti/Si ceramic-coating/metal-adhesion layer/substrate systems through
instrumented compression of cylindrical micropillars was reported in [24]. The
last thirty years have seen intense world-wide activity in modeling size
effects in metal plasticity at the microscale, initiated by the Strain
Gradient Plasticity (SGP) work of Fleck, Hutchinson and co-workers [11]; the
study of length-scale effects in plasticity was initiated earlier by Aifantis
[4] and co-workers. The confined plasticity results in [19], however, have not
been successfully modeled by SGP as pointed out in [19], without rather
drastic modifications to the structure of the theory and introducing extra
fitting parameters [14, 15, 16]. The aim of our work is to report on the
reasonably successful modeling of the experiments in [19] with Mesoscale Field
Dislocation Mechanics theory [5, 1, 7, 2], and to provide a mechanistic
understanding of the observed effects within the idealization of the model. We
also analyze the mechanics of local fields and other interesting results
bearing on historically important observed effects in the large deformation of
elastic-plasic materials.

Conventional plasticity models do not have a material length scale and produce
size-independent response (for homogeneous materials). In SGP theory, the
material response is assumed to depend on both the plastic strain and its
spatial gradient [13, 10, 12], and the work conjugate of the plastic strain
gradient is interpreted to be a ‘microscopic stress,’ of unspecified physical
origin. SGP theories predict a much stronger dependence on the film thickness
under nominal simple shearing conditions compared to experimental observations
[19, 17]. This prompted Kuroda and Needleman [15] to introduce an ad-hoc
modification to the constrained boundary condition specification of SGP theory
where it is assumed that a threshold exists on the magnitude of the plastic
strain gradient at boundaries. Above this threshold, the constrained boundary
condition is released and plastic straining is allowed at the boundaries - as
is clear, this threshold is not a material parameter, and it is not clear what
its validity is, and how it is to be determined, in general modes of loading.
Kuroda et al. [16] use this boundary condition in a finite deformation setting
along with extensive fitting of their model parameters to the data of [19]
(including the classical work-hardening modulus and the initial yield stress)
to produce results in accord with the experimental observations on compression
of ceramic-metal thin film sandwich micropillars. Another effort to address
this shortcoming of SGP theory is that of Dahlberg and Ortiz [9] who introduce
a fractional derivative of the effective plastic strain in the material
response (a fractional SGP theory), whose fractional exponent is intended as a
fitting parameter to recover the experimentally observed scaling of shear
stress with layer thickness. Arora and Acharya [2] used the finite deformation
implementation [7] of MFDM theory, without any special fitting beyond the use
of generic material parameters used in MFDM simulations of polycrystalline
metals, to model constrained shearing of a thin film. They recovered the
observed size-effect trend in [19] corresponding to $45^{\circ}$-oriented thin
film in the micropillar experiments.

Moving beyond the simple thin film-only geometry, in this work we use the MFDM
framework to study the compression experiments of [19] involving metal thin
films sandwiched in ceramic micropillars, within a plane-strain idealization.
As in the experiments, both $45^{\circ}$ and $90^{\circ}$ oriented thin films
are considered. The metal layer is sandwiched between two ceramic blocks above
and below it, and the whole composite is put under compression. The schematic
of pillars with the thin film in the $90^{\circ}$ and $45^{\circ}$
orientations is shown in Fig. 1. As explained in detail in [19], compression
of the micropillar with the film in the $45^{\circ}$ orientation results in a
nominal simple shearing boundary condition imposed on it; compression with the
film in the $90^{\circ}$ orientation results in compression boundary
conditions on the film, with restrained lateral movement at the top and bottom
of the boundaries. The two configurations result in dramatically different
observations of size effects, which we recover.

This paper is organised as follows: the following paragraph contains some
notational details. Section 2 recalls the MFDM governing equations, boundary,
and initial conditions. For details on the physical basis of the model and the
computational framework, the interested reader is referred to [1, 7, 6].
Section 3 shows the scaling obtained for applied nominal stress with thin film
thickness, for both $90^{\circ}$ and $45^{\circ}$ orientations of the film,
and provides a mechanistic explanation for the observations. It is also shown
that a free-standing film under simple shear produces normal stress, as
observed in metals bars under torsion [8], and the Poynting effect in non-
linear elastic solids [21]. Section 4 contains some concluding remarks.

Vectors and tensors are represented by bold face lower and upper-case letters.
The action of a second order tensor ${\mathbold A}$ on a vector ${\mathbold
b}$ is denoted by ${\mathbold A}{\mathbold b}$. The inner product of two
vectors is denoted by ${\mathbold a}\cdot{\mathbold b}$, while the inner
product of two second order tensors is denoted by ${\mathbold A}:{\mathbold
B}$. A rectangular Cartesian coordinate system is invoked for ambient space
and all (vector) tensor components are expressed with respect to the basis of
this coordinate system. $(\cdot)_{,i}$ denotes the partial derivative of the
quantity $(\cdot)$ w.r.t. the $x_{i}$ coordinate of this coordinate system.
${\mathbold e}_{i}$ denotes the unit vector in the $x_{i}$ direction.
Einstein’s summation convection is always implied, unless mentioned otherwise.
The symbols $\mathop{\rm grad}\nolimits$, $\mathop{\rm div}\nolimits$, and
$\mathop{\rm curl}\nolimits$ denote the gradient, divergence, and curl on the
current configuration. For a second order tensor ${\mathbold A}$, vectors
${\mathbold v}$, ${\mathbold a}$, and ${\mathbold c}$, a spatially constant
vector field ${\mathbold b}$, the operations of $\mathop{\rm div}\nolimits$,
$\mathop{\rm curl}\nolimits$, and cross-product of a tensor $(\times)$ with a
vector are defined as follows:

$\displaystyle\left(\mathop{\rm div}\nolimits{\mathbold
A}\right)\cdot{\mathbold b}$ $\displaystyle=\mathop{\rm
div}\nolimits\left({\mathbold A}^{T}{\mathbold
b}\right),\qquad\forall\>\>{\mathbold b}$ $\displaystyle{\mathbold
b}\cdot\left(\mathop{\rm curl}\nolimits{\mathbold A}\right){\mathbold c}$
$\displaystyle=\left[\mathop{\rm curl}\nolimits\left({\mathbold
A}^{T}{\mathbold b}\right)\right]\cdot{\mathbold
c},\qquad\forall\>\>{\mathbold b},{\mathbold c}$ $\displaystyle{\mathbold
c}\cdot\left({\mathbold A}\times{\mathbold v}\right){\mathbold a}$
$\displaystyle=\left[\left({\mathbold A}^{T}{\mathbold
c}\right)\times{\mathbold v}\right]\cdot{\mathbold
a}\qquad\forall\>\>{\mathbold a},{\mathbold c}.$

In rectangular Cartesian coordinates, these are denoted by

$\displaystyle\left(\mathop{\rm div}\nolimits{\mathbold A}\right)_{i}$
$\displaystyle=A_{ij,j},$ $\displaystyle\left(\mathop{\rm
curl}\nolimits{\mathbold A}\right)_{ri}$
$\displaystyle=\varepsilon_{ijk}A_{rk,j},$ $\displaystyle\left({\mathbold
A}\times{\mathbold v}\right)_{ri}$
$\displaystyle=\varepsilon_{ijk}A_{rj}v_{k},$

where $\varepsilon_{ijk}$ are the components of the third order alternating
tensor ${\mathbold X}$.

(a)

(b)

Figure 1: Schematic of idealized plane strain model of micropillar compression
experiments.

##  2 Theory

The governing equations for MFDM [7] are as follows:

$\displaystyle\mathring{\boldsymbol{\alpha}}\equiv(\mathop{\rm
div}\nolimits(\boldsymbol{v}))\boldsymbol{\alpha}+\dot{\boldsymbol{\alpha}}-\boldsymbol{\alpha}{\mathbold
L}^{T}=-\mathop{\rm curl}\nolimits(\boldsymbol{\alpha}\times{\mathbold
V}+{\mathbold L}^{p})$ (3a) $\displaystyle\begin{array}[]{c}{\mathbold
W}=\boldsymbol{\chi}+\mathop{\rm grad}\nolimits{\mathbold
f},\>\boldsymbol{F}^{e}=\boldsymbol{W}^{-1}\\\ \mathop{\rm
curl}\nolimits\boldsymbol{\chi}=-\boldsymbol{\alpha}\\\ \mathop{\rm
div}\nolimits\boldsymbol{\chi}=\boldsymbol{0}\end{array}\Bigg{\\}}$ (3e)
$\displaystyle\mathop{\rm div}\nolimits\left(\mathop{\rm
grad}\nolimits\dot{\boldsymbol{f}}\right)=\mathop{\rm
div}\nolimits\left(\boldsymbol{\alpha}\times\boldsymbol{V}+{\mathbold
L}^{p}-\dot{\boldsymbol{\chi}}-\boldsymbol{\chi}\boldsymbol{L}\right)$ (3f)
$\displaystyle\mathop{\rm div}\nolimits\left[{\mathbold T}({\mathbold
W})\right]=\bigg{\\{}\begin{array}[]{l}\boldsymbol{0}\quad\text{quasistatic}\\\
\rho\dot{{\mathbold v}}\quad\text{dynamic}.\end{array}$ (3i)

Here, $\boldsymbol{F}^{e}$ is the elastic distortion tensor and
$\boldsymbol{W}$ is the inverse of elastic distortion tensor,
$\boldsymbol{\chi}$ is the incompatible part of ${\mathbold W}$, ${\mathbold
f}$ is the plastic position vector, and $\mathop{\rm grad}\nolimits{\mathbold
f}$ is the compatible part of ${\mathbold W}$. $\mathbold{\alpha}$ is the
dislocation density tensor, ${\mathbold v}$ is the material velocity,
${\mathbold L}=\mathop{\rm grad}\nolimits{\mathbold v}$ is the velocity
gradient, ${\mathbold T}$ is the (symmetric) Cauchy stress tensor, while
${\mathbold L}^{p}$ is the additional meso-scale field which represents the
averaged rate of plastic straining due to all dislocations that cannot be
resolved by $\mathbold{\alpha}\times{\mathbold V}$, where both fields
represent space-time running averages. The fields in the MFDM framework are
the running space-time averages of corresponding fields of FDM theory.

To close the above set of equations, constitutive laws for
$\boldsymbol{V},\boldsymbol{L}^{p},\boldsymbol{T}$ are derived by keeping the
mechanical dissipation to be non-negative [1].

In the present work, plastic straining in accordance to
$\mathrm{\textit{J}_{2}}$ plasticity theory is augmented with the regularizing
effects of an idealization of microscopic dislocation core energy to get
${\mathbold L}^{p}$.

Initial conditions for $\mathbold{\alpha},{\mathbold f}$ and boundary
conditions for $\mathbold{\alpha},{\mathbold f},\mathbold{\chi}$ and
${\mathbold v}$ are specfied for well-set evolution.

###  2.1 Boundary Conditions

  * •

The $\boldsymbol{\alpha}$ evolution equation has a convective boundary
condition of the form $(\mathbold{\alpha}\times{\mathbold V}+{\mathbold
L}^{p})\times{\mathbold n}=\boldsymbol{\Phi}$, where $\boldsymbol{\Phi}$ is a
second order tensor valued function of time and position on the boundary
characterising the flux of dislocations at the surface satisfying the
constraint $\boldsymbol{\Phi}{\mathbold n}=\mathbf{0}$.

There are two ways in which the boundary condition is specified: (a)
Constrained: $\boldsymbol{\Phi}({\mathbold x},t)=\mathbf{0}$ at a point
${\mathbold x}$ on the boundary for all times, which ensures that there is no
outflow of dislocations at that point of the boundary, and only parallel
motion along the boundary is allowed. (b) Unconstrained: A less restrictive
boundary condition where $\hat{{\mathbold L}}^{p}\times{\mathbold n}$ is
simply evaluated at the boundary (akin to an outflow condition), along with
the specification of dislocation flux $\mathbold{\alpha}({\mathbold
V}\cdot{\mathbold n})$ on the inflow part of the boundary. Additionally, for
all calculations presented in this paper $(\mathop{\rm
curl}\nolimits\mathbold{\alpha}\times{\mathbold n})=0$ is imposed, a
particular specification of a boundary condition that arises from simple
mathematical modeling of the manifestation of dislocation core energy at the
mesoscale.

  * •

For the incompatibility equation, $\boldsymbol{\chi}\boldsymbol{n}=\mathbf{0}$
is applied on the outer boundary of the domain, which along with the system
(3e) ensures that $\mathbold{\chi}$ vanishes when $\mathbold{\alpha}$ is zero
in the entire domain.

  * •

The $\boldsymbol{f}$ evolution equation requires a Neumann boundary condition
i.e., $(\mathop{\rm grad}\nolimits\dot{{\mathbold f}}){\mathbold
n}=(\mathbold{\alpha}\times{\mathbold V}+{\mathbold
L}^{p}-\dot{\mathbold{\chi}}-\mathbold{\chi}{\mathbold L}){\mathbold n}$ on
the outer boundary of the domain.

  * •

The material velocity boundary conditions are applied based on the loading
type, which is discussed later in Section 3.

###  2.2 Initial Conditions

  * •

The initial condition $\boldsymbol{\alpha}(\boldsymbol{x},0)=\boldsymbol{0}$
is assumed for all sample sizes.

  * •

In general, the initial condition for $\boldsymbol{f}$ is obtained by solving
for $\mathbold{\chi}$ from the incompatibility equation and solving for
$\boldsymbol{f}$ from the equilibrium equation, for prescribed
$\boldsymbol{\alpha}$ on the given initial configuration. We refer to this
scheme as the elastic theory of continuously distributed of dislocations
(ECDD) solve on that configuration. For the inital conditions on
$\mathbold{\alpha}$ considered above, this step is trivial, with ${\mathbold
f}={\mathbold X}$, where ${\mathbold X}$ is the position field on the initial
configuration.

  * •

The model admits an arbitrary specification of $\dot{{\mathbold f}}$ at a
point to uniquely evolve ${\mathbold f}$ using (3f) in time, and this rate is
prescribed to vanish.

###  2.3 Constitutive Relations

Constitutive relations in MFDM are required for the stress ${\mathbold T}$,
the plastic distortion rate ${\mathbold L}^{p}$, and the dislocation velocity
${\mathbold V}$. The details of the thermodynamically consistent constitutive
formulations can be found in Sec. 3.1 of [1]. Table 1 presents the
constitutive relation for Cauchy stress and mesoscopic core energy density for
the material. Tables 2 and 3 show the constitutive relations for plastic
distortion rate and dislocation velocity, respectively. Table 4 shows the
evolution equation for material strength.

Saint-Venant-Kirchhoff Material | $\displaystyle\phi({\mathbold W})=\frac{1}{2\rho^{*}}{\mathbold E}^{e}:\mathbb{C}:{\mathbold E}^{e},\quad{\mathbold T}={\mathbold F}^{e}[\mathbb{C}:{\mathbold E}^{e}]{{\mathbold F}^{e}}^{T}$
---|---
Core energy density | $\displaystyle\Upsilon(\mathbold{\alpha})=\frac{1}{2\rho^{*}}\epsilon\mathbold{\alpha}:\mathbold{\alpha}$
Table 1: Constitutive relations for Cauchy stress and core energy density. $\mathrm{\textit{J}_{2}}$ plasticity | $\displaystyle\hat{{\mathbold L}}^{p}=\hat{\gamma}{\mathbold W}\frac{{\mathbold T}^{{}^{\prime}}}{|{\mathbold T}^{{}^{\prime}}|};\quad\hat{\gamma}=\hat{\gamma}_{0}\left(\frac{|{\mathbold T}^{{}^{\prime}}|}{\sqrt{2}g}\right)^{\frac{1}{m}}$
---|---
$\displaystyle{\mathbold L}^{p}=\hat{{\mathbold
L}}^{p}+l^{2}\hat{\gamma}\mathop{\rm curl}\nolimits\mathbold{\alpha}$
Table 2: Constitutive relations for plastic strain rate ${\mathbold L}^{p}$. $\displaystyle T^{{}^{\prime}}_{ij}=T_{ij}-\frac{T_{mm}}{3}\delta_{ij}$; | $\displaystyle a_{i}=\frac{1}{3}T_{mm}\varepsilon_{ijk}F^{e}_{jp}\alpha_{pk};$ | $\displaystyle c_{l}=\varepsilon_{ijk}T^{{}^{\prime}}_{jr}F^{e}_{rp}\alpha_{pk}$
---|---|---
$\displaystyle{\mathbold d}={\mathbold c}-\left({\mathbold c}-\frac{{\mathbold a}}{|{\mathbold a}|}\right)\frac{{\mathbold a}}{|{\mathbold a}|};$ | $\displaystyle{\mathbold V}=\zeta\frac{{\mathbold d}}{|{\mathbold d}|};$ | $\displaystyle\zeta=\left(\frac{\mu}{g}\right)^{2}\eta^{2}b\hat{\gamma}$
Table 3: Constitutive relations for dislocation velocity ${\mathbold V}$. $\displaystyle\dot{g}=h(\boldsymbol{\alpha},g)\left(|\boldsymbol{F}^{e}\boldsymbol{\alpha}\times\boldsymbol{V}|+\hat{\gamma}\right);$ | $\displaystyle h(\boldsymbol{\alpha},g)=\frac{\mu^{2}\eta^{2}b}{2(g-g_{0})}k_{0}|\boldsymbol{\alpha}|+\Theta_{0}\left(\frac{g_{s}-g}{g_{s}-g_{0}}\right)$
---|---
Table 4: Constitutive relations for material strength $g$.

The physical meanings of the material parameters in our model are: $\mu$ is
the shear modulus, $\hat{\gamma}_{0}$ is the reference strain rate, $m$ is the
material rate sensitivity, $\eta$ is a non-dimensional material constant in
the empirical Taylor relationship for macroscopic strength vs dislocation
density, $b$ is the Burgers vector magnitude of a full dislocation in the
crystalline material, $g_{0}$ is the initial strength (initial yield stress in
shear), $g_{s}$ is the saturation strength, $\Theta_{0}$ is the Stage II
hardening rate, $k_{0}$ (non-dimensional) characterizes the hardening rate due
to geometrically necessary dislocations (GNDs), and $l$ is a material length
related to the gross modeling of mesoscale effects of dislocation core energy.

All parameters in our model, except $k_{0}$ and $l$, are part of the
constitutive structure of well-accepted models of classical plasticity theory.
The parameter $k_{0}$ was introduced in [3]. The length scale $l$ simply
controls the refinement of the GND microstructure and does not play a
physically significant role in our results.

##  3 Results and Discussion

With reference to Fig. 1, all micropillar compression simulations are
performed on initial domain sizes of $5\mu m\times 10\mu m$, containing
polycrystalline Cu thin films in $45^{\circ}$ and $90^{\circ}$ orientations.
Four (initial) thin film thicknesses of $0.5\mu m$, $0.8\mu m$, $1.0\mu m$,
and $1.2\mu m$ are considered in both orientations. Simulations of thin films
with the same thicknesses and of $5\mu m$ width in free-standing
configurations are also performed, under compression and shear loading. The
nominal compression loading rate is $|\hat{e}|=0.001\,s^{-1}$ for both
micropillar and free-standing film configurations, and the simple shear
loading rate is $\hat{\Gamma}=0.001\,s^{-1}$, for the free-standing films. At
any time $t$, the nominal compressive strain is $|e|=|\hat{e}|t$, while the
nominal shear strain is $\Gamma=\hat{\Gamma}t$.

For the micropillar simulations, the interfaces between the thin film and the
ceramic blocks are assumed to be plastically unconstrained. For the free-
standing films, the top/bottom boundaries are plastically constrained, while
the left/right boundaries are plastically unconstrained.

The boundary conditions for material velocity for the micropillar simulations
are as follows:

  * •

$v_{1}=v_{2}=0$ at the bottom boundary of the domain.

  * •

$v_{2}=-|\hat{e}|L$ at the top boundary of the domain, where $L$ is the height
of the pillar in the (undeformed) initial configuration at $t=0$.

  * •

The applied traction in the horizontal direction is zero on the top boundary
of the domain.

For the free-standing thin films under compression loading, the boundary
conditions are as follows:

  * •

$v_{1}=v_{2}=0$ at the bottom boundary of the domain.

  * •

$v_{1}=0$ and $v_{2}=-|\hat{e}|L$ at the top boundary of the domain. The
lateral constraint $v_{1}=0$ on the top and bottom boundaries is imposed to
model the effects of resistance to material flow along the width of the film
due to the ceramic blocks.

For the free-standing thin films under simple shear loading, the boundary
conditions are [2]:

  * •

$v_{1}=v_{2}=0$ at the bottom boundary of the domain.

  * •

$v_{1}=\hat{\Gamma}L$ and $v_{2}=0$ at the top boundary of the domain.

  * •

$v_{1}=\hat{\Gamma}y$ and $v_{2}=0$ at the left and right boundaries of the
domain, where $y$ is the difference between the $x_{2}$ coordinate of any
point on the boundary with respect to the bottom boundary.

We also define the yield strength in compression $(\sigma_{0})$ in terms of
the yield strength in shear $(g_{0})$ for different cases:
$\sigma_{0}=1.14\sqrt{3}g_{0}$ for the free standing film under compression
and for the micropillar sandwich with $90^{\circ}$ orientation, while
$\sigma_{0}=2g_{0}$ for the micropillar sandwich with $45^{\circ}$
orientation.

Parameter | $\hat{\gamma}_{0}$ | $m$ | $\eta$ | $b$ | $g_{0}$ | $g_{s}$ | $\Theta_{0}$ | $k_{0}$ | $l$ | $E$ | $\nu$
---|---|---|---|---|---|---|---|---|---|---|---
($s^{-1}$) |  |  | (Å) | (MPa) | (MPa) | (MPa) |  | ($\mu m$) | (GPa) |
Value | 0.001 | 0.03 | $\frac{1}{3}$ | 4.05 | 17.3 | 161 | 392.5 | 20 | $\sqrt{3}\times 0.1$ | 62.78 | 0.3647
Table 5: Material parameters for metal. Parameter |  $E$ (GPa) | $\nu$
---|---|---
Value | 110.0 | 0.20
Table 6: Material parameters for ceramic.

The parameter values used for all the simulations here are given in Table 5
(identical to [2] for the metal) and Table 6 for the ceramic. These are the
typical material constants for copper and silicate ceramics. Here, $E$ and
$\nu$ are the Young’s modulus and the Poisson’s ratio of the materials.

###  3.1 Size effects

We begin with the size effect results for the free-standing thin films under
simple shear and compression loadings. These simulations are an idealization
of the compression experiments on sandwiched thin films in the $45^{\circ}$
and $90^{\circ}$ orientations within micropillars. These simplified cases
already provide the basic explanation for the different behavior of the
micropillars with thin films in the two orientations. As mentioned earlier,
the constrained boundary conditions reflect the constraints of the ceramic
blocks on the films. Without these constraints the deformation is essentially
homogeneous and no size effect is observed. With the no plastic flow
constraint in effect, but no lateral constraint on material velocity at the
top and bottom boundaries in the case of compression, a weak size effect,
similar in both shear (Fig. 2) and compression is predicted; the shear results
were also obtained in [2, Sec. 4.1], in accord with the experimental trends
for the film in the $45^{\circ}$ orientation. With the lateral b.c. constraint
on material velocity in compression, a strong size effect in compression,
shown in Fig. 3, is observed. As already observed in [19], this is essentially
due to the inhomogeneous lateral material deformation, from the top and bottom
boundaries to the center of the film, induced by the material velocity
boundary conditions on the film; such inhomogeneity is absent in simple shear
loading with the no-plastic flow constraint. In both cases, we observe that
the metal film with the smallest thickness hardens the most, due to higher
gradients in ${\mathbold L}^{p}$ for smaller domain sizes (by scaling
arguments), which are produced due to the constrained plastic flow boundary
condition on the top/bottom boundary.

With this basic understanding of the gross behavior, we then model the
compression of the various micropillar configurations. Interestingly, we again
recover the size-effect trends observed in experiment, but we demonstrate
significant differences in local mechanical fields, in comparison to the free-
standing films, that we describe in Sec. 3.2.

Fig. 4 shows the stress-strain curves and the fitted power law relationship
between the stress and film thickness i.e.,
$(\sigma-\sigma_{0})/\sigma_{0}=a\>H^{m}$. Here, $|m|$ approaches 1.0 at $6\%$
overall strain, and the same value for the exponent is reported in the
experimental work of Meng and co-workers [19] on Cu-CrN as-deposited
micropillars. As shown in Fig. 5, for the micropillar sandwich with the thin
film in the $90^{\circ}$ orientation, the film keeps bulging out with increase
in the applied compressive strain, as similarly observed in the post test
scanning electron microscope images of Cu interlayers in [24]. Fig. 5 also
shows the norm of the Logarithmic strain tensor, $|ln\>(\boldsymbol{V})|$ at
$6\%$ overall compression strain for different film thickness, where
$\boldsymbol{V}$ is the left stretch tensor of the Polar Decomposition of
${\mathbold F}$ from the initial configuration at $t=0$.

(a)  (b)

Figure 2: (a) Stress-strain curve in simple shear of the free-standing film
for different film thickness $(H)$. (b) Power law fitted for the stress vs
film thickness ($(\tau-g_{0})/g_{0}=a\>H^{m}$) at $10\%-40\%$ strain. The
first and second term in the legend denote the strain and magnitude of power
law exponent $(|m|)$, respectively.

(a)  (b)

Figure 3: (a) Stress-strain curve in compression of the free-standing film for
different film thickness $(H)$. (b) Power law fitted for the stress vs film
thickness ($(\sigma-\sigma_{0})/\sigma_{0}=a\>H^{m}$) at $5\%-15\%$ strain.
The first and second term in the legend denote the strain and magnitude of
power law exponent $(|m|)$, respectively.

(a)  (b)

Figure 4: (a) Stress-strain curve in compression of the micropillar sandwich
with $90^{\circ}$ film orientation and for different film thickness $(H)$ .
(b) Power law fitted for the stress vs film thickness
($(\sigma-\sigma_{0})/\sigma_{0}=a\>H^{m}$) at $1\%-6\%$ strain. The first and
second term in the legend denote the strain and magnitude of power law
exponent $(|m|)$, respectively.

(a)

(b)

(c)

Figure 5: Norm of Logarithmic strain tensor $(|ln(\boldsymbol{V})|)$ at 6%
compression strain $(|e|)$ for the micropillar sandwich with $90^{\circ}$ film
orientation for different film thickness $(H)$.

(a)  (b)

Figure 6: (a) Stress-strain response in compression of the micropillar
sandwich with $45^{\circ}$ orientation and for different film thickness $(H)$.
(b) Power law fitted for the stress vs film thickness ($(\tau-
g_{0})/g_{0}=a\>H^{m}$) at $1\%-3.5\%$ strain. The first and second term in
the legend denote the strain and magnitude of power law exponent $(|m|)$,
respectively.

Fig. 6(a) shows the stress-strain response for micropillars with the metal
thin film in the $45^{\circ}$ orientation, and Fig. 6(b) shows the fitted
power-law relationship between the ‘applied shear stress’ (for the thin film
in the given orientation, given by $\tau=\frac{1}{2}\sigma$, where $\sigma$ is
the applied nominal compressive stress) and film thickness i.e., $(\tau-
g_{0})/g_{0}=a\>H^{m}$. The magnitude of the power-law exponent $(|m|)$
decreases from 0.728 to 0.347 as $|e|$ increases from $1.0\%$ to $3.0\%$. The
$|m|$ obtained at $3.5\%$ strain is 0.395, however, if only samples with
$H>=0.8\>\mu m$ are considered, then $|m|$ obtained will be close to zero. Our
model for the plastic straining due to statistical dislocations (i.e.
${\mathbold L}^{p}$) is phenomenological and it is expected that as the
overall length scale over which plastic flow occurs decreases, the
contribution of ${\mathbold L}^{p}$ in the total plastic strain rate for the
model should decrease (physically, there are fewer and fewer sources, but the
phenomenological $J_{2}$ plasticity model assumes that there is an abundant
supply of sources, and all that is required is stress to mobilize them). This
is not accounted for in the current model and at small scales, e.g.
$H=0.5\>\mu m$, there is excessive hardening due to higher gradients in
${\mathbold L}^{p}$.

An eventual reduction in the hardening modulus in the shear stress-strain
response in Fig. 6(a) is observed, approaching a ‘flat’ response on average,
but not the shear stress plateau as experimentally observed by [19]. Our model
does not have a failure model for the interfacial regions which would be
capable of demonstrating such response under the intense local shearing
observed in the experiments and in the simulations. Kuroda et al. [16] are
able to demonstrate the observed plateau for the film in the $45^{\circ}$
orientation in the micropillar; a small value of the hardening modulus is
employed in their calculations, a material description that would be hard-
pressed to reproduce the stress-strain behaviour of polycrystalline copper for
macroscopic pillar and film sizes.

###  3.2 Local fields in the free-standing film vs. the micropillar sandwich

In this section, we show the significant differences in the local mechanical
response of the free-standing films and the micropillar-film sandwich (which,
nevertheless, produce the same qualitative size effects, as already shown).

When comparing the local fields for the micropillar sandwich with the film in
the $45^{\circ}$ orientation and the free-standing film, the film for the
micropillar case is rotated by $45^{\circ}$ in an anti-clockwise sense, for
the ease of visualization.

In the following, we define $|e_{m}|$ (the subscript $m$ stands for metal) as
the nominal compressive strain in the metal film corresponding to a given
nominal compressive strain ($|e|$) for the entire domain - for the free-
standing film, $|e_{m}|=|e|$, whereas for the micropillar sandwich, it is the
magnitude of the engineering compressive strain calculated from the applied
boundary condition on the top of the pillar and the initial length of the
pillar.

We also define $\Gamma_{m}$ as the nominal shear strain in the metal film for
micropillar sandwich as well as the free-standing film.

####  3.2.1 Micropillar with thin film in the $90^{\circ}$ orientation

Fig. 7 compares the $\rho_{g}:=|\boldsymbol{\alpha}|/b$ (GND) field plots at
$|e_{m}|=15\%$ compressive strain in the metal film for the free-standing film
and the micropillar sandwich. Excess dislocation boundary layers are absent in
the micropillar film configuration. They arise in the free-standing films due
to no-flow boundary conditions. The $\rho_{g}$ patterns observed in the
micropillar film are driven by the deformation inhomogeneity induced by the
the resistance to horizontal material flow at the boundaries of the film with
the ceramic blocks.

(a)

(b)

(c)

Figure 7: $\rho_{g}=|\boldsymbol{\alpha}|/b\>(m^{-2})$ at $|e_{m}|=15\%$ in
the metal film for the free-standing film and the micropillar sandwich with
$90^{\circ}$ orientation. The corresponding $|e|$ for the micropillar is
$1.56\%$.

(a)

(b)

(c)

Figure 8: von Mises stress at $|e_{m}|=15\%$ in the metal film for the free-
standing film and the micropillar sandwich with $90^{\circ}$ orientation. The
corresponding $|e|$ for the micropillar is $1.56\%$.

(a)

(b)

(c)

Figure 9: Hydrostatic stress at $|e_{m}|=15\%$ in the metal film for the free-
standing film and the micropillar sandwich with $90^{\circ}$ orientation. The
corresponding $|e|$ for the micropillar is $1.56\%$.

Fig. 8 shows the von Mises stress and Fig. 9 the hydrostatic stress fields
(both normalised with yield stress) at $|e_{m}|=15\%$ in the metal film for
the free-standing film and the micropillar sandwich. The magnitude of
$\sigma_{v}$ for the free-standing film in Fig. 8 is higher at the top/bottom
boundary of the thin film. The hydrostatic stress is much higher in the case
of free-standing film as compared to the micropillar case. The constraint to
lateral motion at the top and bottom boundaries of the film is softer for the
film in the micropillar sandwich, which results in the generation of less
hydrostatic stress.

####  3.2.2 Micropillar with thin film in the $45^{\circ}$ orientation

(a)

(b)

Figure 10: Nominal shear strain ($\Gamma_{m}$) in the film at $|e|=3.5\%$ for
micropillar with $45^{\circ}$ orientation, for different film thickness $(H)$.
The green solid line denotes an undeformed line, while the cyan solid line
denotes the corresponding deformed line.

Fig. 10 shows the deformed mesh and the nominal shear strain ($\Gamma_{m}$) in
the metal film at $|e|=3.5\%$ compressive strain, for different film
thickness. Under increased compression, the film material near the lateral
boundaries of the micropillar rotate excessively due to the shearing and
(damage-free) constraint of the ceramic blocks. The calculations are stopped
when the film material is close to penetrating the ceramic blocks, as shown in
Fig. 11. Although the magnitude of $|e|$ is not that high, $\Gamma_{m}$ in the
metal film is very large. For instance, at $|e|=3.5\%$, $\Gamma_{m}$ in the
metal film is $102\%$ for $H=0.5\>\mu m$, while it is $39\%$ for $H=1.2\>\mu
m$.

Figure 11: The deformed mesh for the micropillar with $45^{\circ}$ orientation
for $H=0.5\>\mu m$ at $|e|=3.8\%$, when the film material is close to
penetrating into ceramic blocks.

(a)

(b)

(c)

Figure 12: Norm of Logarithmic strain tensor $(|ln(\boldsymbol{V})|)$ at 3.5%
compression strain in the micropillar sandwich with $45^{\circ}$ orientation,
for different film thickness $(H)$. The solid red line shows the boundary of
the undeformed metal thin film.

(a)

(b)

(c)

Figure 13: $\rho_{g}=|\boldsymbol{\alpha}|/b\>(m^{-2})$ at $\Gamma_{m}=30\%$
in the metal film for the micropillar sandwich with $45^{\circ}$ orientation
and for the free-standing film. The corresponding $|e|$ for the micropillar is
$1.15\%$, and the solid red line in (a) shows the boundary of the undeformed
metal thin film.

(a)

(b)

(c)

Figure 14: von Mises stress at $\Gamma_{m}=30\%$ in the metal film for the
micropillar sandwich with $45^{\circ}$ orientation and for the free-standing
film. The corresponding $|e|$ for the micropillar is $1.15\%$, and the solid
red line in (a) shows the boundary of the undeformed metal thin film.

(a)

(b)

(c)

Figure 15: Hydrostatic stress at $\Gamma_{m}=30\%$ in the metal film for the
micropillar sandwich with $45^{\circ}$ orientation and for the free-standing
film. The corresponding $|e|$ for the micropillar is $1.15\%$, and the solid
red line in (a) shows the boundary of the undeformed metal thin film.

Fig. 12 shows field plots of the norm of the Logarithmic strain tensor
$\left(|ln(\boldsymbol{V})|\right)$ in the film within the micropillar, at
$|e|=3.5\%$ nominal compressive strain. As expected, for the smaller
thickness, the magnitude of $|ln(\boldsymbol{V})|$ is higher as the nominal
shear strain in the metal film is higher for smaller thickness (for equal
lateral displacement).

Fig. 13 shows the $\rho_{g}$ field plots at $\Gamma_{m}=30\%$ for the free-
standing film and the micropillar sandwich with $45^{\circ}$ orientation. As
for the $90^{\circ}$ orientation case, there are no boundary layers in the
film within the micropillar here, while they are present for the free-standing
film.

Fig. 14 and 15 respectively show the von Mises stress $(\sigma_{v})$ and the
hydrostatic stress $(\sigma_{h})$ field plots for the free-standing and
micropillar sandwich configurations.

###  3.3 Hydrostatic and normal stress in shear

In small deformation elasticity or metal plasticity, shearing induced normal
stresses are rarely observed. However, such couplings are observed in finite
elastic [21] and plastic [8] deformations. We explore such possibilities in
our simulations of simple shearing of free-standing films and compression of
micropillar sandwich with $45^{\circ}$ film orientation. To do so, we define
the ‘driving’ stress ($\sigma_{2}^{ap}$) in both configurations as the average
reaction stress in the vertical direction on the top boundary of the domain.
The averages of the various components of the Cauchy stress tensor over the
metal film are also defined as ($\sigma_{1}^{avg}$, $\sigma_{2}^{avg}$,
$\sigma_{3}^{avg}$, $\tau^{avg}$), denoting the averages of the ($T_{11}$,
$T_{22}$, $T_{33}$, $T_{12}$) components, respectively. The averaged generated
hydrostatic stress in the film is denoted by $\sigma_{h}^{avg}$.

The values of the driving stress and the averages of generated stress
components for the micropillar sandwich and the free-standing film at
$\Gamma_{m}=30\%$ are reported in Table 7. The state of loading is multi-axial
for both the micropillar sandwich and the free-standing film. However, the
‘average stress state’ in the case of the micropillar sandwich is more
hydrostatic (confined) in nature, while it is predominantly shear driven for
the free-standing film, as also shown in Fig. 15. The magnitude of hydrostatic
stress in the free-standing film is much smaller than in the film in the
micropillar, primarily due to the applied compressive loading on the pillar;
nevertheless, the free-standing film generates compressive normal reactions at
the boundary, for solely applied shear velocity boundary conditions (along
with a constraint to motion of the boundary in the normal direction, of
course).

Specimen | $\sigma_{2}^{ap}/\sigma_{0}$ | $\sigma_{1}^{avg}/\sigma_{0}$ | $\sigma_{2}^{avg}/\sigma_{0}$ | $\sigma_{3}^{avg}/\sigma_{0}$ | $\tau^{avg}/\sigma_{0}$ | $\sigma_{h}^{avg}/\sigma_{0}$
---|---|---|---|---|---|---
Micropillar sandwich | -6.086 | -0.122 | -6.241 | -2.283 | 0.109 | -2.882
Free-standing film | -0.048 | -0.01 | -0.035 | -0.026 | 2.912 | -0.024
Table 7: Driving stress and averaged stress components in the micropillar
sandwich and the free-standing film at $\Gamma_{m}=30\%$.

Despite this significant difference, the averaged hydrostatic stress
normalized by the driving normal stress is of the same order of magnitude for
both cases:

$\displaystyle\sigma_{h}^{avg}/\sigma_{2}^{ap}$
$\displaystyle=0.474,\quad\text{(Micropillar Sandwich)}$ (4a)
$\displaystyle\sigma_{h}^{avg}/\sigma_{2}^{ap}$
$\displaystyle=0.497\quad\text{(Free-Standing Film)}$ (4b)

(considering a normalization by applied stress component normal to the film in
the micropillar increases the result by a factor of 2). This suggests a
unifying ‘collapse’ of data for understanding normal stress effects in nominal
simple shearing, produced by significantly different applied loading
conditions.

Figure 16: Normal and transverse reaction stress on the top boundary of the
free-standing film subjected to simple shearing, for $H=0.5\>\mu m$.

The normal stress observed here arises as a combination of the Swift effect
[8] and the Poynting effect [21]. One way to think about how these normal
stresses are generated in the free-standing film is as follows: the $T_{22}$
component of the Cauchy stress tensor is non-zero, a nonlinear elastic
Poynting effect to begin with. This further generates $T^{{}^{\prime}}_{22}$
and $T^{{}^{\prime}}_{11}$ components of the deviatoric stress tensor, and
these components of stress generate plastic straining in the $L^{p}_{11}$ and
$L^{p}_{22}$ components. However, due to the deformation constraints imposed
by the simple shear boundary conditions in the $x_{1}$-direction on the
lateral sides of the specimen, and the $x_{2}$-direction on the top and bottom
of the domain, further $T_{11}$ and $T_{22}$ components of stress are
generated. All of this combined generates a non-zero hydrostatic stress field
and the idealized Swift effect shown in Fig. 16, in simple shearing of the
free-standing films.

###  3.4 Effect of b.c.s on plastic flow in the micropillar sandwich

Figure 17: Stress-strain curve for the micropillar sandwich with $45^{\circ}$
orientation for $H=0.8\>\mu m$ case and with C and UC metal-ceramic interface.

(a)

(b)

(c)

Figure 18: $\rho_{g}=|\boldsymbol{\alpha}|/b\>(m^{-2})$ at $|e|=3.5\%$ for
$H=0.8\mu m$ with C and UC metal-ceramic interface. The solid red line shows
the boundary of the undeformed metal thin film.

Fig. 17 shows the stress-strain curve and Fig. 18 shows the $(\rho_{g})$ field
plots at $|e|=3.5\%$, for $H=0.8\>\mu m$, and for both the plastically
constrained (C) and unconstrained (UC) metal-ceramic interfaces. The stress-
strain response for the constrained case is marginally harder as compared to
the unconstrained case. For the constrained case, there is a thin boundary
layer in $\rho_{g}$ field at the metal-ceramic interface. The effect of the
interface condition is not prominent here due to strong heterogeneity in
plastic flow across the interface, and hence, having an elastic-plastic
interface with plastic flow unconstrained has more or less the same effect as
a plastically constrained interface. This is borne out also in the simple
shear loading of free-standing films, where the effect of
constrained/unconstrained boundary condition is much more prominent in stress-
strain response, as shown by [2] (refer to Fig. 6(b) in their paper), due to
the absence of constraining elastic material. Similar effects were observed in
the study of mechanical response of multi-crystalline thin films in [20]. It
was observed there that the effect of constrained/unconstrained grain boundary
conditions on the stress-strain response decreased on increasing the
misorientation between adjacent grains, as the plastic flow through a grain
boundary decreases with increase in the misorientation.

##  4 Conclusion

Experimentally observed size effects in confined metal thin films, in
different orientations within micropillars subjected to compression, have been
successfully reproduced in our work, using a computational framework based on
MFDM developed in [7]. This advances the current state-of-the-art of SGP
models addressing this set of experiments, as no extraneous assumptions about
the nature of required boundary conditions involving, in addition, extra
fitting parameters of dubious physical origin, are introduced. The results
obtained in this paper further strengthens the case, beyond [22, 23, 20, 1, 2,
7], for MFDM as an appropriate model for dislocation mediated mesoscale
plasticity. A shortcoming that needs to be addressed is to improve our model
for plastic straining due to statistical dislocations, ${\mathbold L}^{p}$, so
that its effect diminishes at smaller length scales.

## Acknowledgments

This work was supported by the grant NSF OIA-DMR $\\#2021019$.

## References

  * [1] Rajat Arora and Amit Acharya.  Dislocation pattern formation in finite deformation crystal plasticity.  International Journal of Solids and Structures, 184:114–135, 2020\.
  * [2] Rajat Arora and Amit Acharya.  A unification of finite deformation J2 Von-Mises plasticity and quantitative dislocation mechanics.  Journal of the Mechanics and Physics of Solids, 143:104050, 2020\.
  * AB [00] A. Acharya and A. J. Beaudoin.  Grain-size effect in viscoplastic polycrystals at moderate strains.  Journal of the Mechanics and Physics of Solids, 48(10):2213–2230, 2000.
  * Aif [87] E. C. Aifantis.  The physics of plastic deformation.  International journal of plasticity, 3(3):211–247, 1987.
  * AR [06] Amit Acharya and Anish Roy.  Size effects and idealized dislocation microstructure at small scales: predictions of a phenomenological model of mesoscopic field dislocation mechanics: Part I.  Journal of the Mechanics and Physics of Solids, 54(8):1687–1710, 2006.
  * Aro [19] Rajat Arora.  Computational Approximation of Mesoscale Field Dislocation Mechanics at Finite Deformation.  PhD thesis, Carnegie Mellon University, 2019.
  * AZA [20] Rajat Arora, Xiaohan Zhang, and Amit Acharya.  Finite element approximation of finite deformation dislocation mechanics.  Computer Methods in Applied Mechanics and Engineering, 367:113076, 2020.
  * Bil [77] E. W. Billington.  Non-linear mechanical response of various metals. III. Swift effect considered in relation to the stress-strain behaviour in simple compression, tension and torsion.  Journal of Physics D: Applied Physics, 10(4):553, 1977.
  * DO [19] C. F. O. Dahlberg and M. Ortiz.  Fractional strain-gradient plasticity.  European Journal of Mechanics-A/Solids, 75:348–354, 2019.
  * FH [01] N. A. Fleck and J. W. Hutchinson.  A reformulation of strain gradient plasticity.  Journal of the Mechanics and Physics of Solids, 49(10):2245–2271, 2001.
  * FMAH [94] N. A. Fleck, G. M. Muller, M. F. Ashby, and J. W. Hutchinson.  Strain gradient plasticity: theory and experiment.  Acta Metallurgica et materialia, 42(2):475–487, 1994.
  * Gud [04] P. Gudmundson.  A unified treatment of strain gradient plasticity.  Journal of the Mechanics and Physics of Solids, 52(6):1379–1406, 2004.
  * Gur [00] M. E. Gurtin.  On the plasticity of single crystals: free energy, microforces, plastic-strain gradients.  Journal of the Mechanics and Physics of Solids, 48(5):989–1036, 2000.
  * [14] M. Kuroda and A. Needleman.  Nonuniform and localized deformation in single crystals under dynamic tensile loading.  Journal of the Mechanics and Physics of Solids, 125:347–359, 2019\.
  * [15] M. Kuroda and A. Needleman.  A simple model for size effects in constrained shear.  Extreme Mechanics Letters, 33:100581, 2019.
  * KTN [21] M. Kuroda, V. Tvergaard, and A. Needleman.  Constraint and size effects in confined layer plasticity.  Journal of the Mechanics and Physics of Solids, 149:104328, 2021\.
  * MCM [14] Y. Mu, K. Chen, and W. J. Meng.  Thickness dependence of flow stress of Cu thin films in confined shear plastic flow.  MRS Communications, 4(3):129–133, 2014.
  * MHM [14] Y. Mu, J. W. Hutchinson, and W. J. Meng.  Micro-pillar measurements of plasticity in confined Cu thin films.  Extreme Mechanics Letters, 1:62–69, 2014.
  * MZHM [16] Y. Mu, X. Zhang, J. W. Hutchinson, and W. J. Meng.  Dependence of confined plastic flow of polycrystalline Cu thin films on microstructure.  MRS Communications, 6(3):289–294, 2016.
  * PDA [11] Saurabh Puri, Amit Das, and Amit Acharya.  Mechanical response of multicrystalline thin films in mesoscale field dislocation mechanics.  Journal of the Mechanics and Physics of Solids, 59(11):2400–2417, 2011.
  * Poy [09] J. H. Poynting.  On pressure perpendicular to the shear planes in finite pure shears, and on the lengthening of loaded wires when twisted.  Proceedings of the Royal Society of London. Series A, 82(557):546–559, 1909.
  * RA [06] Anish Roy and Amit Acharya.  Size effects and idealized dislocation microstructure at small scales: predictions of a phenomenological model of mesoscopic field dislocation mechanics: Part ii.  Journal of the Mechanics and Physics of Solids, 54(8):1711–1743, 2006.
  * ZAWB [15] Xiaohan Zhang, Amit Acharya, Noel J Walkington, and Jacobo Bielak.  A single theory for some quasi-static, supersonic, atomic, and tectonic scale applications of dislocations.  Journal of the Mechanics and Physics of Solids, 84:145–195, 2015\.
  * ZZM+ [17] X. Zhang, B. Zhang, Y. Mu, S. Shao, C. D. Wick, B. R. Ramachandran, and W. J. Meng.  Mechanical failure of metal/ceramic interfacial regions under shear loading.  Acta Materialia, 138:224–236, 2017.


[20] Pommaret, J.-F.: Partial Differential Equations and Group Theory, Kluwer
(1994).
http://dx.doi.org/10.1007/978-94-017-2539-2
[21] Pommaret, J.-F.: Partial Differential Control Theory, Kluwer, Dordrecht
(2001) (Zbl 1079.93001).
[22] Pommaret, J.-F.: Algebraic Analysis of Control Systems Defined by Partial
Differential Equations, in ”Advanced Topics in Control Systems Theory”,
Springer, Lecture Notes in Control and Information Sciences 311 (2005) Chapter
5, pp. 155-223.
[23] Pommaret, J.-F.: Parametrization of Cosserat Equations, Acta Mechanica,
215 (2010) 43-55.
http://dx.doi.org/10.1007/s00707-010-0292-y
[24] Pommaret, J.-F.: Spencer Operator and Applications: From Continuum
Mechanics to Mathematical Physics, in ”Continuum Mechanics-Progress in
Fundamentals and Engineering Applications”, Dr. Yong Gan (Ed.), ISBN:
978-953-51-0447–6, InTech (2012) Available from:
http://dx.doi.org/10.5772/35607
[25] Pommaret, J.-F.: The Mathematical Foundations of General Relativity
Revisited, Journal of Modern Physics, 4 (2013) 223-239.
https://dx.doi.org/10.4236/jmp.2013.48A022
[26] Pommaret, J.-F.: Macaulay Inverse Systems and Cartan-Kähler Theorem,
https://arxiv.org/abs/1411.7070
[27] Pommaret, J.-F.: Relative Parametrization of Linear Multidimensional
Systems, Multidim. Syst. Sign. Process., 26 (2015) 405-437.
DOI 10.1007/s11045-013-0265-0
[28] Pommaret, J.-F.: Airy, Beltrami, Maxwell, Einstein and Lanczos Potentials
Revisited, Journal of Modern Physics, 7 (2016) 699-728.
https://dx.doi.org/10.4236/jmp.2016.77068
[29] Pommaret, J.-F.: Deformation Theory of Algebraic and Geometric
Structures, Lambert Academic Publisher (LAP), Saarbrucken, Germany (2016). A
short summary can be found in ”Topics in Invariant Theory ”, Séminaire P.
Dubreil/M.-P. Malliavin, Springer Lecture Notes in Mathematics, 1478 (1990)
244-254.
https://arxiv.org/abs/1207.1964
[30] Pommaret, J.-F.: New Mathematical Methods for Physics, Mathematical
Physics Books, Nova Science Publishers, New York (2018) 150 pp.
[31] Pommaret, J.-F.: The Mathematical Foundations of Elasticity and
Electromagnetism Revisited, Journal of Modern Physics, 10 (2019) 1566-1595.
https://doi.org/10.4236/jmp.2019.1013104 (https://arxiv.org/abs/1802.02430 )
[32] Pommaret, J.-F.: Generating Compatibility Conditions and General
Relativity, J. of Modern Physics, 10, 3 (2019) 371-401.
https://doi.org/10.4236/jmp.2019.103025
[33] Pommaret, J.-F.: Differential Homological Algebra and General Relativity,
J. of Modern Physics, 10 (2019) 1454-1486.
https://doi.org/10.4236/jmp.2019.1012097
[34] Pommaret, J.-F.: The Conformal Group Revisited,
https://arxiv.org/abs/2006.03449
[35] Pommaret, J.-F.: A Mathematical Comparison of the Schwarzschild and Kerr
Metrics, Journal of modern Physics, 11, (2020) 1672-1710.
https://dx.doi.org/10.4236/jmp.2020.1110104
[36] Pommaret, J.-F.: Minimum Parametrization of the Cauchy Stress Operator,
Journal of Modern Physics, 12 (2021) 453-482.
https://doi.org/10.4236/jmp.2021.124032
[37] Pommaret, J.-F.: Homological Solution of the Lanczos Problems in
Arbitrary Dimension, Journal of Modern Physics, 12 (2021) 829-858.
https://doi.org/10.4236/jmp.2021.126053
[38]Quadrat, A.: Une Introduction à l’Analyse Algébrique Constructive et à ses
Applications, INRIA Research Report 7354, AT-SOP Project, july 2010. Les Cours
du CIRM, 1 no. 2: Journées Nationales de Calcul Formel (2010), p281-471
(doi:10.5802/ccirm.11).
[39] Rotman, J.J.: An Introduction to Homological Algebra, Pure and Applied
Mathematics, Academic Press (1979).
[40] Schneiders, J.-P.: An Introduction to D-Modules, Bull. Soc. Roy. Sci.
Liège, 63, 223-295 (1994).
[41] Shankar, S.: Controllability and Vector Potential,
https://arxiv.org/abs/1911.01238v2
[42] Shankar, S.: The Canonical Controller for Distributed Systems,
Multidimensional Systems and Signals Processing, 32 (2021) 303-311.
https://arxiv.org/10.1007/s11045-020-00740-1
[43] Spencer, D.C.: Overdetermined Systems of Partial Differential Equations,
Bull. Am. Math. Soc., 75 (1965) 1-114.
[44] Zerz, E.: Topics in Multidimensional Linear Systems Theory, Lecture Notes
in Control and Information Sciences, Springer, LNCIS 256 (2000).
[45] Zerz, E.: An Algebraic Analysis Approach to Linear Time-Varying Systems,
IMA Journal of Mathematical Control and Information, 23 (2005) 113-126.
https://doi.org/10.3182/20050703-6-CZ-1902.00200


# Deep Diffusion Models for Seismic Processing

Ricard Durall  Fraunhofer ITWM  Fraunhofer Center Machine Learning  Ammar
Ghanim  Fraunhofer ITWM  Fraunhofer Center Machine Learning  Mario Fernandez
Fraunhofer ITWM  Fraunhofer Center Machine Learning  École Normale Supérieure
Norman Ettrich  Fraunhofer ITWM  Fraunhofer Center Machine Learning  Janis
Keuper  Fraunhofer ITWM  Fraunhofer Center Machine Learning  IMLA, Offenburg
University

###### Abstract

Seismic data processing involves techniques to deal with undesired effects
that occur during acquisition and pre-processing. These effects mainly
comprise coherent artefacts such as multiples, non-coherent signals such as
electrical noise, and loss of signal information at the receivers that leads
to incomplete traces. In the past years, there has been a remarkable increase
of machine-learning-based solutions that have addressed the aforementioned
issues. In particular, deep-learning practitioners have usually relied on
heavily fine-tuned, customized discriminative algorithms. Although, these
methods can provide solid results, they seem to lack semantic understanding of
the provided data. Motivated by this limitation, in this work, we employ a
generative solution, as it can explicitly model complex data distributions and
hence, yield to a better decision-making process. In particular, we introduce
diffusion models for three seismic applications: demultiple, denoising and
interpolation. To that end, we run experiments on synthetic and on real data,
and we compare the diffusion performance with standardized algorithms. We
believe that our pioneer study not only demonstrates the capability of
diffusion models, but also opens the door to future research to integrate
generative models in seismic workflows.

##  1 Introduction

Deep generative learning has become an important research area in the machine
learning community, being more relevant in many applications. Namely, they are
widely used for image synthesis and various image-processing tasks such as
editing, interpolation, colourization, denoising, and super-resolution.
Recently, diffusion probabilistic models [1, 2] have emerged as a novel,
powerful class of generative learning methods. In a short period of time,
these models have achieved surprisingly high performance [3, 4, 5, 6], and
have even surpassed state-of-the-art algorithms like generative adversarial
networks [7] (GANs) and variational autoencoders [8] (VAEs).

At the same time, the geophysics community has been actively adopting deep-
learning techniques to boost and automate numerous seismic interpretation
tasks including fault picking [9, 10], salt delineation [11, 12], well-to-
seismic tie [13, 14], horizon tracking [15, 16], multiple removal [17, 18],
etc. Nonetheless, to the best of our knowledge, there has not been yet any
work exploring the application of diffusion models to seismic data and thus,
studying their potential advantages to already established deep-learning
approaches in this domain. Driven by this motivation, in this work, we study
the applicability of diffusion models for seismic processing.

Seismic imaging is essential to discover and characterize economically
worthwhile geological reservoirs, such as hydrocarbons accumulations, and to
manage the extraction of the resources stored in them. Unfortunately, recorded
seismic signals at the surface are inevitably contaminated by coherent and
incoherent noise of various nature. The process of removing the noise, while
retaining the primary signal, is called seismic processing. In this paper, we
focus on three relevant, well-known seismic processing tasks: demultiple,
denoising and interpolation. Demultiple and denoising are both removing
unwanted signals from the seismic section; the first gets rid of coherent
noise caused by reverberations of waves between strong reflectors, whereas the
latter removes incoherent noise of miscellaneous causes. The goal of
interpolation is to fill-in gaps in the image caused by limitations during
acquisition. Although at the first glance the nature of these problems might
look different or unrelated, it is possible to formulate a common framework,
in which they can be solved. This is feasible, due to the fact that the
diffusion models, like most of generative models, learn the density
distribution of the input data. In other words, unlike discriminative
approaches which draw boundaries in the data space, the generative approaches
model how data is placed throughout the space [19]. As a result, they are
powerful algorithms that can be independently applied to a large diversity of
problems.

##  2 Background

Generative models for modelling estimate the marginal distribution, denoted as
$p(x)$, over observable variables $x$, e.g., images. In the literature, we can
find different formulations that tackle this problem such as autoregressive
generative models, latent variable models, flow-based models, and energy-based
models.

###  2.1 Latent Variable Models

The main idea of this type of models is to utilize latent variables $z$ to
formulate the joint distribution $p(x,z)$, which describes the marginal
distribution as a function of learnable parameters ${\theta}$ (likelihood).
Mathematically, it can be written as:

$\begin{split}&z\sim p_{\theta}(z)\\\ &x\sim p_{\theta}(x|z)\\\
p_{\theta}(x)=\int_{z}p_{\theta}(x&,z)=\int_{z}p_{\theta}(x|z)p_{\theta}(z).\end{split}$
(1)

Unfortunately, for most of the problems we do not have access to the true
distribution $p(x)$ and hence, we need to fit our model to some empirically
observed subset. One solution is to use Monte Carlo sampling to approximate
the integral over $z$ to try to estimate the model parameters $\theta$.
Nonetheless, this approach does not scale to high dimensions of $z$ and
consequently, we will suffer from issues associated with the curse of
dimensionality. Another solution is to use variational inference, e.g., VAE
[8]. In particular, the lower bound of the log-likelihood function, called the
Evidence Lower BOund (ELBO). The ELBO provides a joint optimization objective,
which simultaneously updates the variational posterior $q_{\phi}(z|x)$ and
likelihood model $p_{\theta}(x|z)$. The objective is written as:

$\begin{split}\mathrm{log}\,p(x)&\geq\mathbb{E}_{z\sim
q_{\phi}(z|x)}[\mathrm{log}\,p_{\theta}(x|z)]\\\
&-\mathrm{KL}[q_{\phi}(z|x)||p(z)],\end{split}$ (2)

where KL stands for the Kullback-Leibler divergence.

Figure 1:  Scheme of the different latent variable models. (Top) Single latent
variable model. (Center) Hierarchical latent variable model. (Bottom)
Diffusion model.

###  2.2 Hierarchical Latent Variable Models

Once defined a single stochastic layer, it is straightforward to derive
hierarchical extensions. For example, let us consider a latent variable model
with two latent variables $z_{1}$ and $z_{2}$. We can define the joint
distribution $p(x,z_{1},z_{2}$) and marginalizing out the latent variables:

$\begin{split}p_{\theta}(x)=&\int_{z_{1}}\int_{z_{2}}p_{\theta}(x,z_{1},z_{2})\\\
=&\int_{z_{1}}\int_{z_{2}}p_{\theta}(x|z_{1})p_{\theta}(z_{1}|z_{2})p_{\theta}(z_{2}).\end{split}$
(3)

Similar to the single latent model, we can derive the variational
approximation (ELBO) to the true posterior as:

$\begin{split}\mathrm{log}\,p(x)&\geq\mathbb{E}_{z_{1}\sim
q_{\phi}(z_{1}|z_{2})}[\mathrm{log}\,p_{\theta}(x|z_{1})]\\\
&-\mathrm{KL}[q_{\phi}(z_{1}|x)||p_{\theta}(z_{1}|x)]\\\
&-\mathrm{KL}[q_{\phi}(z_{2}|z_{1})||p(z_{2})].\end{split}$ (4)

###  2.3 Diffusion Models

Diffusion models belong to the latent variable family as well. In fact, we can
think of them as a specific realization of a hierarchical latent variable
model, where the inference model111Remember that the inference model relates a
set of observable variables to a set of latent variables, e.g., $q(z|x)$. does
not have learnable parameters. Instead, it is constructed so that the final
latent distribution $q(x_{T})$ converges to a standard Gaussian (where $T$ is
the number of latent variables). The objective function of diffusion models is
written as:

$\begin{gathered}\mathrm{log}\,p(x)\geq\\\ \mathbb{E}_{x_{1:T}\sim
q(x_{1:T}|x_{0})}[\mathrm{KL}(q(x_{T}|x_{0})||p_{\theta}(x_{T}))\\\
+\sum_{t=2}^{T}\mathrm{KL}(q(x_{t-1}|x_{t},x_{0})||p_{\theta}(x_{t-1}|x_{t}))\\\
-\mathrm{log}\,p_{\theta}(x_{0}|x_{1})].\end{gathered}$ (5)

Under certain assumptions, this objective can be further simplified, leading
to the following approximation:

$\begin{split}\mathrm{log}\,p(x)&\gtrsim\sum_{t=2}^{T}\mathrm{KL}(q(x_{t-1}|x_{t},x_{0})||p_{\theta}(x_{t-1}|x_{t}))\\\
&=\sum_{t=2}^{T}||\epsilon-\epsilon_{\theta}(\sqrt{\bar{\alpha}_{t}}x_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon,t)||^{2}.\end{split}$
(6)

Note that we drop the expectation for clarity. The exact derivation can be
found in [2].

##  3 Methodology

In this section, we provide a brief overview of diffusion models formulation.
Note that we do not aim at covering the entire derivations. For a more in-
depth, detailed mathematical description, we refer the reader to [2].

###  3.1 Background

On a high level, diffusion models consist of two parts: forward diffusion and
parametrized reverse. The forward diffusion part can be described as a
process, where Gaussian noise $\epsilon$ is gradually applied to the input
image $x_{0}$ until the image becomes entirely unrecognizable from a normal
distribution $x_{T}\sim\mathcal{N}(0,\mathrm{I})$ ($T$ is the number of
transformation steps). That is to say, at each step of this process, the noise
is incrementally added to the data,
$x_{0}\xrightarrow{+\epsilon}x_{1}\xrightarrow{+\epsilon}...\xrightarrow{+\epsilon}x_{T}$.
This procedure together with the Markov assumption222Markov assumption is used
to describe a model that holds the memoryless property of a stochastic
process. leads to a simple parameterization forward process expressed as:

$\begin{split}q(x_{1:T}|x_{0})&=\prod_{t=1}^{T}q(x_{t}|x_{t-1})\\\
&=\prod_{t=1}^{T}\mathcal{N}(x_{t};\sqrt{1-\beta_{t}}x_{t-1},\beta_{t}\mathrm{I}),\end{split}$
(7)

where the variable $\beta$ defines a fixed variance schedule, chosen such that
$q(x_{T}|x_{0})\approx\mathcal{N}(0,\mathrm{I})$.

Figure 2: Denoising diffusion process. While the Markov chain of the forward
diffusion gradually adds noise to the input (dash arrows), the reverse process
removes it stepwise (solid arrows).

The second part, the parametrized reverse process, represents the data
synthesis. Thus, it undoes the forward diffusion process and performs
iterative denoising. To that end, the reverse process is trained to generate
data by converting random noise into realistic data. Formally, this generative
process is defined as a stochastic process, which iteratively removes noise
from the input images using deep neural networks. Starting with the pure
Gaussian noise $p(x_{T})=\mathcal{N}(x_{T},0,\mathrm{I})$, the model learns
the joint distribution $p_{\theta}(x_{0:T})$ as:

$\begin{split}p_{\theta}(x_{0:T})&=p(x_{T})\prod_{t=1}^{T}p_{\theta}(x_{t-1}|x_{t})\\\
&=p(x_{T})\prod_{t=1}^{T}\mathcal{N}(x_{t-1};\mu_{\theta}(x_{t},t),\Sigma_{\theta}(x_{t},t)),\end{split}$
(8)

where the time-dependent parameters of the Gaussian transformations $\theta$
are learned. Note in particular that the Markov formulation asserts that a
given reverse diffusion transformation distribution depends only on the
previous timestep.

###  3.2 Training

A diffusion model is trained by finding the reverse Markov transitions that
maximize the likelihood of the training data. In practice, this process
consists of optimizing the variational lower bound on the log likelihood.
Hereunder the simplified expression derived by [2]:

$\begin{split}\mathrm{log}\,p(x)&\gtrsim\sum_{t=2}^{T}||\epsilon-\epsilon_{\theta}(\sqrt{\bar{\alpha}_{t}}x_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon,t)||^{2},\end{split}$
(9)

where

$\displaystyle\alpha_{t}=1-\beta_{t}\;\mathrm{and}\;\bar{\alpha}_{t}=\prod_{i=1}^{T}\alpha_{i}.$
(10)

Note, ultimately, the deep neural network learns to predict the noise
component $\epsilon$ at any given timestep.

##  4 Experiments

In this section, we validate the flexibility of diffusion models for different
seismic tasks. In particular, we analyse three case studies: demultiple,
denoising and interpolation. To do that, we present an end-to-end deep-
learning approach that can deal (separately) with demultiple, denoising and
interpolation scenarios. Furthermore, we benchmark the results with
alternative paradigms that are currently employed in both academia and
industry domains.

The implementation details are as following: In all our experiment, we train
the diffusion model for 200,000 iterations with a batch size of 32; we set
$\beta$ to follow a linear schedule, and we use a depth of 2000 timesteps for
both the forward process (see Equation 7) and the reverse denoising process
(see Equation 8).

Figure 3: In each reverse step $t$, the model $\epsilon_{\theta}$ is fed with
the semi-denoised multiple-free image $x_{t}$ and the multiple-infested input.
As an output, the network generates the image $x_{t-1}$, which should have
less noise and no multiples.
Figure 4: This figure displays two cropped gathers that contain multiples
(input), and the results after applying the demultiple algorithms. Moreover,
we plot the difference between the input and the output to check the content
that has been removed. Note that we apply a scaling factor of 3 in the
differences to stress the changes.

###  4.1 Architecture

Image diffusion models commonly employ a time-conditional U-net [20],
parametrized as $\epsilon_{\theta}(\circ,t)$, as a neural backbone. This
architecture was initially introduced in [2], where the main motivation for
this topology choice was the requirement for the model to have identical input
and output dimensionality. The architecture consists of a stack of residual
layers and downsampling convolutions, followed by a stack of residual layers
with upsampling convolutions; skip connections connect the layers with the
same spatial size. Furthermore, it uses a global attention layer with a single
head to add a projection of the timestep embedding into each residual block.

###  4.2 Demultiple

Primary seismic reflections are events which have reflected only once, and
they are employed to describe the subsurface interfaces. Multiples, on the
contrary, are events which appear when the signal has not taken a direct path
from the source to the receiver after reflecting on a subsurface boundary. The
presence of multiples in a recorded dataset can trigger erroneous
interpretations, since they do not only interfere with the analysis in the
post-stack domain, e.g., stratigraphic interpretation, but also with the
analysis in the pre-stack domain, e.g., amplitude variation with offset
inversion. Thereby, the demultiple process plays a crucial role in any seismic
processing workflow.

In this first experiment, we follow the approach from [21, 18], and generate
synthetic pairs of multiple-infested and multiple-free gathers. This data
setup allows us to train the model in a supervised manner and therefore, we
can frame the demultiple problem as an image-to-image transformation task,
where the network learns to remove the multiples without removing primary
energy. As in [18], the training dataset is designed to include a rich amount
of features present in real datasets, to maximize transferability to real case
uses. To that end, we employ as a baseline a conditional diffusion models
proposed by [22]. More specifically, we condition our model by concatenating
the semi-denoised multiple-free image $x_{t}$ with the multiple-infested input
(see Figure 3). Ideally, the network should return an improved semi-denoised
multiple-free gather $x_{t-1}$ that after $T$ reverse steps should converge
into a noise- and multiple-free gather $x_{0}$.

Once the model is trained, it is crucial to assess the inference capabilities
of the network when working on real data, i.e., generalizability. Nonetheless,
this is not a granted property in deep-learning models due to the distribution
gap between different datasets, e.g., the gap between synthetic and real
datasets [23]. In our experiments, we test the diffusion approach on the
dataset from the Volve field made available under Equinor Open Data Licence.
Furthermore, we compare the outcomes with two other multiple-attenuation
methodologies: one based on Radon-transform [24] and one based on deep
learning [18]. Figure 4 shows an example of such a comparison, where we can
observe how the diffusion solution offers competitive results, despite minimal
hyperparameter tuning involved. For additional results, see Figure 9 in the
Appendix.

Figure 5: Mean and standard deviation of SSIM and SNR metrics calculated on
500 random denoised images. Results from diffusion and FX-Decon scenarios.
Figure 6: This figure shows an example of denoising. The first row contains
the original image and the input image (original with noise). The second row
presents the diffusion and the FX-Decon results. Finally, the third and four
rows display the difference between the results and the original and the input
data, respectively. Figure 7: This figure shows an example of interpolation.
From left to right: the original image, the mask, the input image (original
with mask), the diffusion result and its difference with respect to the
original image.

###  4.3 Denoising

Incoherent noise can be caused by superposition of numerous unwanted signals
from various sources such as ocean waves, wind and electrical instrument noise
among others. Removing such incoherent noise can improve the overall signal-
to-noise ratio and, consequently, increase the certainty of interpretation.
Traditional approaches can be subdivided into two main categories: the
prediction filtering methods and domain transform methods. The first type
assumes linearity and predictability of the signal, and constructs a
predictive filter to suppress the noise [25, 26]. These methods have been
widely adopted by the industry due to their efficiency, although they tend to
under-suppress noise and occasionally suffer from signal leakage [27]. The
second type of methods uses mathematical transformations, e.g., Fourier
transform [28], wavelet transform [29], curvelet transform [30, 31], to steer
the seismic data into domains, where seismic signals and noise can be easier
separated and then leverage the sparse characteristics of seismic data. This
approach, however, often requires a time-consuming transform coefficient
tuning. To cope with this drawback, a new trend based on deep-learning
algorithms has emerged, resulting in optimized solutions that remove
incoherent noise from seismic data as well as speed up the inference time [32,
33].

Similar to the demultiple scenario, we create pairs of images to train our
diffusion model. Nonetheless, this time, the objective is to eliminate
undesired uncorrelated noise, while preserving the inherent characteristics of
the data. To that end, the pairs of training data consist of a real image and
their noisy version. To create the noisy images, we synthetically add Gaussian
noise to the original real images with a variability of the 50% of their
energy. For this second case of study, we train on 1994 BP [34] dataset, from
which we extract random patches (from different shot gathers) that neither
overlap among each other, nor have more than 40% of their content equal to 0.
In this fashion, we try to guarantee certain level of variety in the training
data.

For the testing set, we apply the same conditions as for training.
Additionally, we employ a second dataset (Model94 [34]) to evaluate the
generalization capacity of our system. As for comparison, we use a spectral
filtering technique based on the Fourier transform, namely a complex Wiener
prediction filter called FX-Decon [25, 26], which is dedicated for signal
extraction and non-coherent noise suppression in the frequency domain. To
assess the results, we use structural similarity index (SSIM) and signal-to-
noise ratio (SNR) as quantitative metrics. Figure 5 displays them for each
configuration, i.e., different datasets and methods, and we can observe how
the diffusion model provides the best scores when we test on data coming from
the same dataset as the one used for training. However, as expected, it has a
drastic drop when we test on a new dataset, e.g., Model94. This phenomenon is
mainly caused by the distribution gap between different datasets. On the other
hand, FX-Decon achieves similar performance on both datasets (no drop), as
this method does not involve any learning, i.e., data fitting. Finally, Figure
6 illustrates a denoising example for both algorithms. The difference between
the outputs and the original data (third row in Figure 6) allows us to see
that diffusion model removes some coherent signal, while FX-Decon does not.
Ideally, this should be corrected, but we leave this improvement for future
work. Nevertheless, overall, the diffusion approach leads to less noisy
outputs, as can be noticed in the output image. For additional results, see
Figure 10 in the Appendix.

###  4.4 Interpolation

Seismic data processing algorithms greatly benefit from regularly sampled and
reliable data. However, it is rarely the case where the acquired data is
presented flawless, i.e., complete shot gathers without missing traces.
Frequently, the reason for that are acquisition constraints such as geophones
issues, topography, and economical limitations. As a consequence,
interpolation techniques are a fundamental key for most seismic processing
systems.

Figure 8: Mean and standard deviation of SSIM and SNR metrics calculated on
500 random interpolated images. Results from diffusion and U-net scenarios.

In this last case of study, we evaluate the capacity of our diffusion model to
interpolate missing traces. To that end, we follow the evaluation methodology
introduced by [35], namely, we consider the scenario with irregular missing
traces and with a level of decimation set to 50% (see Figure 7). Regarding the
data for this experiment, we repeat the setup presented in the denoising
section, using 1994 BP dataset for training and testing, and Model94 for
testing on a new dataset. Finally, to have a baseline to compare with, we
implement the so-called “standard” topology from [35], which is essentially a
U-net-like network.

Figure 8 shows the qualitative evaluation of the diffusion approach and of the
U-net baseline. Although results from the latter are superior, the improvement
could be considered marginal given the small metric differences. Furthermore,
both algorithms seem to struggle when inferring on unseen datasets. On the
other hand, besides the quantitative results, the potential that diffusion
models might bring is objectively higher than discriminative models, as the
former are generative models and therefore, can capture more advanced data
properties. For additional results, see Figure 11 in the Appendix.

##  5 Discussion

In this work, we propose a generative framework based on diffusion models to
address several seismic tasks. In particular, our case studies include
demultiple, denoising and interpolation. To solve them, we define the problem
as an image-to-image transformation, where we have an input image that
requires certain modifications so that, the output result belongs to the
target domain. For example, in the demultiple scenario, given a multiple-
infested gather (input domain), our diffusion approach has to identify the
multiples and cancel them out, leading to a multiple-free output gather
(target domain).

The results of our experimental evaluations are fairly encouraging, as they
show competitive performance, when comparing with standardized, customized
algorithms. As we pointed out before, diffusion models for seismic data is an
unexplored field to date and hence, the ultimate goal of this project is not
to outperform these current algorithms in their respective areas, but to
provide a solid analysis of the applicability and flexibility of this novel
framework. Therefore, the main success of our implementation can be regarded
as proof of concept that can be used to adopt generative models, namely
diffusion models, in the geoscience community. We believe that our work can
help to lay the foundation for future research that can benefit both academia
and industry.

##  6 Acknowledgement

This work was developed in the Fraunhofer Cluster of Excellence Cognitive
Internet Technologies. The authors would like to acknowledge the members of
the Fraunhofer ITWM DLSeis consortium (http://dlseis.org) for their financial
support. We are also grateful to Equinor and Volve Licence partners for
releasing Volve seismic field data under an Equinor Open Data Licence.

## References

  * [1] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli, “Deep unsupervised learning using nonequilibrium thermodynamics,” in International Conference on Machine Learning, pp. 2256–2265, PMLR, 2015.
  * [2] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” Advances in Neural Information Processing Systems, vol. 33, pp. 6840–6851, 2020.
  * [3] P. Dhariwal and A. Nichol, “Diffusion models beat gans on image synthesis,” Advances in Neural Information Processing Systems, vol. 34, pp. 8780–8794, 2021.
  * [4] C. Saharia, W. Chan, H. Chang, C. A. Lee, J. Ho, T. Salimans, D. J. Fleet, and M. Norouzi, “Palette: Image-to-image diffusion models,” arXiv preprint arXiv:2111.05826, 2021.
  * [5] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer, “High-resolution image synthesis with latent diffusion models,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 10684–10695, 2022.
  * [6] A. Ramesh, P. Dhariwal, A. Nichol, C. Chu, and M. Chen, “Hierarchical text-conditional image generation with clip latents,” arXiv preprint arXiv:2204.06125, 2022.
  * [7] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” Advances in neural information processing systems, vol. 27, 2014.
  * [8] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” arXiv preprint arXiv:1312.6114, 2013.
  * [9] Y. An, J. Guo, Q. Ye, C. Childs, J. Walsh, and R. Dong, “Deep convolutional neural network for automatic fault recognition from 3d seismic datasets,” Computers & Geosciences, vol. 153, p. 104776, 2021.
  * [10] X. Wu, L. Liang, Y. Shi, and S. Fomel, “FaultSeg3D: using synthetic datasets to train an end-to-end convolutional neural network for 3D seismic fault segmentation,” GEOPHYSICS, vol. 84, no. 3, pp. IM35–IM45, 2019.
  * [11] S. Oh, K. Noh, D. Yoon, S. J. Seol, and J. Byun, “Salt delineation from electromagnetic data using convolutional neural networks,” IEEE Geoscience and Remote Sensing Letters, vol. 16, no. 4, pp. 519–523, 2018.
  * [12] Y. Shi, X. Wu, and S. Fomel, “Saltseg: Automatic 3d salt segmentation using a deep convolutional neural network,” Interpretation, vol. 7, no. 3, pp. SE113–SE122, 2019.
  * [13] P. Nivlet, R. Smith, and N. AlBinHassan, “Automated well-to-seismic tie using deep neural networks,” in SEG Technical Program Expanded Abstracts 2020, pp. 2156–2160, Society of Exploration Geophysicists, 2020.
  * [14] V. Tschannen, A. Ghanim, and N. Ettrich, “Partial automation of the seismic to well tie with deep learning and bayesian optimization,” Computers & Geosciences, vol. 164, p. 105120, 2022.
  * [15] L. Yang and S. Z. Sun, “Seismic horizon tracking using a deep convolutional neural network,” Journal of Petroleum Science and Engineering, vol. 187, p. 106709, 2020.
  * [16] V. Tschannen, M. Delescluse, N. Ettrich, and J. Keuper, “Extracting horizon surfaces from 3d seismic data using deep learning,” Geophysics, vol. 85, no. 3, pp. N17–N26, 2020.
  * [17] A. J. Bugge, A. K. Evensen, J. E. Lie, and E. H. Nilsen, “Demonstrating multiple attenuation with model-driven processing using neural networks,” The Leading Edge, vol. 40, no. 11, pp. 831–836, 2021.
  * [18] R. Durall, A. Ghanim, N. Ettrich, and J. Keuper, “Dissecting u-net for seismic application: An in-depth study on deep learning multiple removal,” arXiv preprint arXiv:2206.12112, 2022.
  * [19] J. M. Tomczak, Deep Generative Modeling.  Springer Nature, 2022.
  * [20] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in International Conference on Medical image computing and computer-assisted intervention, pp. 234–241, Springer, 2015\.
  * [21] A. Breuer, N. Ettrich, and P. Habelitz, “Deep learning in seismic processing: Trim statics and demultiple,” in SEG Technical Program Expanded Abstracts 2020, pp. 3199–3203, Society of Exploration Geophysicists, 2020.
  * [22] C. Saharia, J. Ho, W. Chan, T. Salimans, D. J. Fleet, and M. Norouzi, “Image super-resolution via iterative refinement,” arXiv preprint arXiv:2104.07636, 2021.
  * [23] R. Durall, V. Tschannen, F.-J. Pfreundt, and J. Keuper, “Synthesizing seismic diffractions using a generative adversarial network,” in SEG Technical Program Expanded Abstracts 2020, pp. 1491–1495, Society of Exploration Geophysicists, 2020.
  * [24] D. Hampson, “Inverse velocity stacking for multiple elimination,” in SEG Technical Program Expanded Abstracts 1986, pp. 422–424, Society of Exploration Geophysicists, 1986.
  * [25] N. Gulunay, “Fx decon and the complex wiener prediction filter for random noise reduction on stacked data,” in 56th SEG annual meeting, pp. 2–6, 1986.
  * [26] M. Galbraith, “Random noise attenuation by fx prediction: A tutorial,” SEG Technical Program Expanded Abstracts 1991, pp. 1428–1431, 1991.
  * [27] N. Gülünay, “Signal leakage in f-x deconvolution algorithms,” Geophysics, vol. 82, no. 5, pp. W31–W45, 2017.
  * [28] M. Naghizadeh, “Seismic data interpolation and denoising in the frequency-wavenumber domain,” Geophysics, vol. 77, no. 2, pp. V71–V80, 2012.
  * [29] S. M. Mousavi, C. A. Langston, and S. P. Horton, “Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform,” Geophysics, vol. 81, no. 4, pp. V341–V355, 2016.
  * [30] R. Neelamani, A. I. Baumstein, D. G. Gillard, M. T. Hadidi, and W. L. Soroka, “Coherent and random noise attenuation using the curvelet transform,” The Leading Edge, vol. 27, no. 2, pp. 240–248, 2008.
  * [31] F. J. Herrmann, D. Wang, G. Hennenfent, and P. P. Moghaddam, “Curvelet-based seismic data processing: A multiscale and nonlinear approach,” Geophysics, vol. 73, no. 1, pp. A1–A5, 2008.
  * [32] S. Yu, J. Ma, and W. Wang, “Deep learning for denoisingdeep learning for denoising,” Geophysics, vol. 84, no. 6, pp. V333–V350, 2019.
  * [33] O. M. Saad and Y. Chen, “Deep denoising autoencoder for seismic random noise attenuation,” Geophysics, vol. 85, no. 4, pp. V367–V376, 2020.
  * [34] S. H. Gray and K. J. Marfurt, “Migration from topography: Improving the near-surface image,” Canadian Journal of Exploration Geophysics, vol. 31, no. 1-2, pp. 18–24, 1995.
  * [35] M. Fernandez, R. Durall, N. Ettrich, M. Delescluse, A. Rabaute, and J. Keuper, “A comparison of deep learning paradigms for seismic data interpolation,” in Second EAGE Digitalization Conference and Exhibition, vol. 2022, pp. 1–5, European Association of Geoscientists & Engineers, 2022.

##  7 Appendix

###  7.1 More Results

We provide additional results, where we can visualize the evolution of the
reverse process for all the aforementioned case studies (see Figure 9, Figure
10 and Figure 11). Note that the subindexes of the $x$ indicate the output of
an intermediate step during the inference process, being $x_{1999}$ random
noise and $x_{0}$ the final output of the model.

Figure 9: This figure displays demultiple results at different intermediate
steps for the reverse process. Note that the first two rows show synthetic
data examples, while the last two from the Volve dataset.

Figure 10: This figure displays denoising results at different intermediate
steps for the reverse process. Note that the examples belong to the Model94
dataset.

Figure 11: This figure displays interpolation results at different
intermediate steps for the reverse process. Note that the examples belong to
the Model94 dataset.