# Phase-field fracture irreversibility using the slack variable approach

Ritukesh Bharali  Department of Industrial and Material Science, Chalmers
University of Technology  Fredrik Larsson  Department of Industrial and
Material Science, Chalmers University of Technology  Ralf Jänicke  Institute
of Applied Mechanics, Technische Universität Braunschweig

###### Abstract

In this manuscript, the phase-field fracture irreversibility constraint is
transformed into an equality-based constraint using the slack variable
approach. The equality-based fracture irreversibility constraint is then
introduced in the phase-field fracture energy functional using the Lagrange
Multiplier Method and the Penalty method. Both methods are variationally
consistent with the conventional variational inequality phase-field fracture
problem, unlike the history-variable approach. Thereafter, numerical
experiments are carried out on benchmark problems in brittle and quasi-brittle
fracture to demonstrate the efficacy of the proposed method.

## Highlights

  * •

Effective equality constraint for fracture irreversibility inequality
constraint using slack variable.

  * •

The Lagrange Multiplier Method and the Penalty method are adopted to augment
the energy functional.

  * •

Numerical experiments carried out on both brittle and quasi-brittle fracture
problems.

Keywords:  phase-field fracture, brittle, quasi-brittle, Lagrange multiplier
method, Penalty method, COMSOL

##  1 Introduction

The phase-field model for fracture emerged from the seminal work of [1],
wherein the Griffith fracture criterion was cast into a variational setting.
Later, a numerical implementation of the same was proposed in [2], using the
Ambrosio–Tortorelli regularisation of the Mumford-Shah potential [3]. In this
implementation, the fracture is represented by an auxiliary variable, that
interpolates between the intact and broken material states. Such a formulation
allows the automatic tracking of fractures on a fixed mesh, thereby
eliminating the need for the tedious tracking and remeshing processes observed
with discrete methods. Furthermore, the phase-field model for fracture is also
able to handle topologically complex (branching, kinking and merging)
fractures, and is able to demonstrate fracture initiation without introducing
any singularity [4]. Owing to these advantages, the phase-field model for
fracture has gained popularity in the computational mechanics community in the
past decade.

A thermodynamically consistent formulation of the phase-field fracture model
was proposed by [5], adopting an energetic cracking driving force definition.
Since then, over the past decade, researchers have extended the work to
thermo-mechanical problems at large strain [6], ductile failure [7, 8],
fracture in thin films [9], anisotropic fracture [10], fracture in
fully/partially saturated porous media [11, 12, 13], hydrogen assisted
cracking [14], dissolution-driven stress corrosion cracking [15], fracture and
fatigue in shape-memory alloys [16], and brittle failure of Reisner-Mindlin
plates [17] to cite a few. Finally, the phase-field model for fracture has
also been used for the multi-scale finite element method [18, 19, 20],
asymptotic homogenisation [21], and variationally consistent homogenisation
[22]. For a detailed overview on the phase-field fracture model, the reader is
referred to the comprehensive review works [23, 24, 25].

Despite its growing popularity, the phase-field model for fracture poses
several challenges when it comes to robust and computationally efficient
solution techniques [25]. These include,

  1. (A)

the poor performance of the monolithic Newton-Raphson (NR) method due to the
non-convex energy functional,

  2. (B)

variational inequality arising from fracture irreversibility, and

  3. (C)

extremely fine meshes are required in the fracture zone.

In order to alleviate the problem concerning the poor-performance of the NR
method (A), [26] utilised a linear extrapolation in time111The linear
extrapolation of the phase-field variable in time is a questionable assumption
although it results in a robust solution method. for the phase-field variable
in the momentum balance equation. In [27], a novel line search technique was
developed, while in [28] and [29] Augmented Lagrangian and modified NR
approaches were proposed. More recently, [30] adopted a recursive multilevel
trust region method, that resulted in improved convergence of the NR method.
The development of robust, monolithic solution technique is still an active
area in phase-field fracture research, and this manuscript is contribution
towards this aspect.

The next issue pertains to the variational inequality problem arising out of
the fracture irreversibility condition (B). This aspect has been treated in
different ways in the phase-field fracture literature. They include crack set
irreversibility [2, 31, 32], penalisation [4, 28, 29], and the implicit
History variable based method proposed by [33]. The lattermost method remains
popular despite its non-variational nature, over-estimation of the fracture
surface energy, and the necessity for computationally expensive alternate
minimisation (staggered) solution procedure [4]. In a novel approach, this
manuscript adopts the slack variable method [34, 35], that constructs an
equivalent variational equality problem, while maintaining the variational
structure of the original problem.

Another issue pertaining to the phase-field fracture model is the requirement
of extremely fine meshes in the smeared fracture zone (C). In this context,
[2] and [9] advocates the use of uniformly refined meshes together with
parallel computing while [36, 37, 38] opted for error-controlled adaptive mesh
refinement. Some other approaches include pre-refined meshes when the fracture
path is known [23, 27, 4], adaptive mesh refinement in [26] based on a phase-
field threshold or Kelly error estimates [39], and multi-level hp-FEM in [40].
In this manuscript, uniformly refined meshes and pre-refined meshes based on a
priori knowledge of the fracture path, are used.

The novelty of this manuscript lies in the alternative treatment of the
fracture irreversibility inequality constraint, using the slack variable
approach. The inequality constraint is replaced by an equivalent equality-
based fracture irreversibility constraint. The constraint is then augmented to
the phase-field fracture energy functional using the Lagrange Multiplier
method and the Penalty method. These methods are variationally consistent,
unlike the history-variable approach proposed in [33].

The manuscript is organised as follows: Section 2 introduces the phase-field
model for fracture, its underlying energy functional and the Karush-Kuhn-
Tucker (KKT) optimality conditions. Thereafter, in Section 3, the Method of
Multipliers is adopted to augment the energy functional a slack variable-based
fracture irreversibility criterion. The equivalence of the new formulation
with the original problem in Section 2 is established in terms of KKT
conditions. Furthermore, the Euler-Lagrange equations are presented. In
Section 4, the slack variable-based fracture irreversibility constraint is
introduced in the energy functional using the Penalisation method, followed by
the derivation of the pertinent Euler-Lagrange equations. The numerical
experiments on benchmark problems are presented in Section 5, and Section 6
lays down the concluding remarks of this manuscript.

## Notation

The following notations are strictly adhered to in this manuscript:

  * •

Zero-order tensors are represented using small italicized letters, e.g., $a$.
Bold symbols are used for first and second-order tensors, for instance, stress
$\sigma$ and strain $\epsilon$.

  * •

A function $f$ with its arguments $x$, $y$ is written in the form $f(x,y)$,
whereas a variable $g$ with operational dependencies $p,q$ is written as
$g[p,q]$.

  * •

The Macaulay operator on a variable $x$ is defined as $\langle
x\rangle_{\pm}=\frac{1}{2}(x\pm|x|)$.

##  2 Phase-field fracture model

###  2.1 Energy functional

Let $\hbox{$\Omega$}\in\mathbb{R}^{\text{dim}}$ ($\text{dim}=2,3$) be the
domain occupied by the fracturing solid, shown in Figure 1. Its boundary
$\Gamma$ is decomposed into a Dirichlet boundary
$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$ and a Neumann boundary
$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$, such that
$\hbox{$\Gamma$}=\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\cup\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}$
and
$\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\cap\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}=\emptyset$.
Furthermore, the fracture is represented by an auxiliary variable (phase-
field) $\hbox{$\varphi$}\in[0,1]$ within a diffusive (smeared) zone of width
$\hbox{$l$}>0$.

$\Omega$$l$$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$
Figure 1: A solid $\hbox{$\Omega$}\in\mathbb{R}^{2}$ embedded with a diffused
(smeared) crack. Dirichlet and Neumann boundaries are indicated with
$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$ and
$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$ respectively. Figure adopted from [22].

The energy functional of the phase-field fracture model [27, 4] is stated as,

${}\displaystyle
E(\hbox{$\mathbf{u}$},\hbox{$\varphi$})=\int_{\hbox{$\Omega$}}g(\hbox{$\varphi$})\Psi^{+}(\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])\>\text{d}\hbox{$\Omega$}+\int_{\hbox{$\Omega$}}\Psi^{-}(\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])\>\text{d}\hbox{$\Omega$}-\int_{\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}}\hbox{$\mathbf{t}_{p}^{\scriptsize{u}}$}\cdot\hbox{$\mathbf{u}$}\>\text{d}\hbox{$\Gamma$}+\int_{\hbox{$\Omega$}}\dfrac{\hbox{$G_{c}$}}{c_{w}\hbox{$l$}}\left(w(\hbox{$\varphi$})+\hbox{$l$}|\hbox{\boldmath$\nabla$}\hbox{$\varphi$}|^{2}\right)\>\text{d}\hbox{$\Omega$},$
(1)

with $\mathbf{u}$ and $\varphi$ representing the displacement and the phase-
field respectively. The symbols adopted for the additional functions and
parameters in (1) are presented in Table 1.

Symbol | Function/Parameter
---|---
$g(\hbox{$\varphi$})$ | Degradation function
$\Psi^{+}$ | Fracture driving strain energy density
$\Psi^{-}$ | Residual strain energy density
$\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}]$ |  Symmetric gradient of $\mathbf{u}$
$\mathbf{t}_{p}^{\scriptsize{u}}$ | Prescribed traction
$G_{c}$ | Griffith fracture energy
$c_{w}$ | Normalisation constant
$l$ | Fracture length-scale
$w(\hbox{$\varphi$})$ | Local fracture energy function
Table 1: List of functions and parameters corresponding to the phase-field
fracture energy functional (1).

In this manuscript, $\Psi^{+}$ and $\Psi^{-}$ represent the tensile and the
compressive strain energy density functions respectively. They are functions
of strain tensor $\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}]$,
given by

${}\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}]=(\hbox{$\mathbf{u}$}\otimes\hbox{\boldmath$\nabla$})^{\text{sym}}.$
(2)

The explicit expression for $\Psi^{\pm}$ is obtained following [5] as,

${}\Psi^{\pm}(\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])=\frac{1}{2}\IfEqCase{1}{{1}{\lambda}{2}{\mu}}[]\langle(\mathbf{I}\colon\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])\rangle^{2}_{\pm}+\IfEqCase{2}{{1}{\lambda}{2}{\mu}}[]\hbox{$\hbox{\boldmath$\epsilon$}$}^{\pm}[\hbox{$\mathbf{u}$}]\colon\hbox{$\hbox{\boldmath$\epsilon$}$}^{\pm}[\hbox{$\mathbf{u}$}],$
(3)

where $\IfEqCase{1}{{1}{\lambda}{2}{\mu}}[]$ and
$\IfEqCase{2}{{1}{\lambda}{2}{\mu}}[]$ are the Lamé constants, $\mathbf{I}$ is
a second-order identity tensor, $\mathbf{u}$ is the displacement. Moreover,
the tensile/compressive strain
$\hbox{$\hbox{\boldmath$\epsilon$}$}^{\pm}[\hbox{$\mathbf{u}$}]$ is defined as

$\hbox{$\hbox{\boldmath$\epsilon$}$}^{\pm}[\hbox{$\mathbf{u}$}]=\sum_{i=1}^{\text{dim}}\langle\epsilon_{i}\rangle_{\pm}\mathbf{p}_{i}\otimes\mathbf{p}_{i}$
(4)

where $\epsilon_{i}$ represents the $i^{th}$ eigenvalue of the strain, and
$\mathbf{p}_{i}$ is its corresponding normalised eigenvector. The
tensile/compressive strain energy definition in (3) yield the corresponding
stresses upon taking the derivative w.r.t the strain. They are explicitly
stated in a compact form as,

$\hbox{$\hbox{\boldmath$\sigma$}$}^{\pm}=\frac{\partial\Psi^{\pm}}{\partial\hbox{$\hbox{\boldmath$\epsilon$}$}}=\IfEqCase{1}{{1}{\lambda}{2}{\mu}}[]\langle(\mathbf{I}\colon\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])\rangle_{\pm}\mathbf{I}+2\IfEqCase{2}{{1}{\lambda}{2}{\mu}}[]\hbox{$\hbox{\boldmath$\epsilon$}$}^{\pm}[\hbox{$\mathbf{u}$}].$
(5)

The last integral in the energy functional (1) represents the fracture energy
in a phase-field fracture model. The phase-field fracture model allows
flexibility in the choice of the degradation function $g(\hbox{$\varphi$})$
and the locally dissipated fracture energy function $w(\hbox{$\varphi$})$.
Table 2 presents commonly adopted functions for brittle and quasi-brittle
fracture.

Model | $g(\hbox{$\varphi$})$ | $w(\hbox{$\varphi$})$ | $c_{w}$
---|---|---|---
Brittle-AT1[41] | $(1-\hbox{$\varphi$})^{2}$ | $\varphi$ | $8/3$
Brittle-AT2[31] | $(1-\hbox{$\varphi$})^{2}$ | $\hbox{$\varphi$}^{2}$ | $2$
Quasi-Brittle[42] | $\dfrac{(1-\hbox{$\varphi$})^{p}}{(1-\hbox{$\varphi$})^{p}+a_{1}\hbox{$\varphi$}+a_{1}a_{2}\hbox{$\varphi$}^{2}+a_{1}a_{2}a_{3}\hbox{$\varphi$}^{3}}$ | $2\hbox{$\varphi$}-\hbox{$\varphi$}^{2}$ | $\pi$
Table 2: Degradation function, local fracture energy function and its
normalisation constant for different phase-field fracture models. The
abbreviation ‘AT’ stands for Ambrosio-Tortorelli.

For the Quasi-Brittle model in Table 2, the parameters $p$, $a_{1}$, $a_{2}$
and $a_{3}$ are chosen such that they mimic the different traction-separation
laws in cohesive zone modelling. Following [42], the constant $a_{1}$ is given
by

$a_{1}=\dfrac{4E_{0}\hbox{$G_{c}$}}{\pi\hbox{$l$}f_{t}^{2}}$ (6)

where, the newly introduced parameters $E_{0}$ and $f_{t}$ represent the
Youngs’ modulus and the tensile strength of the material. Next, Table 3
presents the different traction-separation laws and the corresponding values
of $p$, $a_{2}$ and $a_{3}$ [42].

Model | $p$ | $a_{2}$ | $a_{3}$
---|---|---|---
Linear Softening | $2$ | $-0.5$ | $0$
Exponential Softening | $2.5$ | $2^{5/3}-3$ | $0$
Cornelissen et. al.[43] Softening | $2$ | $1.3868$ | $0.6567$
Table 3: Quasi-brittle phase-field fracture models and their parameters, based
on [42].

The prediction of fracture initiation and propagation in the solid occupying
the domain $\Omega$ requires solving the constrained minimisation problem:

######  Problem 1.

Find $\hbox{$\mathbf{u}$}(t)$ and $\hbox{$\varphi$}(t)$, where
${}^{n}\hbox{$\varphi$}=\hbox{$\varphi$}(t_{n})$ is known for $t_{n}<t$ such
that,

$\begin{array}[]{lll}min_{\tiny{\hbox{$\mathbf{u}$}(t),\hbox{$\varphi$}(t)}}E&\text{and}&h(\hbox{$\varphi$})\geq
0,\end{array}$ (7)

where $h(\hbox{$\varphi$})=\hbox{$\varphi$}-{}^{n}\hbox{$\varphi$}$, and the
left superscript $n$ denotes the previous time-step. Moreover, the problem is
augmented with (pseudo) time-dependent Dirichlet boundary conditions
$\hbox{$\mathbf{u}$}^{p}$ on $\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$ and
$\hbox{$\varphi$}^{p}$ on
$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$, and/or Neumann boundary
conditions $\mathbf{t}_{p}^{\scriptsize{u}}$ on
$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$ and
$\mathbf{t}_{p}^{\scriptsize{\hbox{$\varphi$}}}$ on
$\hbox{$\Gamma$}_{N}^{\scriptsize{\hbox{$\varphi$}}}$. Furthermore, the
boundary $\Gamma$ is split as
$\hbox{$\Gamma$}=\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\cup\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}$,
$\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\cap\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}=\emptyset$
and
$\hbox{$\Gamma$}=\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$}\cup\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{\hbox{$\varphi$}}}$}$,
$\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$}\cap\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{\hbox{$\varphi$}}}$}=\emptyset$
respectively.  $\blacksquare$

###  2.2 Karush-Kuhn-Tucker conditions

The KKT triple ($\hbox{$\mathbf{u}$},\hbox{$\varphi$},\Lambda$) in appropriate
spaces for the Problem 1 requires

$\displaystyle\big{[}E(\hbox{$\mathbf{u}$},\hbox{$\varphi$})-\int_{\hbox{$\Omega$}}\Lambda
h(\hbox{$\varphi$})\,\text{d}\hbox{$\Omega$}\big{]}^{\prime}_{\hbox{$\mathbf{u}$}}$
$\displaystyle=0,$ (8a)
$\displaystyle\big{[}E(\hbox{$\mathbf{u}$},\hbox{$\varphi$})-\int_{\hbox{$\Omega$}}\Lambda
h(\hbox{$\varphi$})\,\text{d}\hbox{$\Omega$}\big{]}^{\prime}_{\hbox{$\varphi$}}$
$\displaystyle=0,,$ (8b) $\displaystyle h(\hbox{$\varphi$})$
$\displaystyle\geq 0,$ (8c) $\displaystyle\Lambda$ $\displaystyle\geq 0,$ (8d)
$\displaystyle\Lambda h(\hbox{$\varphi$})$ $\displaystyle=0$ (8e)

to hold, where $\Lambda$ is a multiplier. Equations (8a) and (8b) represent
the Euler-Lagrange equations pertaining to the stationary condition, while
equations (8c-8e) are the primal and dual feasibility conditions, and the
complementary slackness respectively. Formally, the latter set of equations
enforces fracture irreversbility. To the same end, alternative, but equivalent
formulations have been adopted in the phase-field fracture literature. For
instance, [26] adopted the primal-dual active set strategy proposed by [44],
[28] utilised an Augmented Lagrangian penalisation approach with the Moreau-
Yoshida indicator function, and [4] opted for a simple penalisation approach.
The authors would like to emphasise that the rather popular history-variable
approach in [5] is not variationally consistent.

##  3 Lagrange Multiplier Method (LMM)

###  3.1 Fracture irreversibility and modified energy functional

In order to enforce fracture irreversibility
$h(\hbox{$\varphi$})=\hbox{$\varphi$}-{}^{n}\hbox{$\varphi$}\geq 0$, a slack
variable is defined as

$\hbox{$\theta$}^{2}=h(\hbox{$\varphi$})=\hbox{$\varphi$}-{}^{n}\hbox{$\varphi$}.$
(9)

It is observed that $\hbox{$\theta$}^{2}$ admits value greater than or equal
to zero, thereby fulfilling the fracture irreversibility criterion. Next, the
constrained minimisation Problem 1 is reformulated as:

######  Problem 2.

Find $\hbox{$\mathbf{u}$}(t)$, $\hbox{$\varphi$}(t)$ and $\hbox{$\theta$}(t)$
such that,

$\begin{array}[]{lll}\text{min}_{\scriptsize{\hbox{$\mathbf{u}$}(t),\hbox{$\varphi$}(t)},\hbox{$\theta$}(t)}\>\>E\\\
\text{subjected to}\>\>h(\hbox{$\varphi$})-\hbox{$\theta$}^{2}=0.\end{array}$
(10)

with suitable (pseudo) time-dependent Dirichlet and/or Neumann boundary
conditions, as mentioned in Problem 1.  $\blacksquare$

Augmenting the equality constraint (10) in (1) via the Lagrange multiplier
$\Lambda$ results in the modified energy functional,

${}\tilde{E}(\hbox{$\mathbf{u}$},\hbox{$\varphi$},\theta,\Lambda)=E(\hbox{$\mathbf{u}$},\hbox{$\varphi$})-\Lambda\big{(}h(\hbox{$\varphi$})-\hbox{$\theta$}^{2}\big{)}.$
(11)

###  3.2 Stationary conditions

The stationary conditions for the Problem 2 are given by

$\displaystyle\tilde{E}^{\prime}_{\hbox{$\mathbf{u}$}}(\hbox{$\mathbf{u}$},\hbox{$\varphi$},\hbox{$\theta$},\Lambda;\hbox{$\delta\mathbf{u}$})$
$\displaystyle=0,$ (12a)
$\displaystyle\tilde{E}^{\prime}_{\hbox{$\varphi$}}(\hbox{$\mathbf{u}$},\hbox{$\varphi$},\hbox{$\theta$},\Lambda;\hbox{$\delta\varphi$})$
$\displaystyle=0,$ (12b)
$\displaystyle\tilde{E}^{\prime}_{\hbox{$\theta$}}(\hbox{$\mathbf{u}$},\hbox{$\varphi$},\hbox{$\theta$},\Lambda;\delta\hbox{$\theta$})$
$\displaystyle=0,$ (12c)
$\displaystyle\tilde{E}^{\prime}_{\Lambda}(\hbox{$\mathbf{u}$},\hbox{$\varphi$},\hbox{$\theta$},\Lambda;\delta\Lambda)$
$\displaystyle=0.$ (12d)

For equivalence of these conditions with those presented for Problem 1, the
reader is referred to Proposition 3.1 in [45].

###  3.3 Euler-Lagrange equations

The Euler-Lagrange equations for Problem 2 are (12a)-(12d), obtained upon
taking the first variation of the energy functional (11) w.r.t its solution
variables $\mathbf{u}$, $\varphi$, $\theta$ and $\Lambda$. This results in,

######  Problem 3.

Find ($\mathbf{u}$, $\varphi$, $\theta$, $\Lambda$)
$\in\mathbb{U}\times\mathbb{V}\times\mathbb{W}\times\mathbb{A}$ with

$\displaystyle R^{\mathbf{u}}$
$\displaystyle=\int_{\hbox{$\Omega$}}\bigg{(}g(\hbox{$\varphi$})\hbox{$\hbox{\boldmath$\sigma$}^{+}$}[\hbox{$\mathbf{u}$}]+\hbox{$\hbox{\boldmath$\sigma$}^{-}$}[\hbox{$\mathbf{u}$}]\bigg{)}\colon\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\delta\mathbf{u}$}]\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\hbox{$\delta\mathbf{u}$}\in\mathbb{U}^{0}$ (13a)
$\displaystyle R^{\varphi}$
$\displaystyle=\int_{\hbox{$\Omega$}}\bigg{(}\dfrac{\hbox{$G_{c}$}}{c_{w}\hbox{$l$}}w^{\prime}(\hbox{$\varphi$})+g^{\prime}(\hbox{$\varphi$})\Psi^{+}(\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])-\Lambda\bigg{)}\hbox{$\delta\varphi$}\>\text{d}\hbox{$\Omega$}+\int_{\hbox{$\Omega$}}\dfrac{2\hbox{$G_{c}$}\hbox{$l$}}{c_{w}}\hbox{\boldmath$\nabla$}\hbox{$\varphi$}\cdot\hbox{\boldmath$\nabla$}\hbox{$\delta\varphi$}\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\hbox{$\delta\varphi$}\in\mathbb{V}^{0}$ (13b)
$\displaystyle R^{\theta}$
$\displaystyle=\int_{\hbox{$\Omega$}}2\Lambda\theta\delta\theta\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\delta\theta\in\mathbb{W}$ (13c) $\displaystyle
R^{\Lambda}$
$\displaystyle=-\int_{\hbox{$\Omega$}}\big{(}h(\hbox{$\varphi$})-\hbox{$\theta$}^{2}\big{)}\delta\Lambda\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\delta\Lambda\in\mathbb{A}.$ (13d)

with trial function spaces

$\displaystyle\mathbb{U}:=\\{\hbox{$\mathbf{u}$}\in[H^{1}(\hbox{$\Omega$})]^{\text{dim}}|\hbox{$\mathbf{u}$}=\hbox{$\mathbf{u}$}^{\text{{p}}}\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\\},$
(14a)
$\displaystyle\mathbb{V}:=\\{\hbox{$\varphi$}\in[H^{1}(\hbox{$\Omega$})]^{1}|\hbox{$\varphi$}=\hbox{$\varphi$}^{\text{{p}}}\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$}\\},$
(14b) $\displaystyle\mathbb{W}:=\\{\theta\in[L_{2}(\hbox{$\Omega$})]\\},$
(14c) $\displaystyle\mathbb{A}:=\\{\Lambda\in[L_{2}(\hbox{$\Omega$})]\\},$
(14d)

and test function spaces

$\displaystyle\mathbb{U}^{0}:=\\{\hbox{$\mathbf{u}$}\in[H^{1}(\hbox{$\Omega$})]^{\text{dim}}|\hbox{$\mathbf{u}$}=0\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\\},$
(15a)
$\displaystyle\mathbb{V}^{0}:=\\{\hbox{$\varphi$}\in[H^{1}(\hbox{$\Omega$})]^{1}|\hbox{$\varphi$}=0\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$}\\},$
(15b)

and pertinent (pseudo) time-dependent Neumann condition

$\mathbf{t}:=\left(g(\hbox{$\varphi$})\hbox{$\hbox{\boldmath$\sigma$}^{+}$}+\hbox{$\hbox{\boldmath$\sigma$}^{-}$}\right)\cdot\mathbf{n}=\hbox{$\mathbf{t}_{p}^{\scriptsize{u}}$}\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}.$
(16)

$\blacksquare$

######  Remark 1.

In the event $\Lambda$ and $\theta$ are zero, their stiffness contribution is
zero and the stiffness matrix is singular. In order to avoid this situation a
fictitious stiffness is introduced to guarantee solvability.

##  4 Penalty (Pen.) method

###  4.1 The energy functional

In an alternative approach, the squared slack variable equation (9) is
introduced in the energy functional (1) as a quadratic term via penalisation.
This results in the modified energy functional,

${}\hat{E}(\hbox{$\mathbf{u}$},\hbox{$\varphi$},\theta)=E(\hbox{$\mathbf{u}$},\hbox{$\varphi$})+\dfrac{\eta}{2}\int_{\hbox{$\Omega$}}\bigg{(}h(\hbox{$\varphi$})-\hbox{$\theta$}^{2}\bigg{)}^{2}\,\text{d}\hbox{$\Omega$},$
(17)

where, $\eta$ is a penalty parameter.

###  4.2 Euler-Lagrange equations

The Euler-Lagrange equations are obtained upon taking the first variation of
the energy functional (17) w.r.t its solution variables $\mathbf{u}$,
$\varphi$ and $\theta$. This results in:

######  Problem 4.

Find ($\mathbf{u}$, $\varphi$, $\theta$)
$\in\mathbb{U}\times\mathbb{V}\times\mathbb{W}$ with

$\displaystyle R^{\mathbf{u}}$
$\displaystyle=\int_{\hbox{$\Omega$}}\bigg{(}g(\hbox{$\varphi$})\hbox{$\hbox{\boldmath$\sigma$}^{+}$}[\hbox{$\mathbf{u}$}]+\hbox{$\hbox{\boldmath$\sigma$}^{-}$}[\hbox{$\mathbf{u}$}]\bigg{)}\colon\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\delta\mathbf{u}$}]\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\hbox{$\delta\mathbf{u}$}\in\mathbb{U}^{0}$ (18a)
$\displaystyle R^{\varphi}$
$\displaystyle=\int_{\hbox{$\Omega$}}\bigg{(}\dfrac{\hbox{$G_{c}$}}{c_{w}\hbox{$l$}}w^{\prime}(\hbox{$\varphi$})+g^{\prime}(\hbox{$\varphi$})\Psi^{+}(\hbox{$\hbox{\boldmath$\epsilon$}$}[\hbox{$\mathbf{u}$}])+\eta(h(\hbox{$\varphi$})-\hbox{$\theta$}^{2})\bigg{)}\hbox{$\delta\varphi$}\>\text{d}\hbox{$\Omega$}+\int_{\hbox{$\Omega$}}\dfrac{2\hbox{$G_{c}$}\hbox{$l$}}{c_{w}}\hbox{\boldmath$\nabla$}\hbox{$\varphi$}\cdot\hbox{\boldmath$\nabla$}\hbox{$\delta\varphi$}\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\hbox{$\delta\varphi$}\in\mathbb{V}^{0}$ (18b)
$\displaystyle R^{\theta}$
$\displaystyle=\int_{\hbox{$\Omega$}}-2\eta\theta(h(\hbox{$\varphi$})-\hbox{$\theta$}^{2})\delta\theta\>\text{d}\hbox{$\Omega$}=0\>\>\>$
$\displaystyle\forall\>\delta\theta\in\mathbb{W}$ (18c)

with trial function spaces

$\displaystyle\mathbb{U}:=\\{\hbox{$\mathbf{u}$}\in[H^{1}(\hbox{$\Omega$})]^{\text{dim}}|\hbox{$\mathbf{u}$}=\hbox{$\mathbf{u}$}^{\text{{p}}}\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\\},$
(19a)
$\displaystyle\mathbb{V}:=\\{\hbox{$\varphi$}\in[H^{1}(\hbox{$\Omega$})]^{1}|\hbox{$\varphi$}=\hbox{$\varphi$}^{\text{{p}}}\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$}\\},$
(19b) $\displaystyle\mathbb{W}:=\\{\theta\in[L_{2}(\hbox{$\Omega$})]|\\},$
(19c)

and test function spaces

$\displaystyle\mathbb{U}^{0}:=\\{\hbox{$\mathbf{u}$}\in[H^{1}(\hbox{$\Omega$})]^{\text{dim}}|\hbox{$\mathbf{u}$}=0\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{u}}$}\\},$
(20a)
$\displaystyle\mathbb{V}^{0}:=\\{\hbox{$\varphi$}\in[H^{1}(\hbox{$\Omega$})]^{1}|\hbox{$\varphi$}=0\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{D}^{\scriptsize{\hbox{$\varphi$}}}$}\\},$
(20b)

and pertinent (pseudo) time-dependent Neumann condition

$\mathbf{t}:=\left(g(\hbox{$\varphi$})\hbox{$\hbox{\boldmath$\sigma$}^{+}$}+\hbox{$\hbox{\boldmath$\sigma$}^{-}$}\right)\cdot\mathbf{n}=\hbox{$\mathbf{t}_{p}^{\scriptsize{u}}$}\>\>\text{on}\>\>\hbox{$\hbox{$\Gamma$}_{N}^{\scriptsize{u}}$}.$
(21)

$\blacksquare$

######  Remark 2.

In the event $h(\hbox{$\varphi$})-\hbox{$\theta$}^{2}$ and $\theta$ are zero,
their stiffness contribution is zero and the stiffness matrix is singular. In
order to avoid this situation a fictitious stiffness is introduced to
guarantee solvability.

##  5 Numerical Study

The numerical experiments on benchmark problems are presented in this section.
They include both brittle and quasi-brittle fracture. In the brittle fracture
domain, the Single Edge Notched Specimen under Tension (SENT), and under Shear
(SENS), and the Notched concrete specimen with a hole [23] is presented. For
quasi-brittle fracture, the benchmark problems include the concrete three-
point bending problem [46], and the Winkler L-shaped panel experiment [47].
For each problem, the geometry as well as the material properties and loading
conditions are presented in the respective sub-section. The load-displacement
curves and the phase-field topology at the final step of analysis is also
presented.

For the penalty method, the penalty parameter $\eta$ is set to $10^{6}$
[N/mm2] throughout the study. This choice is motivated by similarity of the
load-displacement curves obtained with those from the Lagrange Multiplier
Method. Furthermore, for both the aforementioned methods, a fully coupled
monolithic solution technique is adopted using the Newton-Raphson solver in
COMSOL Multiphysics. The iterations are terminated when the error measure
defined as the weighted Euclidean norm of the solution update,

$err:=\sqrt{\dfrac{1}{M}}\sqrt{\sum_{j=1}^{M}\dfrac{1}{N_{j}}\sum_{i=1}^{N}\bigg{(}\dfrac{|E_{i,j}|}{W_{i,j}}\bigg{)}^{2}},$
(22)

is less than $1e-4$. In the above equation, $M$ indicates the number of
fields, $j$ is the number of degree of freedom for each field $j=1,2,..M$, $E$
represents the absolute update (say,
$\hbox{$\mathbf{u}$}^{(m)}$-${\hbox{$\mathbf{u}$}}^{(m-1)}$ for displacement,
$m$ being the iteration count), and $W_{i,j}=max(|\mathbf{U}_{i,j}|,S_{j})$.
The entire solution vector is represented with $\mathbf{U}$ and $S_{j}$ refers
to scaling of solution variables222Scaling of solution variables prevents
possible ill-conditioning of the stiffness matrix..

###  5.1 Brittle: Single Edge Notched specimen under Tension (SENT)

A unit square (in mm) embedded with a horizontal notch, midway along the
height is considered, as shown in Figure 3. The length of the notch (shown in
red) is equal to half of the edge length of the plate. The notch is modelled
explicitly in the finite element mesh. A quasi-static loading is applied at
the top boundary in the form of prescribed displacement increment
$\Delta\hbox{$u$}=1e-5$[mm] for the first $450$ steps, following which it is
changed to $1e-6$[mm]. Furthermore, the bottom boundary remains fixed. The
additional model parameters are presented in Table 4.

$\Delta\hbox{$u$}$ Figure 2: SENT experiment

$\Delta\hbox{$u$}$ Figure 3: SENS experiment

Parameters | Value
---|---
Model | Brittle-AT2
$\lambda$ | 121.154 [GPa]
$\mu$ | 80.769 [GPa]
$G_{c}$ | 2700 [N/m]
$l$ | 1.5e-2 [mm]
Element size | $\hbox{$l$}/2$

Table 4: Parameters

$0$$1$$2$$3$$4$$5$$6$$\cdot 10^{-3}$$0$$0.2$$0.4$$0.6$$0.8$$1$Displacement
[mm]Load [kN]LMMMiehe [5]Pen.Ambati [23] (a)

$0$$0.5$$1$$\varphi$ (b)

Figure 4: Figure (a) presents the load-displacement curves for the single edge
notched specimen under tension. Here, LMM and Pen. refer to the Lagrange
Multiplier Method and the Penalty method respectively. Figure (b) shows the
distribution of the phase-field variable at the final step of the analysis.

Figure 4(a) presents the load-displacement curves obtained using the Lagrange
Multiplier Method (LMM) and Penalty (Pen.) method, along with those from the
literature [5, 23]. Both methods yield a similar behaviour compared to [5], in
terms of the peak load and the post-peak behaviour. Moreover, the phase-field
fracture topology at failure in Figure 4(b) is also identical to those
reported in the aforementioned literature.

###  5.2 Brittle: Single Edge Notched specimen under Shear (SENS)

In order to perform a shear test, the SEN specimen is loaded horizontally
along the top edge as shown in Figure 3. The material properties remain same
as presented in Table 4. A quasi-static loading is applied to the top boundary
in the form of prescribed displacement increment $\Delta\hbox{$u$}=1e-4$[mm]
for the first $90$ steps, following which it is changed to $1e-5$[mm].
Furthermore, the bottom boundary remains fixed, and roller support is
implemented in left and right edges thereby restricting the vertical
displacement.

Figure 5(a) presents the load-displacement curves obtained using the Lagrange
Multiplier Method and the Penalty method, along with those from the literature
[5]. The former methods result in $9\%$ higher peak load estimation as
compared to [5]. However, the phase-field fracture topology in the final step
of the analysis is consistent with those reported in the aforementioned
literature.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$\cdot
10^{-2}$$0$$0.2$$0.4$$0.6$$0.8$$1$Displacement [mm]Load [kN]LMMMiehe [5]Pen.
(a)

$0$$0.5$$1$$\varphi$ (b)

Figure 5: Figure (a) presents the load-displacement curves for the single edge
notched specimen under shear. Here, LMM and Pen. refer to the Lagrange
Multiplier Method and the Penalty method respectively. Figure (b) shows the
distribution of the phase-field variable at the final step of the analysis.

###  5.3 Brittle: Notched concrete specimen with a hole

A notched concrete specimen with a hole, shown in Figure 7 is considered in
this study. The experimental and numerical analysis of the same has been
carried out earlier in [23]. The specimen has dimensions 65 x 120 [mm2], the
hole being 20 [mm] in diameter located at (36.5 [mm], 51[mm]). Moreover, a
notch, 10 [mm] in length is located at 65 [mm] height from the bottom of the
plate. The model parameters used in the simulation is presented in Table 5,
while the experimentally observed fracture pattern is shown in Figure 7. As
shown in Figure 7, the plate is loaded via the upper pin (in grey) with
displacement increment $\Delta\hbox{$u$}=1e-3$ mm, while the lower pin (in
grey) is remains fixed.

$\Delta\hbox{$u$}$ Figure 6: Geometry and constraints

Figure 7: Experimental results [23]

Symbol | Value
---|---
Model | Brittle-AT2
$\lambda$ |  1.94 [kN/mm2]
$\mu$ |  2.45 [kN/mm2]
$G_{c}$ | 2.28e-3 [kN/mm]
$l$ | 0.25 [mm]
Element size | $\hbox{$l$}/2$

Table 5: Material properties [48]

Figure 8(a) presents the load-displacement curves obtained using the Lagrange
Multiplier Method and the Penalty method, along with those from the literature
[23, 48]. The curves corresponding to the Lagrange Multiplier Method and the
Penalty method coincide, and predict peak loads close to those obtained in
[23] and [48]. However, the load-displacement curves are significantly
different from [23]. The reason for this behaviour is not investigated in this
manuscript. However, it is important to note that in this work, the upper and
lower pins are assumed to be rigid connectors. As such, the loading or fixed
boundary conditions is applied w.r.t. the centre of the pin. Furthermore, [23]
opted for $\hbox{$l$}=0.1$ [mm], in this manuscript as well as in [48],
$\hbox{$l$}=0.25$ [mm]. Finally, from Figure 8(b), it observed that the final
phase-field fracture topology does match the experimentally observed fracture
patterns in Figure 7.

$0$$0.5$$1$$1.5$$2$$0$$0.2$$0.4$$0.6$$0.8$$1$Displacement [mm]Load
[kN]LMMNavidtehrani [48]Pen.Ambati [23] (a)

$0$$0.5$$1$$\varphi$ (b)

Figure 8: Figure (a) presents the load-displacement curves for the notched
specimen with hole test, obtained from the current implementation, [23], and
[48]. Here, LMM and Pen. refer to the Lagrange Multiplier Method and the
Penalty method respectively. Figure (b) shows the distribution of the phase-
field variable at the final step of the analysis in a section of the specimen.

###  5.4 Quasi-brittle: Concrete three-point bending

A three-point bending experiment on a notched concrete beam reported in [46]
is considered here. The beam has dimensions $450\times 100$ [mm2], and has a
notch $5\times 50$ [mm2]. A schematic of the beam along with the loading
conditions is presented in Figure 9. Displacement-based load increments of
$\Delta\hbox{$u$}=1e-3$ [mm] is enforced throughout the simulation. The model
parameters are presented in Table 6.

$\Delta\hbox{$u$}$ Figure 9: Three point bending test

Parameters | Value
---|---
Model | Quasi-brittle
Softening |  Cornellisen et. al. [43]
$E_{0}$ | 2e4 [MPa]
$\nu$ | 0.2 [-]
$f_{t}$ | 2.4 [MPa]
$G_{c}$ | 0.113 [N/mm]
$l$ | 2.5 [mm]
Element size | $\hbox{$l$}/5$

Table 6: Parameters for three point bending test

Figure 10(a) presents the load-displacement curves obtained using the Lagrange
Multiplier Method and the Penalty method, along with those obtained from the
experiments [46]. It is observed that both methods yield identical curves.
These curves are close to the experimental upper bound (shaded region).
Furthermore, the phase-field fracture topology at the final step in the
analysis in Figure 10(b) is also identical to that reported in [46].

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.5$$1$$1.5$Displacement [mm]Load [kN]LMMPen.
(a) Load-displacement plot

$0$$0.5$$1$$\varphi$ (b)

Figure 10: Figure (a) presents the load-displacement curves for the concrete
three-point bending test. Here, LMM and Pen. refer to the Lagrange Multiplier
Method and the Penalty method respectively. The experimental range is
represented by the shaded area. Figure (b) shows the distribution of the
phase-field variable at the final step of the analysis in a section of the
beam.

###  5.5 Quasi-brittle: Winkler L-panel

The L-shaped panel studied by [47, 49] is considered in this sub-section.
Figure 11 shows the geometry as well as the loading conditions. The longer
edges of the panel are 500 [mm] and the smaller edges are 250 [mm]. The
loading is applied on the edge marked in blue, 30 [mm] in length, and is in
the form of displacement-based increments of $\Delta\hbox{$u$}=1e-3$ [mm]. The
additional model parameters are presented in Table 7.

$\Delta\hbox{$u$}$ Figure 11: Winkler L-panel

Parameters | Value
---|---
Model | Quasi-brittle
Softening |  Cornellisen et. al. [43]
$E_{0}$ | 2.0e4 [MPa]
$\nu$ | 0.18 [-]
$f_{t}$ | 2.5 [MPa]
$G_{c}$ | 0.130 [N/mm]
$l$ | 5 [mm]
Element size | $\hbox{$l$}/5$

Table 7: Parameters for L-shaped panel test [49]

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$2$$4$$6$$8$Displacement [mm]Load [kN]LMMPen. (a)
Load-displacement plot

$0$$0.5$$1$$\varphi$ (b)

Figure 12: Figure (a) presents the load-displacement curves for the concrete
Winkler L-panel test. Here, LMM and Pen. refer to the Lagrange Multiplier
Method and the Penalty method respectively. The experimental range is
represented by the shaded area. Figure (b) shows the distribution of the
phase-field variable at the final step of the analysis in a section of the
beam.

Figure 12(a) presents the load-displacement curves obtained using the Lagrange
Multiplier Method and the Penalty method, along with those obtained from the
experiments [47, 49]. It is observed that the curves obtained using both
methods are identical. These curves are also close to the experimental range
(shaded region). Furthermore, the phase-field fracture topology at the final
step in the analysis in Figure 12(b) is also identical to that reported in
[47, 49].

##  6 Concluding Remarks

An alternative treatment of the phase-field fracture irreversibility
constraint is presented in this manuscript. The fracture irreversibility
constraint $h(\hbox{$\varphi$})\geq 0$ is transformed into an equivalent
equality constraint. This is carried out upon introducing a slack variable
$\theta$, and equating $h(\hbox{$\varphi$})$ to $\theta^{2}$. By definition,
$\theta^{2}$ admits zero or positive values, thereby ensuring the fracture
irreversibility constraint. The constraint is then introduced in the phase-
field fracture energy functional using the Lagrange Multiplier Method and the
Penalty method. The slack variable approach preserves the variational nature
of the phase-field fracture problem, unlike the history variable approach in
[33]. Furthermore, numerical experiments are conducted on benchmark brittle
(AT2) and quasi-brittle fracture problems to demonstrate the efficacy of the
proposed methods.

The scope for future studies include the extension towards complex multi-
physics problems, large-scale problems, that maybe require special
preconditioning techniques. Other future research directions may be towards
accelerating the convergence rate using Anderson acceleration [50] or Large
time increment method [51, 52].

##  7 Software Implementation and Data Availability

The numerical study in Section 5 is carried out using the equation-based
modelling approach in the software package COMSOL Multiphysics 5.6. The source
files would be made available in Github repository of the corresponding author
(https://github.com/rbharali/PhaseFieldFractureCOMSOL/).

## Acknowledgements

The financial support from the Swedish Research Council for Sustainable
Development (FORMAS) under Grant 2018-01249 and the Swedish Research Council
(VR) under Grant 2017-05192 is gratefully acknowledged. The authors would also
like to thank Prof. Laura de Lorenzis (ETH Zürich) for insightful comments and
suggestions.

## References

  * [1] G.A. Francfort and J.-J. Marigo “Revisiting brittle fracture as an energy minimization problem” In _Journal of the Mechanics and Physics of Solids_ 46.8, 1998, pp. 1319–1342 DOI: 10.1016/S0022-5096(98)00034-9
  * [2] B. Bourdin, G.A. Francfort and J-J. Marigo “Numerical experiments in revisited brittle fracture” In _Journal of the Mechanics and Physics of Solids_ 48.4, 2000, pp. 797–826 DOI: https://doi.org/10.1016/S0022-5096(99)00028-9
  * [3] David Mumford and Jayant Shah “Optimal approximations by piecewise smooth functions and associated variational problems” In _Communications on Pure and Applied Mathematics_ 42.5, 1989, pp. 577–685 DOI: 10.1002/cpa.3160420503
  * [4] T. Gerasimov and L. De Lorenzis “On penalization in variational phase-field models of brittle fracture” In _Computer Methods in Applied Mechanics and Engineering_ 354, 2019, pp. 990–1026 DOI: https://doi.org/10.1016/j.cma.2019.05.038
  * [5] C. Miehe, F. Welschinger and M. Hofacker “Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations” In _International Journal for Numerical Methods in Engineering_ 83.10, 2010, pp. 1273–1311 DOI: 10.1002/nme.2861
  * [6] Christian Miehe, Lisa Marie Schänzel and Heike Ulmer “Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids” In _Computer Methods in Applied Mechanics and Engineering_ 294, 2015, pp. 449–485 DOI: 10.1016/j.cma.2014.11.016
  * [7] C. Miehe, M. Hofacker, L.. Schänzel and F. Aldakheel “Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids” In _Computer Methods in Applied Mechanics and Engineering_ 294, 2015, pp. 486–522 DOI: 10.1016/j.cma.2014.11.017
  * [8] M. Ambati, T. Gerasimov and L. De Lorenzis “Phase-field modeling of ductile fracture” In _Computational Mechanics_ 55.5, 2015, pp. 1017–1040 DOI: 10.1007/s00466-015-1151-4
  * [9] A. Mesgarnejad, B. Bourdin and M.M. Khonsari “A variational approach to the fracture of brittle thin films subject to out-of-plane loading” In _Journal of the Mechanics and Physics of Solids_ 61.11, 2013, pp. 2360–2379 DOI: https://doi.org/10.1016/j.jmps.2013.05.001
  * [10] Thanh Tung Nguyen, Julien Réthoré and Marie-Christine Baietto “Phase field modelling of anisotropic crack propagation” In _European Journal of Mechanics - A/Solids_ 65, 2017, pp. 279–288 DOI: https://doi.org/10.1016/j.euromechsol.2017.05.002
  * [11] Christian Miehe and Steffen Mauthe “Phase field modeling of fracture in multi-physics problems. Part III. Crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media” In _Computer Methods in Applied Mechanics and Engineering_ 304 Elsevier B.V., 2016, pp. 619–655 DOI: 10.1016/j.cma.2015.09.021
  * [12] T. Cajuhi, L. Sanavia and L. De Lorenzis “Phase-field modeling of fracture in variably saturated porous media” In _Computational Mechanics_ 61.3 Springer Berlin Heidelberg, 2018, pp. 299–318 DOI: 10.1007/s00466-017-1459-3
  * [13] A. Mikelić, M.. Wheeler and T. Wick “Phase-field modeling through iterative splitting of hydraulic fractures in a poroelastic medium” In _GEM - International Journal on Geomathematics_ 10 Springer International Publishing, 2019 DOI: 10.1007/s13137-019-0113-y
  * [14] Emilio Martínez-Pañeda, Alireza Golahmar and Christian F. Niordson “A phase field formulation for hydrogen assisted cracking” In _Computer Methods in Applied Mechanics and Engineering_ 342, 2018, pp. 742–761 DOI: 10.1016/j.cma.2018.07.021
  * [15] Chuanjie Cui, Rujin Ma and Emilio Martínez-Pañeda “A phase field formulation for dissolution-driven stress corrosion cracking” In _Journal of the Mechanics and Physics of Solids_ 147, 2021, pp. 104254 DOI: https://doi.org/10.1016/j.jmps.2020.104254
  * [16] Marlini Simoes and Emilio Martínez-Pañeda “Phase field modelling of fracture and fatigue in Shape Memory Alloys” In _Computer Methods in Applied Mechanics and Engineering_ 373, 2021, pp. 113504 DOI: https://doi.org/10.1016/j.cma.2020.113504
  * [17] G. Kikis, M. Ambati, L. De Lorenzis and S. Klinkel “Phase-field model of brittle fracture in Reissner–Mindlin plates and shells” In _Computer Methods in Applied Mechanics and Engineering_ 373, 2021, pp. 113490 DOI: https://doi.org/10.1016/j.cma.2020.113490
  * [18] R.U. Patil, B.K. Mishra and I.V. Singh “An adaptive multiscale phase field method for brittle fracture” In _Computer Methods in Applied Mechanics and Engineering_ 329, 2018, pp. 254–288 DOI: 10.1016/j.cma.2017.09.021
  * [19] R.U. Patil, B.K. Mishra, I.V. Singh and T.Q. Bui “A new multiscale phase field method to simulate failure in composites” In _Advances in Engineering Software_ 126, 2018, pp. 9–33 DOI: https://doi.org/10.1016/j.advengsoft.2018.08.010
  * [20] R.U. Patil, B.K. Mishra and I.V. Singh “A multiscale framework based on phase field method and XFEM to simulate fracture in highly heterogeneous materials” In _Theoretical and Applied Fracture Mechanics_ 100, 2019, pp. 390–415 DOI: https://doi.org/10.1016/j.tafmec.2019.02.002
  * [21] F. Fantoni, A. Bacigalupo, M. Paggi and J. Reinoso “A phase field approach for damage propagation in periodic microstructured materials” In _International Journal of Fracture_ , 2019 DOI: 10.1007/s10704-019-00400-x
  * [22] Ritukesh Bharali, Fredrik Larsson and Ralf Jänicke “Computational homogenisation of phase-field fracture” In _European Journal of Mechanics - A/Solids_ 88, 2021, pp. 104247 DOI: https://doi.org/10.1016/j.euromechsol.2021.104247
  * [23] Marreddy Ambati, Tymofiy Gerasimov and Laura De Lorenzis “A review on phase-field models of brittle fracture and a new fast hybrid formulation” In _Computational Mechanics_ 55.2, 2014, pp. 383–405 DOI: 10.1007/s00466-014-1109-y
  * [24] Jian-Ying Wu et al. “Chapter One - Phase-field modeling of fracture” 53, Advances in Applied Mechanics Elsevier, 2020, pp. 1–183 DOI: https://doi.org/10.1016/bs.aams.2019.08.001
  * [25] Laura De Lorenzis and Tymofiy Gerasimov “Numerical implementation of phase-field models of brittle fracture” In _Modeling in Engineering Using Innovative Numerical Methods for Solids and Fluids_ Springer, 2020, pp. 75–101
  * [26] Timo Heister, Mary F. Wheeler and Thomas Wick “A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach” In _Computer Methods in Applied Mechanics and Engineering_ , 2015 DOI: 10.1016/j.cma.2015.03.009
  * [27] T. Gerasimov and L. De Lorenzis “A line search assisted monolithic approach for phase-field computing of brittle fracture” In _Computer Methods in Applied Mechanics and Engineering_ 312, 2016, pp. 276–303 DOI: 10.1016/j.cma.2015.12.017
  * [28] Thomas Wick “An Error-Oriented Newton/Inexact Augmented Lagrangian Approach for Fully Monolithic Phase-Field Fracture Propagation” In _SIAM Journal on Scientific Computing_ 39.4, 2017, pp. B589–B617 DOI: 10.1137/16m1063873
  * [29] Thomas Wick “Modified Newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation” In _Computer Methods in Applied Mechanics and Engineering_ , 2017 DOI: 10.1016/j.cma.2017.07.026
  * [30] Alena Kopaničáková and Rolf Krause “A recursive multilevel trust region method with application to fully monolithic phase-field models of brittle fracture” In _Computer Methods in Applied Mechanics and Engineering_ 360, 2020, pp. 112720 DOI: https://doi.org/10.1016/j.cma.2019.112720
  * [31] B. Bourdin “Numerical implementation of the variational formulation for quasi-static brittle fracture” In _Interfaces and Free Boundaries_ 9, 2007, pp. 411–430 DOI: 10.4171/IFB/171
  * [32] Siobhan Burke, Christoph Ortner and Endre Süli “An Adaptive Finite Element Approximation of a Variational Model of Brittle Fracture” In _SIAM Journal on Numerical Analysis_ 48.3, 2010, pp. 980–1012 DOI: 10.1137/080741033
  * [33] Christian Miehe, Martina Hofacker and Fabian Welschinger “A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits” In _Computer Methods in Applied Mechanics and Engineering_ 199.45-48 Elsevier B.V., 2010, pp. 2765–2778 DOI: 10.1016/j.cma.2010.04.011
  * [34] FA Valentine “The problem of Lagrange with differentiable inequality as added side conditions. 407–448” In _Contributions to the Calculus of Variations_ 37, 1933
  * [35] D P Bertsekas “Nonlinear Programming” In _Journal of the Operational Research Society_ 48.3 Taylor & Francis, 1997, pp. 334–334 DOI: 10.1057/palgrave.jors.2600425
  * [36] Siobhan Burke, Christoph Ortner and Endre Süli “An adaptive finite element approximation of a variational model of brittle fracture” In _SIAM Journal on Numerical Analysis_ 48.3 SIAM, 2010, pp. 980–1012
  * [37] Marco Artina, Massimo Fornasier, Stefano Micheletti and Simona Perotto “Anisotropic mesh adaptation for crack detection in brittle materials” In _SIAM Journal on Scientific Computing_ 37.4 SIAM, 2015, pp. B633–B659
  * [38] Thomas Wick “Goal functional evaluations for phase-field fracture using PU-based DWR mesh adaptivity” In _Computational Mechanics_ 57.6 Springer, 2016, pp. 1017–1035
  * [39] D.. Kelly, J.. De S.., O.. Zienkiewicz and I. Babuska “A posteriori error analysis and adaptive processes in the finite element method: Part I—error analysis” In _International Journal for Numerical Methods in Engineering_ 19.11, 1983, pp. 1593–1619 DOI: https://doi.org/10.1002/nme.1620191103
  * [40] S. Nagaraja et al. “Phase-field modeling of brittle fracture with multi-level hp-FEM and the finite cell method” In _Computational Mechanics_ 63.6, 2019, pp. 1283–1300 DOI: 10.1007/s00466-018-1649-7
  * [41] Kim Pham, Hanen Amor, Jean-Jacques Marigo and Corrado Maurini “Gradient damage models and their use to approximate brittle fracture” In _International Journal of Damage Mechanics_ 20.4 SAGE Publications Sage UK: London, England, 2011, pp. 618–652
  * [42] Jian-Ying Wu “A unified phase-field theory for the mechanics of damage and quasi-brittle failure” In _Journal of the Mechanics and Physics of Solids_ 103 Elsevier, 2017, pp. 72–99
  * [43] H Cornelissen, D Hordijk and H Reinhardt “Experimental determination of crack softening characteristics of normalweight and lightweight” In _Heron_ 31.2, 1986, pp. 45–46
  * [44] M. Hintermüller, K. Ito and K. Kunisch “The Primal-Dual Active Set Strategy as a Semismooth Newton Method” In _SIAM Journal on Optimization_ 13.3, 2002, pp. 865–888 DOI: 10.1137/S1052623401383558
  * [45] R.. Tapia “A stable approach to Newton’s method for general mathematical programming problems inRn” In _Journal of Optimization Theory and Applications_ 14.5, 1974, pp. 453–476 DOI: 10.1007/BF00932842
  * [46] Jan Rots “Computational modeling of concrete fracture”, 1988
  * [47] Bernhard Josef Winkler “Traglastuntersuchungen von unbewehrten und bewehrten Betonstrukturen auf der Grundlage eines objektiven Werkstoffgesetzes für Beton” Innsbruck University Press, 2001
  * [48] Yousef Navidtehrani, Covadonga Betegón and Emilio Martínez-Pañeda “A Unified Abaqus Implementation of the Phase Field Fracture Method Using Only a User Material Subroutine” In _Materials_ 14.8, 2021 DOI: 10.3390/ma14081913
  * [49] Jörg F Unger, Stefan Eckardt and Carsten Könke “Modelling of cohesive crack growth in concrete structures with the extended finite element method” In _Computer methods in applied mechanics and engineering_ 196.41-44 Elsevier, 2007, pp. 4087–4100
  * [50] Donald G Anderson “Iterative procedures for nonlinear integral equations” In _Journal of the ACM (JACM)_ 12.4 ACM New York, NY, USA, 1965, pp. 547–560
  * [51] Dan Givoli, Ritukesh Bharali and Lambertus J Sluys “LATIN: A new view and an extension to wave propagation in nonlinear media” In _International Journal for Numerical Methods in Engineering_ 112.2 Wiley Online Library, 2017, pp. 125–156
  * [52] Ritukesh Bharali “An Adaptive Generic Incremental LATIN method”, 2017


# Mie scattering theory: A review of physical features and limitations

Yuriy A. Akimov<EMAIL_ADDRESS>Institute of High Performance
Computing (IHPC), Agency for Science, Technology and Research (A*STAR),
1 Fusionopolis Way, #16-16 Connexis, Singapore 138632

###### Abstract

Mie theory is the classical problem for modeling of light scattering by
spherical particles. In this paper, we perform a spherical harmonic analysis
of its solution for the induced fields to reveal the physics underlying the
resonant behavior predicted for scattering and absorption. We disclose the two
distinct groups of current-sourced and current-free scattered fields, whose
interference makes light-matter interaction resonant in Mie theory for every
orbital index. Being a model, the current-free scattered fields naturally
limit applicability of the classical solution. We discuss those limitations
and demonstrate the ways for further refinement of the theory for description
of the excitation source, sphere interface and scatterer localization.

##  I Introduction

Mie theory is the classical problem for light scattering on a spherical
particle [1, 2, 3]. It is named after German physicist Gustav Mie, who first
published his paper on scattering of an electromagnetic wave on a sphere in
1908 [4]. Since then, this theory played a key part in shaping of our
knowledge about light-matter interaction. It bridged the theories of Rayleigh
scattering for small particles and Rayleigh-Gans-Debye scattering for larger
ones [5], enabling calculations for particles of arbitrary size. It explained
frequency-selective scattering and size-dependent color change, which occur in
solutions of colloid nanoparticles. Later, the theory was successfully
extended to multilayer spheres [6] and non-plane incident fields emitted by
dipoles [7] or swift electrons [8]. Now, Mie theory is extensively used for
design of optical resonances in plasmonic and photonic structures at the
nanoscale [9, 10, 11, 12].

Despite the long history of success, Mie theory keeps a number of fundamental
questions unanswered. For instance, why Mie scattering and absorption behave
differently with change of particle size? Why does resonant scattering
monotonically increase with particle size, while resonant absorption does not
[13]? Mie theory explains these behaviors by resonant nature of the Mie
scattering coefficients. But how do they become resonant? What is the physics
underlying their resonances? If the resonances in Mie coefficients arise from
interference of the dipole and higher-order currents induced in the sphere as
was suggested in Refs. [14, 15, 16], then how do similar resonances appear in
deeply subwavelength metal particles [17, 18] whose higher-order currents are
negligible? How can a dipole current along make the Mie scattering
coefficients resonant? Other questions arise from the way the light-matter
interaction is treated in Mie theory. The theory successfully predicts static
flows for the electromagnetic energy, but not for the momentum. This suggests
that Mie theory provides a truncated description for light-matter interaction.
Then, how limited is that description? All these questions highlight lack of
our understanding of the Mie scattering. To address them, we perform a
rigorous spherical harmonic analysis of Mie theory revealing the physics
underlying its solution for scattering and absorption of light.

##  II Mie theory

Mie theory describes scattering and absorption of harmonic plane waves of
angular frequency $\omega$ incident on a spherical particle. Following
Maxwell’s equations, electromagnetic fields in domains with uniform dielectric
permittivity $\varepsilon$ in the absence of any external charges and currents
can be completely described with the transverse magnetic (TM) and transverse
electric (TE) fields as follows [9, 10]:

$\displaystyle\displaystyle{\bf H}=\left({\bf H}^{\rm TM}-\frac{{\rm
i}}{k_{0}}\sqrt{\frac{\varepsilon_{0}}{\mu_{0}}}\nabla\times{\bf E}^{\rm
TE}\right)\mathop{\rm e}\nolimits^{-{\rm i}\omega t},$ (1)
$\displaystyle\displaystyle{\bf E}=\left({\bf E}^{\rm TE}+\frac{{\rm
i}}{k_{0}\varepsilon}\sqrt{\frac{\mu_{0}}{\varepsilon_{0}}}\nabla\times{\bf
H}^{\rm TM}\right)\mathop{\rm e}\nolimits^{-{\rm i}\omega t},$ (2)

where ${\bf H}^{\rm TM}$ and ${\bf E}^{\rm TE}$ are the amplitudes of the
governing fields for the TM and TE polarizations,
$k_{0}=\omega\sqrt{\varepsilon_{0}\mu_{0}}$ is the vacuum wavenumber,
$\varepsilon_{0}$ and $\mu_{0}$ are the electric and magnetic constants. In
the spherical geometry, the governing TM and TE fields can be decomposed over
the vector spherical harmonics of different orbital and azimuthal indices
$l,m$ to separate their radial and angular dependencies [9, 10]:

$\displaystyle\displaystyle{\bf H}^{\rm
TM}=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}{H}_{lm}\left({\bf
e}_{\phi}\frac{{\rm\partial}Y_{lm}}{{\rm\partial}\theta}-\frac{{\bf
e}_{\theta}}{\sin\theta}\frac{{\rm\partial}Y_{lm}}{{\rm\partial}\phi}\right),$
(3) $\displaystyle\displaystyle{\bf E}^{\rm
TE}=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}{E}_{lm}\left({\bf
e}_{\phi}\frac{{\rm\partial}Y_{lm}}{{\rm\partial}\theta}-\frac{{\bf
e}_{\theta}}{\sin\theta}\frac{{\rm\partial}Y_{lm}}{{\rm\partial}\phi}\right),$
(4)

where $(r,\theta,\phi)$ are the spherical coordinates, and ${\bf e}_{r},{\bf
e}_{\theta},{\bf e}_{\phi}$ are the respective unit vectors. In this
separation, the angular dependence is fully given by the scalar spherical
harmonics

$Y_{lm}(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}\,P_{l}^{m}(\cos\theta)\mathop{\rm
e}\nolimits^{{\rm i}m\phi},$ (5)

where $P_{l}^{m}(\cos\theta)$ are the associated Legendre polynomials, while
the radial dependence is given by the scalar functions ${H}_{lm}$ and
${E}_{lm}$ generally independent of each other and obeying the wave equations

$\displaystyle\displaystyle\frac{{\rm d}^{2}H_{lm}}{{\rm
d}r^{2}}+\frac{2}{r}\frac{{\rm d}H_{lm}}{{\rm
d}r}-\left[\frac{l(l+1)}{r^{2}}-k_{0}^{2}\varepsilon\right]H_{lm}=0,$ (6)
$\displaystyle\displaystyle\frac{{\rm d}^{2}E_{lm}}{{\rm
d}r^{2}}+\frac{2}{r}\frac{{\rm d}E_{lm}}{{\rm
d}r}-\left[\frac{l(l+1)}{r^{2}}-k_{0}^{2}\varepsilon\right]E_{lm}=0.$ (7)

For a spherical particle of radius $R$, Mie theory defines the radial
functions in a piecewise manner [9]:

$\displaystyle\displaystyle{E}_{lm}=\left\\{\begin{array}[]{ll}E_{lm}^{\rm
int},&r<R,\\\\[3.0pt] E_{lm}^{\rm inc}+E_{lm}^{\rm
sca},&r>R,\end{array}\right.$ (10)
$\displaystyle\displaystyle{H}_{lm}=\left\\{\begin{array}[]{ll}H_{lm}^{\rm
int},&r<R,\\\\[3.0pt] H_{lm}^{\rm inc}+H_{lm}^{\rm
sca},&r>R,\end{array}\right.$ (13)

where the incident, scattered and internal fields are given by

$\displaystyle\displaystyle H_{lm}^{\rm inc}=\widetilde{H}_{lm}j_{l}(k_{e}r),$
(14) $\displaystyle\displaystyle E_{lm}^{\rm
inc}=\widetilde{E}_{lm}j_{l}(k_{e}r),$ (15) $\displaystyle\displaystyle
H_{lm}^{\rm sca}=-a_{l}\widetilde{H}_{lm}h_{l}^{(1)}(k_{e}r),$ (16)
$\displaystyle\displaystyle E_{lm}^{\rm
sca}=-b_{l}\widetilde{E}_{lm}h_{l}^{(1)}(k_{e}r),$ (17)
$\displaystyle\displaystyle H_{lm}^{\rm
int}=d_{l}\widetilde{H}_{lm}j_{l}(k_{i}r),$ (18) $\displaystyle\displaystyle
E_{lm}^{\rm int}=c_{l}\widetilde{E}_{lm}j_{l}(k_{i}r)$ (19)

with $k_{i,e}=k_{0}\varepsilon_{i,e}^{1/2}$ being the wavenumbers of
transverse fields in the internal ($r<R$) and external ($r>R$) homogeneous
domains denoted with the indices $i$ and $e$, accordingly. The radial
distributions of the scattered and internal fields in Mie theory are chosen to
be the spherical Hankel functions of the first kind $h_{l}^{(1)}(k_{e}r)$ and
the spherical Bessel functions $j_{l}(k_{i}r)$ providing finite solutions to
Eqs. (6) and (7) in the two domains, while the distributions of the incident
fields are taken from the spherical decomposition of the $x$-polarized
transverse electromagnetic (TEM) plane wave propagating in the $z$-direction:

$\displaystyle\displaystyle{{\bf E}}^{\rm inc}={{\bf e}}_{x}E_{0}\mathop{\rm
e}\nolimits^{{\rm i}(k_{e}r\cos\theta-\omega t)},$ (20)
$\displaystyle\displaystyle{{\bf H}}^{\rm
inc}=\sqrt{\frac{\varepsilon_{0}\varepsilon_{e}}{\mu_{0}}}{{\bf
e}}_{y}E_{0}\mathop{\rm e}\nolimits^{{\rm i}(k_{e}r\cos\theta-\omega t)},$
(21)

where $E_{0}$ is the complex amplitude of the incident electric field. If the
external medium is fully transparent with ${\rm Im}~{}\varepsilon_{e}=0$, then
the decomposition results in the following spherical harmonics amplitudes for
the incident fields [9]:

$\displaystyle\displaystyle\widetilde{H}_{lm}=-{\rm
i}^{l}\sqrt{\frac{\pi(2l+1)}{l(l+1)}}\sqrt{\frac{\varepsilon_{0}\varepsilon_{e}}{\mu_{0}}}m\delta_{m,\pm
1}E_{0},$ (22) $\displaystyle\displaystyle\widetilde{E}_{lm}=-{\rm
i}^{l+1}\sqrt{\frac{\pi(2l+1)}{l(l+1)}}\delta_{m,\pm 1}E_{0}.$ (23)

These amplitudes are used in induced fields (16)–(19) for normalization
together with dimensionless coefficients $a_{l}$, $b_{l}$, $c_{l}$ and
$d_{l}$, known as Mie coefficients. The latter are given by the boundary
conditions implying the continuity of the tangential components of ${\bf
E}^{\rm TE}$ and ${\bf H}^{\rm TM}$ on the sphere interface $r=R$ [1, 2, 3]:

$\displaystyle\displaystyle
a_{l}=\frac{q_{i}\psi_{l}(q_{i})\psi_{l}^{\prime}(q_{e})-q_{e}\psi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}{q_{i}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{e}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(24) $\displaystyle\displaystyle
b_{l}=\frac{q_{e}\psi_{l}(q_{i})\psi_{l}^{\prime}(q_{e})-q_{i}\psi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}{q_{e}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(25) $\displaystyle\displaystyle
c_{l}=\frac{q_{i}\psi_{l}(q_{e})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{e})}{q_{e}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(26) $\displaystyle\displaystyle
d_{l}=\frac{q_{i}\psi_{l}(q_{e})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{e})}{q_{i}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{e}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(27)

where $\psi_{l}(q_{i})=q_{i}j_{l}(q_{i})$ and
$\xi_{l}(q_{e})=q_{e}h_{l}^{(1)}(q_{e})$ are the Riccati–Bessel functions,
with $q_{i}=k_{i}R$ and $q_{e}=k_{e}R$.

With the scattered fields determined in the spherical system, Mie theory
enables integration of the power scattered and absorbed by the particle and
calculates the total cross-sections of scattering and absorption for the
incident plane waves [1, 2, 3]:

$\displaystyle\displaystyle\sigma_{\rm
sca}=\frac{2\pi}{k_{e}^{2}}\sum\limits_{l=1}^{\infty}(2l+1)(|a_{l}|^{2}+|b_{l}|^{2}),$
(28) $\displaystyle\displaystyle\sigma_{\rm
abs}=\frac{2\pi}{k_{e}^{2}}\sum\limits_{l=1}^{\infty}(2l+1)[{\rm
Re}(a_{l}+b_{l})-(|a_{l}|^{2}+|b_{l}|^{2})],~{}~{}~{}$ (29)

where the contributions by the TM and TE spherical vector harmonics are
represented with the scattering coefficients $a_{l}$ and $b_{l}$, accordingly.

##  III Features

###  III.1 Two types of scattered fields

Figure 1: Contributions of source-free scattered fields of different orbital
indices and polarizations to the governing fields ${\bf H}^{\rm TM}$ and ${\bf
E}^{\rm TE}$ under $q_{e}=1$. The contributions are normalized by the
amplitudes of the plane incident fields $E_{0}$ and
$H_{0}=\sqrt{\varepsilon_{0}\varepsilon_{e}/\mu_{0}}E_{0}$.

As Mie theory is linear to electromagnetic fields, one can split the incident
fields into two parts:

$\displaystyle\displaystyle H_{lm}^{\rm inc}=H_{lm}^{\rm inc,1}+H_{lm}^{\rm
inc,2},$ (30) $\displaystyle\displaystyle E_{lm}^{\rm inc}=E_{lm}^{\rm
inc,1}+E_{lm}^{\rm inc,2},$ (31)

and study how each part separately induces the scattered and internal fields:

$\displaystyle\displaystyle H_{lm}^{\rm sca}=H_{lm}^{\rm sca,1}+H_{lm}^{\rm
sca,2},$ (32) $\displaystyle\displaystyle E_{lm}^{\rm sca}=E_{lm}^{\rm
sca,1}+E_{lm}^{\rm sca,2},$ (33) $\displaystyle\displaystyle H_{lm}^{\rm
int}=H_{lm}^{\rm int,1}+H_{lm}^{\rm int,2},$ (34) $\displaystyle\displaystyle
E_{lm}^{\rm int}=E_{lm}^{\rm int,1}+E_{lm}^{\rm int,2}.$ (35)

Particularly, one can decompose the incident fields into spherically outgoing
and spherically incoming fields:

$\displaystyle\displaystyle H_{lm}^{{\rm
inc},j}=\frac{1}{2}\widetilde{H}_{lm}h_{l}^{(j)}(k_{e}r),$ (36)
$\displaystyle\displaystyle E_{lm}^{{\rm
inc},j}=\frac{1}{2}\widetilde{E}_{lm}h_{l}^{(j)}(k_{e}r),$ (37)

where $h_{l}^{(j)}(k_{e}r)$ are the spherical Hankel functions of the kind
$j=1,2$. Then, the respective internal and scattered fields can be written as
follows:

$\displaystyle\displaystyle H_{lm}^{{\rm
sca},j}=-a_{lm}^{(j)}\widetilde{H}_{lm}h_{l}^{(1)}(k_{e}r),$ (38)
$\displaystyle\displaystyle E_{lm}^{{\rm
sca},j}=-b_{lm}^{(j)}\widetilde{E}_{lm}h_{l}^{(1)}(k_{e}r),$ (39)
$\displaystyle\displaystyle H_{lm}^{{\rm
int},j}=d_{lm}^{(j)}\widetilde{H}_{lm}j_{l}(k_{i}r),$ (40)
$\displaystyle\displaystyle E_{lm}^{{\rm
int},j}=c_{lm}^{(j)}\widetilde{E}_{lm}j_{l}(k_{i}r).$ (41)

The new Mie coefficients appear to be

$\displaystyle\displaystyle a_{l}^{(1)}=b_{l}^{(1)}=\frac{1}{2},\quad
c_{l}^{(1)}=d_{l}^{(1)}=0$ (42)

for the spherically outgoing incident fields and

$\displaystyle\displaystyle
a_{l}^{(2)}=\frac{1}{2}\frac{q_{i}\psi_{l}(q_{i})\zeta_{l}^{\prime}(q_{e})-q_{e}\zeta_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}{q_{i}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{e}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(43) $\displaystyle\displaystyle
b_{l}^{(2)}=\frac{1}{2}\frac{q_{e}\psi_{l}(q_{i})\zeta_{l}^{\prime}(q_{e})-q_{i}\zeta_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}{q_{e}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(44) $\displaystyle\displaystyle
c_{l}^{(2)}=\frac{1}{2}\frac{q_{i}\psi_{l}(q_{e})\zeta_{l}^{\prime}(q_{e})-q_{i}\zeta_{l}(q_{e})\psi_{l}^{\prime}(q_{e})}{q_{e}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(45) $\displaystyle\displaystyle
d_{l}^{(2)}=\frac{1}{2}\frac{q_{i}\zeta_{l}(q_{e})\xi_{l}^{\prime}(q_{e})-q_{i}\xi_{l}(q_{e})\zeta_{l}^{\prime}(q_{e})}{q_{i}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{e}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}$
(46)

for the spherically incoming incident fields, where
$\zeta_{l}(q_{e})=q_{e}h_{l}^{(2)}(q_{e})$ are the spherically incoming
Riccati-Bessel functions.

One can see the unusual scattering of spherically outgoing fields. Their
scattered fields appear fully compensating the incident fields, $H_{lm}^{\rm
sca,1}=-H_{lm}^{\rm inc,1}$, $E_{lm}^{\rm sca,1}=-E_{lm}^{\rm inc,1}$, and
accompanied by the zero internal fields, $H_{lm}^{\rm int,1}=0$, $E_{lm}^{\rm
int,1}=0$. This type of scattering corresponds to the trivial solution of
Maxwell’s equations with the zero total fields, ${\bf E}=0$, ${\bf H}=0$,
inside and outside of the sphere. Such a solution exhibit nonzero scattered
fields under the zero internal current induced in the scatterer. Thus, the
scattered fields $H_{lm}^{{\rm sca},1}$ and $E_{lm}^{{\rm sca},1}$ are source-
free. In contrast to them, the scattered fields $H_{lm}^{{\rm sca},2}$ and
$E_{lm}^{{\rm sca},2}$ are accompanied by nonzero internal fields,
$H_{lm}^{\rm int,1}\neq 0$, $E_{lm}^{\rm int,1}\neq 0$, highlighting their
current-sourced nature. In other words, Mie theory deals with two
fundamentally different scattered fields: (i) source-free and (ii) current-
sourced ones.

Figure 2: Amplitude ratios of the current-sourced and source-free scattered
fields as functions of ${\rm Im}~{}\varepsilon_{i}/\varepsilon_{e}$ for
different orbital indices and polarizations under $q_{e}=1$. The labels show
the values of ${\rm Re}~{}\varepsilon_{i}/\varepsilon_{e}=-5$, 1 and 5 kept
fixed for the respective curves.

The paradox of two kinds of scattered fields comes from the formulation of Mie
theory. It extends the conventional electrodynamic formulation of current-
sourced fields and introduces free sourceless fields. Namely, the incident
fields used in Mie theory are source-free [1, 2, 3]. As a result, all the
fields in Mie theory, including the scattered ones, are composed of the
current-sourced and source-free electromagnetic fields. The former are excited
by the internal currents induced in the scatterer, while the latter appear as
the trivial response of Maxwell’s equations to the introduced free incident
fields.

Noteworthy that the source-free scattered fields are independent of the sphere
size and material, so the effect of scatterer in Mie theory is solely defined
by the current-sourced fields. The typical distribution of the source-free
scattered fields with different orbital indices $l$ is shown in Fig. 1. The
current-sourced scattered fields have similar distributions, but with
different complex amplitudes. For non-absorbing particles with ${\rm
Im}~{}\varepsilon_{i}=0$, the amplitudes of the current-sourced scattered
fields differ from those of the source-free scattered fields by phase only:
$|a_{l}^{(2)}|=|a_{l}^{(1)}|$ for the TM polarization and
$|b_{l}^{(2)}|=|b_{l}^{(1)}|$ for the TE one. For absorbing particles with
${\rm Im}~{}\varepsilon_{i}>0$, the complex amplitudes of the source-free and
current-sourced scattered fields differ by both amplitudes and phases, as
shown in Fig. 2.

###  III.2 Super-radiating and non-radiating states

Figure 3: Super-radiating (dashed lines) and non-radiating (solid lines)
states of the TM and TE polarizations for the orbital indices $l=1,2,3$.
Figure 4: Normalized scattering cross-section with respective contributions of
the TM and TE fields with $l=1,2$ for (a)–(c) a silver particle of $R=300$ nm
and (d)–(f) a silicon particle of $R=200$ nm embedded in silicon dioxide. The
limits in (b), (c) and (e), (f) are given by $\sigma_{{\rm sca},l}^{\rm sr}$.

As the scattered fields in Mie theory are not fully current-sourced, their
effect on $\sigma_{\rm sca}$ and $\sigma_{\rm abs}$ can be understood from the
interference of the two types of fields:

$\displaystyle\displaystyle a_{l}=a_{l}^{(1)}+a_{l}^{(2)},$ (47)
$\displaystyle\displaystyle b_{l}=b_{l}^{(1)}+b_{l}^{(2)},$ (48)

rather than from the classical analysis of the current induced in the particle
[2, 14] that is obviously inapplicable to the source-free electromagnetic
fields. Following this interference, spherical particles support a number of
super-radiating and non-radiating states for each polarization and orbital
index $l$, as shown in Fig. 3.

According to Eq. (28), the super-radiating states exhibit either $a_{l}=1$ or
$b_{l}=1$. They appear at those size parameters $q_{e}$ and relative
dielectric permittivities $\varepsilon_{i}/\varepsilon_{e}$, which result in
the constructive interference of the source-free and current-sourced scattered
fields: $a_{l}^{(2)}=a_{l}^{(1)}$ for the TM polarization and
$b_{l}^{(2)}=b_{l}^{(1)}$ for the TE one. Contrary, the non-radiating states,
also known as anapoles, exhibit $a_{l}=0$ or $b_{l}=0$ caused by the
destructive interference of the two kinds of scattered fields:
$a_{l}^{(2)}=-a_{l}^{(1)}$ for the TM polarization and
$b_{l}^{(2)}=-b_{l}^{(1)}$ for the TE one.

As the super-radiating and non-radiating states require purely real dielectric
permittivity, those states are not achievable for real materials with finite
frequency dispersion and dissipation. Nonetheless, they define the limits,
when every polarization and orbital index contribute

$\displaystyle\displaystyle\sigma_{l}^{\rm
sca,sr}=\frac{2\pi}{k_{e}^{2}}(2l+1),\quad\sigma_{{\rm abs},l}^{\rm sr}=0,$
(49) $\displaystyle\sigma_{{\rm sca},l}^{\rm nr}=0,\qquad\sigma_{{\rm
abs},l}^{\rm nr}=0$ (50)

to the total cross-sections $\sigma_{\rm sca}$ and $\sigma_{\rm abs}$. Figure
4 shows the typical interference fringes in scattering cross-section of silver
and silicon spheres with the upper limits given by $\sigma_{{\rm sca},l}^{\rm
sr}$.

Figure 5: Normalized scattering cross-section with respective contributions of
the TM and TE fields with $l=1,2$ for (a)–(c) a silver particle of $R=25$ nm
and (d)–(f) a silicon particle of $R=40$ nm embedded in silicon dioxide. The
limits in (b) and (e) are given by $\sigma_{{\rm sca},l}^{\rm sr}$.

Noteworthy that the density of interference fringes observed in scattering
spectrum is heavily affected by the size of spheres. Figure 3 suggests that by
moving towards lower $q_{e}$, we can reduce the number of available supper-
radiating states and, hence, the density of scattering fringes for every $l$.
Furthermore, if we go down to $q_{e}\leq 2$, we can filter out contributions
with high orbital indices, whose current-sourced and source-free scattered
fields can be kept destructively interfering each other. Such cases are shown
in Fig. 5 for silver and silicon nanospheres with the dominant dipolar ($l=1$)
contributions. They demonstrate one TM maximum for the metallic particle and
two TE and one TM maxima for the dielectric one, which arise from the
interference of the dipolar scattered fields that runs in a more constructive
way.

Figure 6: Super-absorbing states of the TM and TE polarizations for the
orbital indices $l=1,2,3$.  Figure 7: Normalized absorption cross-section with
respective contributions of the TM and TE fields with $l=1,2$ for (a)–(c) a
silver particle of $R=300$ nm and (d)–(f) a silicon particle of $R=200$ nm
embedded in silicon dioxide. The limits in (b), (c) and (e), (f) are given by
$\sigma_{{\rm abs},l}^{\rm sa}$.

In addition to the mentioned maxima, Fig. 5 demonstrates several scattering
minima caused by the interference of the dipolar scattered fields that runs in
a more destructive way. The low-frequency minimum common to both the silver
and silicon cases arises from the fundamental non-radiating state given by
$q_{e}=0$. This state is independent of the polarization, orbital index and
scatterer material, as shown in Fig. 3. Therefore, the associated regime,
known as the Rayleigh scattering [5, 3], is common to all structures. It is
observed for $q_{e}\ll 1$ with $|a_{1}|\gg|b_{1}|,|a_{2}|,|b_{2}|,...$, where
$|a_{1}|\propto q_{e}^{3}$. The effect of another fundamental non-radiating
state given by $\varepsilon_{i}=\varepsilon_{e}$ is seen in Fig. 5(a) for the
silver particle around $\hbar\omega=3.9$ eV. This state is independent of the
polarization, orbital index and scatterer size, as shown in Fig. 3, but
requires a real-valued $\varepsilon_{i}$. As a result, the associated minimum
of scattering is common to the structures that exhibit low optical loss [17,
18].

###  III.3 Super-absorbing states

Another limiting case important for understanding of Mie theory is the maximum
achievable absorption. Following Eq. (29), the super-absorbing states must
exhibit either $a_{l}=1/2$ or $b_{l}=1/2$. These states require the complex
values of $\varepsilon_{i}/\varepsilon_{e}$ optimized for every size parameter
$q_{e}$, as shown in Fig. 6. The above conditions also can be considered as
those when the current-sourced scattered fields vanish: $a_{l}^{(2)}=0$ for
the TM polarization or $b_{l}^{(2)}=0$ for the TE one. In other words, the
super-absorbing states appear exactly the source-free scattering states.
Contrary to the super-scattering states, there are multiple TM and TE super-
absorbing states with ${\rm Re~{}}\varepsilon_{i}/\varepsilon_{e}>0$ and only
one TM state with ${\rm Re~{}}\varepsilon_{i}/\varepsilon_{e}<0$.

The super-absorbing states define the upper limits for absorption
contributions

$\displaystyle\sigma_{{\rm abs},l}^{\rm sa}=\frac{\pi}{2k_{e}^{2}}(2l+1)$ (51)

by every polarization and orbital index. At those states, the contribution to
the scattering and absorption cross-section appears exactly the same, equal to
the quarter of the super-radiating limit:

$\displaystyle\sigma_{{\rm sca},l}^{\rm sa}=\sigma_{{\rm abs},l}^{\rm
sa}=\frac{1}{4}\sigma_{{\rm sca},l}^{\rm sr}.$ (52)

Figure 8: Normalized absorption cross-section with respective contributions of
the TM and TE fields with $l=1,2$ for (a)–(c) a silver particle of $R=25$ nm
and (d)–(f) a silicon particle of $R=40$ nm embedded in silicon dioxide. The
limits in (b) and (e) are given by $\sigma_{{\rm abs},l}^{\rm sa}$.

Contrary to the super-radiating states, the super-absorbing ones are
achievable for dispersive materials with finite dissipation under the proper
size and material design. The typical absorption fringes given by the
interference of the source-free and current-sourced scattered fields and
bounded by $\sigma_{{\rm abs},l}^{\rm sa}$ are shown in Fig. 7 for silver and
silicon particles.

Similar to scattering, the fringes in absorption are strongly affected by the
size of spheres. We can reduce the density of fringes for every orbital $l$ by
moving towards lower $q_{e}$, as suggested by Fig. 6. If we go down to
$q_{e}\leq 2$, we can filter out contributions with high orbital indices,
whose current-sourced and source-free scattered fields destructively interfere
each other. Such cases are shown in Fig. 8 for silver and silicon nanospheres
with the dominant dipolar ($l=1$) contributions, where the respective TM and
TE contributions are maximized for the particle absorption. Note that at the
limiting cases given by $\sigma_{{\rm abs},l}^{\rm sa}$ the corresponding
scattered fields are purely source-free.

##  IV Limitations

In our analysis of the solution for induced fields, we already saw the
modeling essence of Mie theory. Below, we will consider the major assumptions
made in the theory for description of the excitation source, sphere interface
and scatterer localization and discuss the ways to relax the limitations.

###  IV.1 Excitation source

Figure 9: Ratios $\eta_{\rm sca}$ of the modified scattering cross-section
contributions given by Eq. (59) to the classical ones for the TM and TE
scattered fields with $l=1,2,3$ under $q_{e}=1$ and
$\varepsilon_{i}/\varepsilon_{e}=5+{\rm i}~{}0.1$.

Use of external incident fields for modeling of source of excitation is a
common approach employed in many scattering problems [3], including Mie
theory. It allows us to exclude the excitation current from consideration and
truncate the problem to a smaller domain. At the same time, it makes the
description of electromagnetic interaction between the excitation current and
the scatterer electrodynamically incomplete. Indeed, incident fields (14) and
(15) can be treated as the far-fields of an external current located at $r\geq
R_{\rm ext}$ far from the scatterer, $R_{\rm ext}\gg R$, whereas the scattered
fields written in the form of Eqs. (16) and (17) do not support existence of
any current at $r>R$. To account for the external excitation current, the
scattered fields must be given by linear combinations of $h^{(1)}(k_{e}r)$ and
$h^{(2)}(k_{e}r)$:

$\displaystyle\displaystyle H_{lm}^{\rm
sca}=-a_{l}\widetilde{H}_{lm}[h_{l}^{(1)}(k_{e}r)+\alpha_{l}h_{l}^{(2)}(k_{e}r)],$
(53) $\displaystyle\displaystyle E_{lm}^{\rm
sca}=-b_{l}\widetilde{E}_{lm}[h_{l}^{(1)}(k_{e}r)+\beta_{l}h_{l}^{(2)}(k_{e}r)]$
(54)

for $R<r<R_{\rm ext}$, where $\alpha_{l}$ and $\beta_{l}$ are the complex
interaction coefficients given by the source currents located at $r\geq R_{\rm
ext}$. With this modification, Mie scattering coefficients change to

$\displaystyle\displaystyle
a_{l}=\frac{q_{i}\psi_{l}(q_{i})\psi_{l}^{\prime}(q_{e})-q_{e}\psi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}{q_{i}\psi_{l}(q_{i})A_{l}^{\prime}(q_{e})-q_{e}A_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(55) $\displaystyle\displaystyle
b_{l}=\frac{q_{e}\psi_{l}(q_{i})\psi_{l}^{\prime}(q_{e})-q_{i}\psi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})}{q_{e}\psi_{l}(q_{i})B_{l}^{\prime}(q_{e})-q_{i}B_{l}(q_{e})\psi_{l}^{\prime}(q_{i})},$
(56)

where

Figure 10: Ratios $\eta_{\rm abs}$ of the modified absorption cross-section
contributions given by Eq. (60) to the classical ones for the TM and TE
scattered fields with $l=1,2,3$ under $q_{e}=1$ and
$\varepsilon_{i}/\varepsilon_{e}=5+{\rm i}~{}0.1$.

$\displaystyle\displaystyle
A_{l}(q_{e})=\xi_{l}(q_{e})+\alpha_{l}\zeta_{l}(q_{e}),$ (57)
$\displaystyle\displaystyle
B_{l}(q_{e})=\xi_{l}(q_{e})+\beta_{l}\zeta_{l}(q_{e}).$ (58)

Finally, the scattering and absorption cross-sections appear

$\displaystyle\displaystyle\sigma_{\rm
sca}=\frac{2\pi}{k_{e}^{2}}\sum\limits_{l=1}^{\infty}(2l+1)[|a_{l}|^{2}(1-|\alpha_{l}|^{2})+|b_{l}|^{2}(1-|\beta_{l}|^{2})],$
(59) $\displaystyle\displaystyle\sigma_{\rm
abs}=\frac{2\pi}{k_{e}^{2}}\sum\limits_{l=1}^{\infty}(2l+1){\rm
Re}[a_{l}(1-\alpha_{l})+b_{l}(1-\beta_{l})-$
$\displaystyle\displaystyle|a_{l}|^{2}(1-|\alpha_{l}|^{2})-|b_{l}|^{2}(1-|\beta_{l}|^{2})].$
(60)

As we can see, the account for external current changes the super-radiating,
non-radiating and super-absorbing states and eventually affects the total
scattering and absorption cross-sections, especially for low orbital indices,
as shown in Figs. 9 and 10.

###  IV.2 Sphere interface

Although the original Mie theory focuses of energy characteristics of the
source-scatterer interaction, it can be used to study momentum characteristics
as well. For instance, we can use Mie theory for calculation of the optical
force that acts on the sphere from the source side:

$\displaystyle\displaystyle{\bf F}^{\rm opt}=\frac{1}{2}{\rm
Re}\int(\rho^{*}{\bf E}^{\rm inc}+\mu_{0}{\bf J}^{*}\times{\bf H}^{\rm
inc}){\rm d}V.$ (61)

Here, $\rho$ is the volume density of the surface charge induced on the sphere
interface:

$\displaystyle\displaystyle\rho=\varepsilon_{0}\nabla\cdot{\bf
E}=\widetilde{\rho}(\theta,\phi)\delta(r-R)\exp(-{\rm i}\omega t)$ (62)

with $\widetilde{\rho}(\theta,\phi)$ being the angular distribution and
$\delta(r-R)$ being the radial Dirac delta-function distribution. The density
of induced current ${\bf J}$ in Eq. (61) is given by

$\displaystyle\displaystyle{\bf J}=-{\rm
i}\omega\varepsilon_{0}(\varepsilon-1){\bf E}.$ (63)

The optical force is widely used for modeling of particle dynamics in colloid
solutions. However, when it comes to calculation of the total force acting on
the sphere:

$\displaystyle\displaystyle{\bf F}^{\rm tot}=\frac{1}{2}{\rm
Re}\int(\rho^{*}{\bf E}+\mu_{0}{\bf J}^{*}\times{\bf H}){\rm d}V,$ (64)

or the compression force of the surface charges:

$\displaystyle{\bf F}^{\rm com}={\bf F}^{\rm tot}-{\bf F}^{\rm opt},$ (65)

Mie theory naturally fails, as its solution for electric field ${\bf E}$ is
discontinuous at $r=R$ and, hence, undefined at the location of the induced
surface charge. It makes calculation of the Coulomb force in ${\bf F}^{\rm
tot}$ mathematically impossible. To resolve this issue, a continuous
transition layer with radially varying $\varepsilon(r)$ between the sphere and
the embedding medium is required. Introduction of such a layer disperses the
induced charge and makes it volume, while the electric field distribution
becomes continuous. This enables description of the forces acting on the
sphere and other characteristics of momentum transfer.

Obviously, introduction of a transition layer modifies scattering and
absorption of light. In radially inhomogeneous media, Eqs. (6) and (7) for
radial distributions change to

$\displaystyle\displaystyle\frac{{\rm d}^{2}(rH_{lm})}{{\rm
d}r^{2}}-\left[\frac{l(l+1)}{r^{2}}-k_{0}^{2}\varepsilon(r)\right](rH_{lm})=$
$\displaystyle\displaystyle\hskip
113.81102pt\frac{1}{\varepsilon(r)}\frac{{\rm d}\varepsilon(r)}{{\rm
d}r}\frac{{\rm d}(rH_{lm})}{{\rm d}r},$ (66)
$\displaystyle\displaystyle\frac{{\rm d}^{2}(rE_{lm})}{{\rm
d}r^{2}}-\left[\frac{l(l+1)}{r^{2}}-k_{0}^{2}\varepsilon(r)\right](rE_{lm})=0.$
(67)

By introducing the relative TM wave impedance $\tilde{Z}_{l}$ and TE wave
admittance $\tilde{Y}_{l}$ normalized with the free-space impedance $Z_{0}$:

$\displaystyle\displaystyle\tilde{Z}_{l}=\frac{Z_{l}}{Z_{0}}=-\frac{{\rm
i}}{k_{0}\varepsilon(r)}\frac{1}{rH_{lm}}\frac{{\rm d}(rH_{lm})}{{\rm d}r},$
(68) $\displaystyle\displaystyle\tilde{Y}_{l}=Z_{0}Y_{l}=\frac{{\rm
i}}{k_{0}}\frac{1}{rE_{lm}}\frac{{\rm d}(rE_{lm})}{{\rm d}r},$ (69)

Eqs. (66) and (67) can be rewritten for $\tilde{Z}_{l}$ and $\tilde{Y}_{l}$:

$\displaystyle\displaystyle\frac{{\rm d}\tilde{Z}_{l}}{{\rm d}r}-{\rm
i}k_{0}\left[1-\varepsilon(r)\tilde{Z}_{l}^{2}-\frac{l(l+1)}{k_{0}^{2}r^{2}\varepsilon(r)}\right]=0,$
(70) $\displaystyle\displaystyle\frac{{\rm d}\tilde{Y}_{l}}{{\rm d}r}+{\rm
i}k_{0}\left[\varepsilon(r)-\tilde{Y}_{l}^{2}-\frac{l(l+1)}{k_{0}^{2}r^{2}}\right]=0.$
(71)

Solving them with $\tilde{Z}_{l}$ and $\tilde{Y}_{l}$ given by internal fields
(18) and (19) at the inner boundary $r=R_{i}$ of the transition layer, we can
get the distributions of $\tilde{Z}_{l}$ and $\tilde{Y}_{l}$ across the entire
layer. Then, continuity of the wave impedance and admittance at the outer
boundary $r=R_{e}$ gives us the modified Mie scattering coefficients:

$\displaystyle\displaystyle
a_{l}=\frac{q_{i}\tilde{R}_{e}\psi_{l}(q_{i})\psi_{l}^{\prime}(q_{e})-q_{e}\tilde{R}_{i}\psi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})+F_{l}\psi_{l}(q_{e})}{q_{i}\tilde{R}_{e}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{e}\tilde{R}_{i}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})+F_{l}\xi_{l}(q_{e})},~{}~{}~{}$
(72) $\displaystyle\displaystyle
b_{l}=\frac{q_{e}\tilde{R}_{i}\psi_{l}(q_{i})\psi_{l}^{\prime}(q_{e})-q_{i}\tilde{R}_{e}\psi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})+G_{l}\psi_{l}(q_{e})}{q_{e}\tilde{R}_{i}\psi_{l}(q_{i})\xi_{l}^{\prime}(q_{e})-q_{i}\tilde{R}_{e}\xi_{l}(q_{e})\psi_{l}^{\prime}(q_{i})+G_{l}\xi_{l}(q_{e})},~{}~{}~{}$
(73)

with $\tilde{R}_{i}=2R_{i}/(R_{i}+R_{e})$,
$\tilde{R}_{e}=2R_{e}/(R_{i}+R_{e})$, $q_{i}=k_{i}R_{i}$, $q_{e}=k_{e}R_{e}$
and the factors

$\displaystyle\displaystyle F_{l}=-2{\rm
i}\Delta\tilde{Z}_{l}\frac{q_{i}q_{e}}{k_{0}(R_{i}+R_{e})}\psi_{l}(q_{i}),$
(74) $\displaystyle\displaystyle G_{l}=2{\rm
i}\Delta\tilde{Y}_{l}\frac{k_{0}R_{e}R_{i}}{R_{i}+R_{e}}\psi_{l}(q_{i})$ (75)

accounting for the the transition layer impedance
$\Delta\tilde{Z}_{l}=\tilde{Z}_{l}(R_{e})-\tilde{Z}_{l}(R_{i})$ and admittance
$\Delta\tilde{Y}_{l}=\tilde{Y}_{l}(R_{e})-\tilde{Y}_{l}(R_{i})$.

For thin layers, $\Delta\tilde{Z}_{l}$ and $\Delta\tilde{Y}_{l}$ are
proportional to thickness of the transition, $\Delta R=R_{e}-R_{i}$.
Therefore, effects of $\Delta\tilde{Z}_{l}$ and $\Delta\tilde{Y}_{l}$ on
scattering coefficients are often weak, except the case when
$\Delta\tilde{Z}_{l}$ become resonant. The latter happens when the transition
layer has a point $r=r_{0}$ where ${\rm Re}~{}\varepsilon(r_{0})=0$. In this
case, the last term in Eqs. (70) exhibits a complex pole that contributes

$\Delta\tilde{Z}_{l}^{\rm res}=-{\rm
i}\int^{r_{0}+0}_{r_{0}-0}\frac{l(l+1)}{k_{0}r^{2}\varepsilon(r)}{\rm d}r$
(76)

to $\Delta\tilde{Z}_{l}$ with a non-zero real part

${\rm Re}~{}\Delta\tilde{Z}_{l}^{\rm
res}\approx\pi\frac{l(l+1)}{k_{0}r_{0}^{2}}\left(\frac{{\rm d}{\rm
Re}~{}\varepsilon}{{\rm d}r}\right)^{-1}_{r=r_{0}}.$ (77)

Thus, presence of the resonant point $r_{0}$ inside the transition enhances
its intrinsic absorption for the TM spherical harmonics, especially with high
$l$. This, however, requires a dielectric permittivity $\varepsilon(r)$
changing sign of its real part inside the transition layer. The above
condition is fulfilled for ${\rm Re}~{}\varepsilon_{i}/\varepsilon_{e}<0$. In
such structures, the induced electric fields experience the localized
resonance on longitudinal polarization eigenmodes, such as phonons, plasmons
or excitons, for every orbital index $l$. The enhanced fields increase
superficial absorption in the sphere and change the volume contributions to
$\sigma_{\rm sca}$ and $\sigma_{\rm abs}$ [17].

Figure 11: Ratios $\eta_{\rm sca}$ of the modified scattering cross-section
contributions given by scattering coefficients (72) and (73) to the classical
ones for the TM and TE scattered fields with $l=1,2,3$ under
$k_{0}(R_{i}+R_{e})/2=1$, $\varepsilon_{e}=2.25$, ${\rm
Im}~{}\varepsilon_{i}=0.1$ for linear profiles of $\varepsilon(r)$.

The typical effects of transitions on the modified contributions of the TM and
TE spherical harmonics to $\sigma_{\rm sca}$ and $\sigma_{\rm abs}$ are shown
in Figs. 11 and 12 for linear profiles of $\varepsilon(r)$. Governed by the
surface-to-volume ratio, continuous transitions have stronger effect at larger
$\Delta R$ and smaller $R_{i}$.

###  IV.3 Scatterer localization

Noteworthy that neither the account for a current source nor the inclusion of
a transition layer considered above solves the issue of source-free scattered
fields. In both the cases, the spherically outgoing part of the incident
fields can still annihilate with the induced scattered fields following the
trivial solution of Maxwell’s equations. This annihilation unphysical from the
electrodynamics point of view comes mathematically possible owing to
differential equations (6) and (7) that govern distributions of both the
incident and scattered fields in the external domain. To understand how this
annihilation can be avoided, we recall the fundamental difference in
definition of the incident and scattered fields.

Conventional solution to any scattering problem consists of two steps. The
first step is consideration of the problem without the scatterer to get the
field distribution for incident radiation. In Mie theory, it corresponds to
finding the source-free incident fields supported by a uniform external
medium. The second step is solution of the problem with the scatterer included
for the given incident fields found at the first step. In other words,
incident and scattered fields are defined in two structurally different
systems of (i) the boundless external medium and (ii) the boundless external
medium with the scatterer. In Mie theory, these two systems have the common
domain at $r>R$, where both the incident and scattered fields fulfill the same
equations and, therefore, can completely compensate each other. To avoid this,
the systems (i) and (ii) should not have any common exterior domain. This
requires the scatterer to be extended over the entire space and vanish only at
$r\rightarrow\infty$. This means, the scatterer should be boundless. It still
can have a major localization in some part of space, but never be fully
confined in any finite domain.

Figure 12: Ratios $\eta_{\rm abs}$ of the modified absorption cross-section
contributions given by scattering coefficients Eqs. (72) and (73) to the
classical ones for the TM and TE scattered fields with $l=1,2,3$ under
$k_{0}(R_{i}+R_{e})/2=1$, $\varepsilon_{e}=2.25$, ${\rm
Im}~{}\varepsilon_{i}=0.1$ for linear profiles of $\varepsilon(r)$.

In other words, the origin of source-free scattered fields observed in Mie
theory is the finite localization of the scatterer at $r<R$. Continuous
extension of the scatterer beyond $r=R$ with weak contributions at $r>R$
solves the issue of electrodynamically unphysical scattered fields. However,
solution of the scattering problem for a radially inhomogeneous sphere
significantly differs from the classical Mie theory for a piecewise
homogeneous sphere, so the scattered fields can no longer be described with
Eqs. (16) and (17). Therefore, the classical scattering coefficients cannot be
modified to account for an infinitely extended scatterer.

##  V Conclusion

In conclusion, we have performed the spherical harmonic analysis of the
scattered fields given by Mie theory and revealed the two groups of source-
free and current-sourced fields it deals with. These fields have been shown to
experience distinct dynamics and influence by the sphere. We have demonstrated
how by interfering with each other, these fields make scattering and
absorption resonant and, thus, lead to all the spectral features observed in
Mie theory. We have discussed the limitations caused by the use of source-free
fields, as well as considered the ways for further refinement of the classical
theory for description of the excitation source, sphere interface and
scatterer localization. The limitations and solutions discussed have been
shown to be common to most of scattering problems and equally applicable to
non-spherical scatterers.

## References

  * [1] J. A. Stratton, Electromagnetic Theory, (Wiley-IEEE Press, 2007).
  * [2] J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, 1999).
  * [3] C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, 3rd ed. (Wiley-VCH, 1998).
  * [4] G. Mie, Annalen der Physik. 330, 377–445 (1908).
  * [5] J. W. Strutt, Phil. Mag. 41, 107–120 (1871).
  * [6] A. L. Aden and M. Kerker, J. Appl. Phys. 22, 1242–1246 (1951).
  * [7] R. Ruppin, J. Chem. Phys. 76, 1681–1684 (1982).
  * [8] F. J. Garcia de Abajo, Phys. Rev. B 59, 3095–3107 (1999).
  * [9] Plasmonic Nanoelectronics and Sensing (H.-S. Chu, E.-P. Li, eds). (Cambridge University Press, 2014).
  * [10] Nanophotonics and Plasmonics: An Integrated View (C. E. Png, Y. Akimov, eds) (CRC Press, 2017).
  * [11] T. Liu, R. Xu, P. Yu, Z. Wang and J. Takahara, Nanophotonics 9, 1115–1137 (2020).
  * [12] Light, Plasmonics and Particles (M. P. Mengüç, M. Francoeur, eds) (Elsevier, 2023).
  * [13] Y. A. Akimov, W. S. Koh and K. Ostrikov, Opt. Express 17, 10195–10205 (2009).
  * [14] P. Grahn, A. Shevchenko and M. Kaivola, New J. Phys. 14, 093033 (2012).
  * [15] A. E. Miroshnichenko et al., Nat. Commun. 6, 8069 (2015).
  * [16] B. Luk’yanchuk, R. Paniagua-Dominguez, A. I. Kuznetsov, A. E. Miroshnichenko and Y. S. Kivshar, Phil. Trans. R. Soc. A 375, 20160069 (2017).
  * [17] Y. A. Akimov, Plasmonics 7, 495–500 (2012).
  * [18] K. Kolwas and A. Derkachova, J. Quant. Spectrosc. Radiat. Transf. 114, 45–55 (2013).