# Navigating with Stability: Local Minima, Patterns, and Evolution in a Gradient Damage Fracture Model

Mourad Terzi, Oguz Umut Salman, Damien Faurie, Andrés Alessandro Leon Baldelli

► To cite this version:

Mourad Terzi, Oguz Umut Salman, Damien Faurie, Andrés Alessandro Leon Baldelli. Navigating with Stability: Local Minima, Patterns, and Evolution in a Gradient Damage Fracture Model. 2024. [hal-04696844](https://hal HAL Id: cnrs.hal.science/hal-04696844v1)

Preprint submitted on 13 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Navigating with Stability: Local Minima, Patterns, and Evolution in a Gradient Damage Fracture Model

M. M. Terzi¹, O. U. Salman¹,², D. Faurie¹, and A. A. León Baldelli³

¹LSPM, CNRS UPR3407, Université Sorbonne Paris Nord, 93400, Villateneuse, France

²Lund University, Department of Mechanical Engineering Sciences, Lund, Sweden

³CNRS, Institut Jean Le Rond d'Alembert, Sorbonne University, UMR 7190, 75005, Paris, France

September 13, 2024

## Abstract

In phase-field theories of brittle fracture, crack initiation, growth and path selection are investigated using non-convex energy functionals and a stability criterion. The lack of convexity with respect to the state poses difficulties to monolithic solvers that aim to solve for kinematic and internal variables, simultaneously. In this paper, we inquire into the effectiveness of quasi-Newton algorithms as an alternative to conventional Newton-Raphson solvers. These algorithms improve convergence by constructing a positive definite approximation of the Hessian, bargaining improved convergence with the risk of missing bifurcation points and stability thresholds. Our study focuses on one-dimensional phase-field fracture models of brittle thin films on elastic foundations. Within this framework, in the absence of irreversibility constraint, we construct an equilibrium map that represents all stable and unstable equilibrium states as a function of the external load, using well-known branch-following bifurcation techniques. Our main finding is that quasi-Newton algorithms fail to select stable evolution paths without exact second variation information. To solve this issue, we perform a spectral analysis of the full Hessian, providing optimal perturbations that enable quasi-Newton methods to follow a stable and potentially unique path for crack evolution. Finally, we discuss the stability issues and optimal perturbations in the case when the damage irreversibility is present, changing the topological structure of the set of admissible perturbations from a linear vector space to a convex cone.

## 1 Introduction

Numerous physical phenomena in materials science, such as crystal plasticity, phase transitions, twinning [18], and fracture [49, 6, 7], can be described by non-linear energy functionals at the mesoscale. The configurational variables within these energy functionals evolve under external loading, navigating equilibrium states. These states correspond to critical points of the energy functional, satisfying both boundary conditions and an optimality criterion. This optimization process is achieved through incremental minimization along a loading program. Outcomes of such optimization are fields (e.g., displacement, strain, stress), energy components, and order parameters depending on the model considered. The corresponding microstructures are crucial to understand and improve the mechanical behavior of materials.

Functionals of the type $\Psi(u)$, non-convex in their argument $u$ (a displacement field) are frequently employed in theories such as quasi-continuum methods, the multi-well Landau-type theory of weak or reconstructive phase transformations, twinning, and crystal plasticity [75, 48, 20, 26, 18, 67, 3, 5]. On the other hand, a second type of functionals, denoted as $\Psi(u, \alpha)$, find application in phase-field theories. Here, the scalar phase-field variable $\alpha$ is an internal variable that elucidates the substance's state, encompassing aspects like crystal structure, symmetry, lattice orientation, [26, 63, 38] or serving as a damage parameter in the variational phase-field theory of fracture [49, 65].

In both cases, one deals with the problem of finding configurations that satisfy $\min_u \Psi(u)$, $\min_{u,\alpha} \Psi(u, \alpha)$, or at least some necessary conditions for energy optimality and, potentially, satisfying constraints. The multiple minima of the functionals and the multitude of equilibrium states accessible during loading
<table><caption>Table 1: Notation and Definitions. The subscripted *t* in $H_t^1(0, 1)$ indicates a *t*-parametrised boundary datum</caption><thead><tr><th>Description</th><th>Symbol</th><th>Remarks</th></tr></thead><tbody><tr><td>Energy stiff</td><td>Ψ</td><td>-</td></tr><tr><td>Energy compliant</td><td>~Ψ</td><td>-</td></tr><tr><td>Load</td><td>ε̃<sub>t</sub></td><td></td></tr><tr><td>Homogeneous energy, stiff</td><td>Ψ<sup>hom</sup></td><td></td></tr><tr><td>Homogeneous energy, compliant</td><td>~Ψ<sup>hom</sup></td><td></td></tr><tr><td>Material functions</td><td>E(α), w(α)</td><td></td></tr><tr><td>Material length</td><td>l̄</td><td>dimensional, [m]</td></tr><tr><td>Film thickness</td><td>h</td><td>dimensional, [m]</td></tr><tr><td>Characteristic size</td><td>L</td><td>dimensional, [m]</td></tr><tr><td>Young modulus, substrate</td><td>Ẽ<sub>2d</sub></td><td></td></tr><tr><td>Young modulus, film</td><td>E<sub>2d</sub></td><td></td></tr><tr><td>State, stiff</td><td>u, α</td><td></td></tr><tr><td>State, compliant</td><td>u, α, v</td><td></td></tr><tr><td>Equilibrium state, stiff</td><td>y<sub>t</sub> = (u<sub>t</sub>, α<sub>t</sub>)</td><td></td></tr><tr><td>Equilibrium state, compliant</td><td>y<sub>t</sub> = (u<sub>t</sub>, α<sub>t</sub>, v<sub>t</sub>)</td><td></td></tr><tr><td>Homogeneous equilibrium state</td><td>y<sup>hom</sup></td><td></td></tr><tr><td>State &amp; perturbations, stiff</td><td>(u, α), (w, β)</td><td></td></tr><tr><td>State &amp; perturbations, compliant</td><td>(u, α, ~u), (v, β, ~v)</td><td></td></tr><tr><td>Admissible perturbations, stiff</td><td>V</td><td>H<sup>1</sup>(0, 1) × H<sup>1</sup>(0, 1)</td></tr><tr><td>Admissible perturbations, compliant</td><td>Ṽ</td><td>H<sup>1</sup>(0, 1) × H<sup>1</sup>(0, 1) × H<sup>1</sup>(0, 1)</td></tr><tr><td>State space, stiff</td><td>X<sub>t</sub></td><td>H<sup>1</sup><sub>t</sub>(0, 1) × H<sup>1</sup>(0, 1)</td></tr><tr><td>State space, compliant</td><td>X̃<sub>t</sub></td><td>H<sup>1</sup><sub>t</sub>(0, 1) × H<sup>1</sup>(0, 1)</td></tr><tr><td>Homogeneous space</td><td>V<sub>0</sub></td><td>H<sup>1</sup><sub>0</sub>(0, 1) × H<sup>1</sup>(0, 1)</td></tr><tr><td>Cone of admissible perturbations</td><td>K<sub>0</sub><sup>+</sup></td><td>H<sup>1</sup><sub>0</sub>((0, 1)) × {w ∈ H<sup>1</sup>((0, 1)) : w(x) ≥ 0 a.e. x ∈ (0, 1)}</td></tr><tr><td>Dual Cone</td><td>K<sup>*</sup></td><td>{y ∈ H<sup>1</sup>(0, 1) : &lang;x, y&rang; ≤ 0, ∀x ∈ K<sub>0</sub><sup>+</sup>}</td></tr><tr><td>Eigenvalue, eigenfunction</td><td>(λ, w<sup>*</sup>)</td><td></td></tr><tr><td>First critical load</td><td>ε̃<sub>1</sub><sup>c</sup></td><td></td></tr><tr><td>Critical loads</td><td>ε̃<sub>*</sub><sup>c</sup></td><td></td></tr><tr><td>n-th critical load</td><td>ε<sub>n</sub><sup>c</sup></td><td></td></tr><tr><td>FEM Residual vectors</td><td>R, R̃</td><td></td></tr><tr><td>FEM base functions and derivatives</td><td>N<sub>i</sub>, N'<sub>i</sub></td><td></td></tr><tr><td>Eigenvalue</td><td>λ<sub>t</sub></td><td></td></tr><tr><td>FEM state</td><td>X</td><td></td></tr><tr><td>FEM new state</td><td>X'</td><td></td></tr><tr><td>FEM current equilibrium state</td><td>X<sup>*</sup></td><td></td></tr><tr><td>FEM perturbation</td><td>P</td><td></td></tr></tbody></table>

spawn many possible evolution paths. One can expect that under a quasi-static loading protocol, the system navigates among metastable states which are continuous branches of equilibria. These branches can bifurcate and intersect as well as terminate at points where the state's stability is lost. At an instability threshold, the system restabilizes in a dissipative manner through a state transition, whether smoothly or suddenly. In this quasi-static setting, how does the system choose a new locally stable equilibrium branch with lower energy? During an isolated switching event, new equilibrium branches can be determined using a steepest descent or a continuation algorithm. The path selection may suffer from indeterminacy, however, because the energy functionals in our focus are strongly nonlinear, they may lack convexity in their arguments and thus exhibit multiple local minima, or none at all. Consequently, conducting stability and bifurcation analyses becomes crucial to distinguish among the various potential solutions or evolution paths, those that are physically relevant.

Bifurcation and stability of equilibrium configurations in dynamic systems without constraints has led to a systematic investigation of local blow-up behaviors at bifurcations points in terms of linearised (canonical) representations, allowing for easier classification and analysis of the bifurcation types [37].
For systems of ODEs the criterion of bifurcation (from a fundamental solution) amounts to the study of the existence of solutions different than the fundamental one in an arbitrary neighborhood of the control parameters. Conditions of failure of the implicit function theorem [37] describe the scenario under which a system of equations can realise more than one smooth solution. Less clear is the picture in presence of nonconvexities and nonlinear constraints associated with internal variables, where quasi-static evolutionary problems defined by optimality conditions take the form of variational inequalities defining the trajectories of a system in phase space. In these scenarios, as noted in the seminal work [33], the study of bifurcation and stability is not equivalent to the existence of solutions infinitesimally near critical points in arbitrary neighborhoods of the control parameters [8].

In our context, a bifurcation condition along the system's evolution parametrized by the control parameter(s) is associated with the uniqueness of a field of vectors tangent to the trajectory in phase space. Brittle thin film system show pattern formation and complex branch transitions. In this work, we aim at characterising the stability (or observability) of static solutions (at a given control parameter) as well as to describe the evolutionary paths stemming from an initial condition. Conditions for uniqueness of the evolution path (or its non-bifurcation) reduce to the uniqueness of solutions to a boundary value problem defined for the *rates of evolution*, or equivalently, the positive definiteness of its bilinear operator in a vector space.

Stability is a conceptually different notion when constraints play a role. The loss of such a property for a stationary points of an energy functional is of paramount importance in materials science and engineering. Illustrative in this sense are Euler buckling [10], wrinkling in thin films [36], homogeneous nucleation of dislocations in a crystal [15, 59, 3, 51, 5], buckling of lattice structures [19, 9], nucleation of cracks in soft solids or in pantographic structures [62, 65], plastic avalanches in crystals or amorphous materials [86, 81, 85].

The absence of analytical solutions in strongly non-linear settings requires resorting to numerical methods for computing and predicting equilibrium configurations that correspond to the minima of an energy functional. The minimization process involves discretizing the continuum fields onto a computational grid using methods such as finite elements, finite differences, or spectral techniques [64, 46, 35]. Afterwards, an iterative solver is employed to seek equilibrium energy states, with options including the Newton-Raphson method [82], fixed-point iteration [17, 39, 73], line-search-based descent algorithms like steepest descent or conjugate gradient [72, 22], quasi-Newton methods such as the highly-efficient Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) approach [45] which involves approximating the Hessian matrix, or more recent advancements like the fast inertial relaxation engine (FIRE) [31]. These solvers iteratively refine solutions starting from an initial guess provided as part of the solution procedure. Despite their widespread application, there remains a lack of clear understanding regarding the performance of these algorithms and their effectiveness in locating local minima.

The fracture phenomenon of thin films bonded to substrates unveils a myriad of complex crack patterns, as evidenced by numerical studies [6, 2, 34, 65, 7], resembling those observed in natural contexts such as sand or dried mud [29], and even in biological structures like animal skin [61] and bark [16, 70]. These observations hold particular relevance in the domain of stretchable and flexible electronics [25, 28] including self-healing metal thin films on a flexible substrates [77]. Within this study, centered on the numerical computation of evolutionary solutions for phase-field fracture models, we consider two one-dimensional phase-field fracture models of a brittle membrane on two types of substrates: one stiff, one compliant. The first model describes a brittle thin film deposited on a stiff (rigid) substrate, while the second model involves a compliant yet unbreakable substrate that can undergo non-uniform deformations. The finite stiffness of the substrate in the second scenario leads to nontrivial qualitative differences in terms of uniqueness of the evolution path, associated with the loss of stability of the unfractured solution [6, 41, 7, 32, 87].

Despite the one-dimensional setting we adopt here which allows for analytical predictions, these models reveal a complex landscape of equilibrium states with multiple local minima. In the absence of an irreversibility constraint, bifurcation points from homogeneous solution can easily be calculated analytically and numerically, by employing continuation techniques. An *equilibrium map* can be constructed in this setup, allowing all inhomogeneous solutions connected to the homogeneous branch to be identified along with their stability. This enables us to monitor the solutions returned by various numerical optimization techniques and assess their observability. Our findings indicate that under quasi-static conditions, line-search-based descent algorithms not relying on full Hessian can fail to detect expected branch-switching events and may return solutions that persist on unstable branches, thus lacking physical relevance. We propose a remedy to this situation which involves utilizing information from the Hessian of the functional when it becomes singular. To discuss this scenario we distinguish two settings, namely
i) that in which damage is reversible and all small perturbations are admissible, and ii) the case where damage is subject to an irreversibility constraint which forbids healing. In the former scenario negative variations of damage are allowed and indeed may occur - if convenient from an energetic viewpoint. In the second setting, instead, we consider damage as a unilateral irreversible process stemming from an irreducible one-directional pointwise growth constraint. We emphasize that we consider the reversible case as a prototypical study, rather than for its general physical relevance. Not only because it allows us to construct an equilibrium map, but also because it allows to highlight on physical grounds the mathematical differences between the reversible and irreversible cases. The study of the reversible setup may still be relevant in certain phase-field damage models where irreversibility is imposed only on crack sets that exceed a given damage threshold, referred to as relaxed crack-set irreversibility [14, 42, 23], or in models with softening elastic energy without irreversibility constraint [78, 66, 65, 4].

The rest of the paper is organized as follows. In Section 2, we present one-dimensional phase-field fracture models with both rigid and compliant elastic foundations. In Section 3 we focus on the analysis of linear and non-linear stability regarding trivial solutions. In Section 4, we construct the equilibrium map in the reversible setup, discuss the selection of equilibrium branches using various numerical optimization algorithms, and explore how irreversibility affects the stability of solutions. In the final Section 5 we summarize our results.

**Notation.** We employ standard notation for scalar Sobolev spaces defined on the unit interval, such as $H^1(0, 1)$, derivatives of one-dimensional fields, and matrix indices. We indicate the $L^2(0, 1)$-inner product of functions $u, v$ by $\langle u, v \rangle = \int_0^1 uvdx$. Subscripted $t$ means $t$-parametrised quantities, superscripted ($k$) means $k$-th iterate of an iterative algorithm. We indicate with boldface letters finite element matrices and vectors. We use the prime sign to indicates spatial derivatives. Throughout the investigation we consider a 2-layer structure. Quantities related to the substrate (fields, material parameters, energies) are denoted with a tilde. We use American English spelling throughout the text.

## 2 Material, Structure, and Evolution

Two one-dimensional fracture models for a brittle thin film bonded to substrates with different mechanical properties provide a framework to investigate the evolution and stability of crack patterns under external loading, in a simple scenario where multiplicity of solutions, equilibrium bifurcations, and stability transitions interplay.

**Material Model** We consider a one-dimensional isotropic and homogeneous brittle material modelled by a state function $W(e, \alpha, \alpha')$ which, at any point $x$, depends on the local membrane strain $e(x)$ (associated with in-plane displacements $u(x)$, namely $e(x) = u'(x)$), the local damage $\alpha(x)$, and the local gradient of the damage $\alpha'(x)$. Here, the damage variable $\alpha$ is a scalar field driving material softening, bounded between 0 and 1, where 0 indicates the undamaged material and 1, the cracked material. Thus, at points where $\alpha = 0$ the material is elastic with a stiffness $E$ (its Young modulus), at points where $\alpha = 1$ the material has a crack and zero residual stiffness, whereas for intermediate damage values the material's stiffness is $0 < E\alpha < E$. The state function $W$ is defined as

$$W(e, \alpha, \alpha') := \frac{1}{2} E \mathbf{a}(\alpha) e^2 + \mathbf{w}(\alpha) + \frac{\ell^2}{2} \alpha'^2, \quad (1)$$

where $\mathbf{w}(\alpha)$ can be interpreted as the energy dissipated during an homogeneous damaging process. It is combined with a term proportional to the square of its gradient which controls the energy cost of spatial damage variations. In the first summand, $\mathbf{a}(\alpha)$ is the function that describes the material softening. For physical consistency, $\mathbf{a}(\alpha)$ is a non-negative function that is one when $\alpha = 0$ and monotonically decreases as $\alpha \to 1$, reaching zero for $\alpha = 1$. On the other hand, $\mathbf{w}(\alpha)$ is a non-negative, zero only if $\alpha = 0$, and monotonically increases with $\alpha$, reaching $\mathbf{w}(1) = 1$. The damage-dependent stress is $\sigma(\alpha) := E\alpha(\alpha)e$. The parameter $\ell$ is a characteristic length that controls the competition between localization and homogeneous damage, effectively controlling the width of damage localizations, the peak stress of the material in one-dimensional traction experiments, and - more in general - structural size effects. Specifically, both functions $\mathbf{a}(\alpha)$ and $\mathbf{w}(\alpha)$ are chosen to be quadratic, namely

$$\mathbf{a}(\alpha) = (1 - \alpha)^2, \quad \mathbf{w}(\alpha) = \mathbf{w}_1 \alpha^2, \quad (2)$$

This modelling choice is common (yet not unique) in phase-field fracture models (cf. [14, 52, 53]). In the current context, it allows damage to evolve for an arbitrarily small value of the load.
**Structural Model** The structure under consideration is a multilayer composite constituted by a brittle thin film made of the material identified by the state function W, attached to an underlying substrate which is either rigid or elastically compliant. The thin film is a one-dimensional membrane with thickness h and length L with L ≫ h, subject to a combination of imposed inelastic tensile strains $\bar{\epsilon}_t$, imposed displacements by the substrate, and loadings at the boundary. The structure's reference configuration is the interval (0, L). The substrate is modelled as a one-dimensional elastic foundation whose displacement field is $v(x)$ and strain field is $e(v) = v'(x)$. The displacement field which is elastically compatible to an homogeneous strain in the substrate is the linear function $v(x, t) = \bar{\epsilon}(t)/2(2x - 1)$, where $\bar{\epsilon}(t)$ is the applied tensile strain and t plays the role of a loading time parameter. The film is subjected to $v(x, t)$ as a tensile imposed load and to given (compatible) displacements at its free ends $x = \{0, L\}$, so that for all $t$, $u(0) = v(0, t)$ and $u(L) = v(L, t)$. Our first model describes a brittle thin film deposited on a stiff, non-deformable substrate. This model assumes the substrate is rigid, meaning, $\tilde{u}$ is a given. The elastic interaction is modelled by a distributed linear elastic foundation of stiffness K, thus the total energy of the structure is a functional $\Psi$ constructed by considering the energy of the thin film and the energy associated with the mechanical coupling between the film and the substrate. In nondimensional form, it is given by

$$ \Psi(\alpha, u) = \int_{0}^{1} \left[ \frac{1}{2} a(\alpha)(u')^{2} + \frac{1}{2\Lambda^{2}}(u - \tilde{u})^{2} + w(\alpha) + \frac{\ell^{2}}{2}(\alpha')^{2} \right] dx, \quad (3) $$

where $\Lambda^2 = \frac{E_{eff}}{K}$, $E_{eff}$ being the effective stiffness of the two dimensional membrane. Our second model involves the same brittle thin film but a compliant elastic substrate that undergoes deformation alongside the film. Unlike for the rigid substrate, the substrate's deformation is an additional unknown which affects the overall energy landscape of the system, incorporating an extra term accounting for the strain energy of the substrate. The state of this structure is identified by the triplet $y := (u, \alpha, v)$, and the nondimensional energy of the compliant system $\tilde{\Psi}(\alpha, u, v)$ reads

$$ \tilde{\Psi}(\alpha, u, v) = \int_{0}^{1} \left[ \frac{1}{2} a(\alpha)(u')^{2} + w(\alpha) + \frac{\ell^{2}}{2}(\alpha')^{2} + \frac{1}{2\Lambda^{2}}(u - v)^{2} + \frac{\rho}{2}(v')^{2} \right] dx. \quad (4) $$

Nondimensional constants appearing in the expression of the energies can be related the three-dimensional parameters of a physical mechanical system by an appropriate de-scaling, remarking that the reduced elastic model (3) can be obtained as an asymptotic limit starting from a thin multilayer system in 3d elasticity, as shown in [43]. Thus, denoting by $E_{2d}$ the two-dimensional stiffness of the thin film membrane (per unit depth), and by a tilde the stiffness and thickness of the substrate, the nondimensional quantities appearing in the expression of the energy depend on the geometric the elastic properties of the system as follows

$$ \ell := \frac{\bar{\ell}}{L}, \quad \Lambda := \frac{E_{2d} h w_1}{\tilde{E}_{2d} L^2}, \quad \rho := \frac{E_{2d} w_1}{\tilde{E}_{2d} h}, \quad (5) $$

while spatial variables and physical displacements are respectively normalized with respect to the film's length L and the displacement scale $u_0 := \frac{w_1 L}{\tilde{E}_{2d} h}$.

For brevity, we denote by $y := (u, \alpha)$ (and $y := (u, \alpha, v)$) the mechanical state of the stiff (respectively, compliant) model system, and by $H^1(0, 1)$ the (Sobolev) space of scalar real functions defined on the unit interval which are square integrable and have square integrable first derivatives. The energy $\Psi(y)$ is well-defined for pairs belonging in the Cartesian vector product space $V := H^1(0, 1) \times H^1(0, 1)$. Similarly, the energy $\tilde{\Psi}(y)$ is well-defined for triplets $(u, \alpha, \tilde{u})$ in $\tilde{V} := V \times H^1(0, 1)$. In both cases, we extend the value of the total energy to $+\infty$ whenever $\alpha < 0$ or $\alpha > 1$ without renaming the functional.

**Evolutionary model** Assuming small enough loading rates, the evolution problem of the structure can be cast in an energetic variational formulation as an incremental rate-independent quasi-static process driven by an energy-minimality principle. This allows identifying sequences of equilibrium configurations as critical states, to ascertain their stability, and determining the system's transition trajectories between different equilibrium states. Optimality conditions defining this problem are derived from the intuitive idea that a state is observable only if it is stable, and in turn, a state is stable only if it is a local minimum of the energy among admissible state perturbations, at fixed load.

Differently from dissipative evolutions driven by an energy gradient flow whereby a system reaches equilibrium conditions through a gradient descent process parametrized by an internal timescale, the quasi-static evolution we consider is a sequence of attained equilibrium states as a necessary condition for local energy minimality. The rate-indipendency further implies that the system does not exhibit
internal timescales. As such, energy minimisation is performed at any given value of the load, and at each increment of the external load its configuration evolves subject to imposed boundary conditions and possible internal constraints.

Specifically, we consider an evolution during the loading interval $t \in [0, T]$ as a time-parametrized mapping $t \mapsto y_t$ such that, for all $t \in [0, T]$ the realised (observed) state of the system $y_t$ is a local energy minimum among all admissible state perturbations, or with respect to all admissible competitor states. In practice, given an initial condition $y_0$ at $t = 0$ we seek a state $y_t := (u, \alpha)_t$ (and $\tilde{y}_t := (u, \alpha, v)_t$ respectively, for the compliant model) such that, for a given value of the control parameter $t$, it satisfies time-dependent kinematic boundary conditions on the displacement variable and is locally energy minimal. For definiteness, denoting by $X_t = \{v \in H^1(0, 1) : v = \tilde{u}(x, t), \text{ for } x = 0 \text{ and } x = 1\} \times H^1(0, 1)$ the affine vector space of (kinematically) admissible states for the stiff model, (respectively $\tilde{X}_t = X_t \times H^1(0, 1)$, for the compliant model) we seek

$$y_t \in X_t: \quad \Psi(y_t) \le \Psi(y), \quad \forall \text{admissible competitors } y, \qquad (6)$$

and similarly for the compliant model substituting $\Psi$ with $\tilde{\Psi}$ and $X_t$ with $\tilde{X}_t$.

**Admissibility of competitors and perturbations** The admissibility of state competitors explicitly depends on the loading parameter through the kinematic boundary conditions (on displacement) and on internal constraints, namely whether damage (and hence the softening material behaviour) evolves in a reversible or irreversible manner. In the first case, as damage can evolve freely within the interval $[0, 1]$, all admissible states in $X_t$ (respectively, in $\tilde{X}_t$) are also admissible competitors. In the second case, the damage field is subject to the pointwise irreversibility constraint $\dot{\alpha}(x) \ge 0, \forall x \in [0, 1]$, requiring that the damage can only increase or stay constant. As a consequence, irreversibility restricts the admissible set of competitors to the set $K_{\alpha_t}^+ := H^1(0, 1) \times \{\beta \in H^1(0, 1) : \beta \ge \alpha_t\}$ (respectively, $\tilde{K}_{\alpha_t}^+ := K_{\alpha_t}^+ \times H^1(0, 1)$). Remark that in the definition of the competitor space $\alpha_t$ (the damage field at time $t$) is unknown at time $t$. In the irreversible case the set of admissible competitor states depends explicitly on the entire history of the evolution through the current damage field. To draw the attention to the consequences of irreversibility on the system's transitions between different equilibrium states we develop the global variational inequality (6) (and the analogous for the compliant model) by expanding the energy around the state $y_t$. Admissible perturbations in the fully reversible case they belong to the linear space $X_0$ associated with $X_t$ (respectively, $\tilde{X}_t$), whereas in the irreversible case they are bound to the closed, pointed, convex cone $K_0^+$ (respectively, $\tilde{K}_0^+$).

An energy expansion in the vicinity the state $y_t$ reads

$$\Psi(y) - \Psi(y_t) = \delta\Psi(y_t)(y - y_t) + \frac{1}{2}(y - y_t)^T \delta^2\Psi(y_t)(y - y_t) + o(\|y - y_t\|^2),$$

(and analogously for the compliant model), which allows to write first and second order necessary (and sufficient) conditions for optimality.

In the subsequent section, we exploit the energy-stability inequality (6) for both mechanical models introduced above, exploring the equilibrium configurations and identifying the conditions under which the system transitions between different equilibrium branches. This analysis will involve both analytical and numerical methods, with a focus on understanding equilibrium bifurcations, energy transitions, and stability properties of these phase-field fracture models.

## 3 Linear and nonlinear stability

Solutions to the incremental evolutionary problem are sought by solving first and second order necessary conditions for optimality encoded in the global variational inequality (6) which reduces to a linear stability problem in a linear vector space for the case of fully reversible damage, and constitutes an instance of a nonlinear stability problem in a convex cone in the presence of irreversibility.

**Linear stability in the reversible case - stiff substrate** Equilibrium equations of the system are obtained as first order necessary conditions for energy minimality, satisfying imposed displacements at the ends of the film.
The first order variation of the energy functional $\Psi$ in the direction $z := (w, \beta)$ is given by the following linear form

$$ \delta\Psi(u, \alpha)(w, \beta) = \int_{0}^{1} \left[ \mathbf{a}(\alpha)u'u' + \frac{1}{\Lambda^2}(u-u_0)w + \left( \frac{1}{2}u'^2\mathbf{a}'(\alpha) + \mathbf{w}'(\alpha) \right)\beta + \ell^2\alpha'\beta' \right] dx, \quad (7) $$

where $z \in X_0$ is a test function (an admissible perturbation) in the linear space associated with $X_t$. First order minimality conditions are local conditions that associated with the stationarity or the energy functional, that is, its first order variation should vanish for all admissible test functions, namely

$$ \delta\Psi(u, \alpha)(w, \beta) = 0, \quad \forall z \in X_0. \quad (8) $$

By using standard arguments of the calculus of variations, localizing the integral and choosing $\beta = 0$ first, and then $v = 0$ leads to establishing the strong form of local equilibrium conditions which couple the mechanical equilibrium and the damage criterion, respectively given by

$$ \left\{ \begin{array}{ll} 2(1-\alpha)\alpha'u' + (1-\alpha)^2u'' - \frac{1}{\Lambda^2}(u-\tilde{u}) & = 0, \quad x \in (0,1) \\ -\ell^2\alpha'' - (1-\alpha)(u')^2 + 2\alpha & = 0, \quad x \in (0,1). \end{array} \right. \quad (9) $$

This differential system is equipped with essential kinematic boundary conditions $u(0) = \tilde{u}(0, t)$, $u(1) = \tilde{u}(1, t)$. Conversely, conditions on the internal field $\alpha'(0) = \alpha'(1) = 0$ naturally follow by minimality. Notice that the choice of boundary conditions for displacements compatible with the substrate's deformation implies that the pair $y^{\text{hom}}(\bar{\epsilon}) := (u_h(x), \alpha_h)(\bar{\epsilon})$ given by $u_h(x) = \tilde{u}(x, t)$ and $\alpha_h$ a load-dependent constant to identify, is always (the unique homogenous) solution to the first order equilibrium equations. This makes it immediate to identify the fundamental homogeneous solution branch $t \mapsto y^{\text{hom}}$ and to decouple the elasticity problem from the evolution of damage. Therefore, the solution to (9) such that $u''(x) = \alpha'(x) = 0, \forall x \in (0,1)$ is the homogeneous branch

$$ u_h(x) = \frac{\bar{\epsilon}_t}{2}(2x-1), \quad \alpha_h = \frac{\bar{\epsilon}_t^2}{2+\bar{\epsilon}_t^2}. \quad (10) $$

Notice that, first order conditions identify the critical load threshold that activates the damaging process. The quantity $\bar{\epsilon}_*^c$, the value of the load for which the undamaged elastic state becomes unstable, is determined injecting the elastic solution $u_h$ in the energy and computing $\delta\Psi(u_h, 0)(0, \beta) = 0$. According to our mechanical energy model, the critical load is given by

$$ \bar{\epsilon}_*^c = \sqrt{2 \frac{w'(0)}{a'(0)}}, \quad (11) $$

which, by our choice of material model (whereby $w(\alpha)$ is quadratic, cf (2)) we get that $\bar{\epsilon}_*^c = 0$ so that damage necessarily starts as soon as $t > 0$. The effective total energy along the homogeneous branch reads

$$ \Psi^{\text{hom}}(\bar{\epsilon}) := \Psi(y^{\text{hom}}(\bar{\epsilon})) = \frac{\bar{\epsilon}^2}{2 + \bar{\epsilon}^2}. \quad (12) $$

An equilibrium configuration $y_t := (u, \alpha)_t$ is a state such that the first variation $\delta\Psi(y_t)(z)$ vanishes for all admissible test fields in the vector space $X_0$. To assess the incremental stability of the homogeneous solution in the reversible (linear) case, we examine the positivity of the second variation, requiring

$$ \delta^2\Psi(y_t)(y - y_t, y - y_t) > 0, \quad \forall y \in X_0, \quad (13) $$

The second directional derivative of the energy is given by the following bilinear form

$$ \delta^2\Psi(u, \alpha)(v, \beta) = \int_0^1 \left[ (1-\alpha)^2 v'^2 + \frac{1}{\Lambda^2}v^2 \right] dx + \int_0^1 \left[ -4(1-\alpha)u'v'\beta + (2+u'^2)v^2 + \ell^2\beta'^2 \right] dx, \quad (14) $$

which is well defined for perturbations $(v, \beta) \in X_0$. We extract information on the onset of instability seeking a solution in Fourier series of the fields $v$ and $\beta$ in (14), such that

$$ v(x) = \sum_{n=1}^{\infty} a_n \sin(n\pi x + \phi_n), \quad \beta(x) = \sum_{n=1}^{\infty} b_n \cos(n\pi x + \psi_n). $$
<figure><img src="image_1.png" /><figcaption>Figure 1: Computation of the determinant det $\mathcal{H}$ for the homogeneous solution $\alpha_h$, utilizing parameter values $\ell = 0.16$ and $\Lambda = 0.34$. The computation is performed for two scenarios: (a) rigid foundation and (b) compliant foundation with the additional parameter $\rho = 0.5$. The closed loops observed indicate a re-entry behavior for large deformations. Black dots mark bifurcation from the homogeneous solution while red and green dots indicate critical strains $\bar{\epsilon}^*$ and $\bar{\epsilon}^{**}$, representing the lower and upper stability limits, respectively. The critical wave number $n_c$ for the lower limit $n_c(\bar{\epsilon}^*)$ differs from that for the upper limit $n_c(\bar{\epsilon}^{**})$.</figcaption></figure>

Then we observe that, thanks to boundary conditions, $\psi_n = \phi_n = 0$ for all natural $n$. The stability condition (13) takes the form:

$$[a_n \ b_n] \mathcal{H} \begin{bmatrix} a_n \\ b_n \end{bmatrix} = [a_n \ b_n] \begin{pmatrix} a(\alpha)(n\pi)^2 + \Lambda^{-2} & \frac{\partial a}{\partial \alpha} \bar{\epsilon}(n\pi) \\ \frac{\partial a}{\partial \alpha} \bar{\epsilon}(n\pi) & \frac{\partial^2 a}{\partial \alpha^2} \bar{\epsilon}^2 + \frac{\partial h^2}{\partial \alpha^2} + \ell^2(n\pi)^2 \end{pmatrix} \begin{bmatrix} a_n \\ b_n \end{bmatrix} > 0. \quad (15)$$

By substituting the homogeneous solution $\alpha_h$ into (15), we compute det $\mathcal{H}$ as a function of $\bar{\epsilon}$ and $n$, which is depicted in Figure 1. The calculations are performed for the parameter values $\ell = 0.16$ and $\Lambda = 0.34$. The figure represents, for a given load $\bar{\epsilon}$, the wave number $n(\bar{\epsilon}) \in \mathbb{N}$ of possible energy-decreasing damage bifurcations. For an increasing loading history $\bar{\epsilon}_t \nearrow$, the wave number is non-monotonic. In Fig. 1(a), the locus det $\mathcal{H} = 0$ forms closed loops, indicative of an elastic background's influence. Notably, a re-entry behavior of the affine configuration (i.e., the homogeneous configuration re-stabilizes at large deformation) is discernible, marked by the emergence of two critical strains denoted as $\bar{\epsilon}^*$ and $\bar{\epsilon}^{**}$ (with $\bar{\epsilon}^* < \bar{\epsilon}^{**}$), representing the lower and upper stability limits for the homogeneous state. These critical points are highlighted by red and green dots in Fig. 1(a). The critical wave number $n_c$ for the lower limit $n_c(\bar{\epsilon}^*)$ differs from the critical wave number for upper limit $n_c(\bar{\epsilon}^{**})$. Finally, we remark that closed-form analytical solutions can be provided for the critical wave number and critical strains and the parametric dependence of the corresponding bifurcation thresholds can be obtained, as detailed in [65].

**Linear stability in the reversible case - compliant substrate** Denoting by $z := (v, \beta, \tilde{v}) \in \tilde{X}_0$ a test function for the state triplet $y_t := (u, \alpha, \tilde{u})_t$ at time $t$, the first order variation of the energy functional $\tilde{\Psi}$ is given by the following linear form

$$\delta \tilde{\Psi}(u, \alpha, \tilde{u})(v, \beta, \tilde{v}) = \int_0^1 [(1-\alpha)^2 u'v' + \frac{1}{\Lambda^2}(u-v)(v-\tilde{v}) + \rho \tilde{u}'\tilde{v}' + \frac{1}{2}u'^2(a'(\alpha) + w'(\alpha))\beta + \ell^2\alpha'\beta'] dx, \quad (16)$$

The Euler-Lagrange equations read

$$\begin{cases} 2(1-\alpha)\alpha'u' + (1-\alpha)^2 u'' - \frac{1}{\Lambda^2}(u-\tilde{u}) = 0, \\ -\ell^2\alpha'' - (1-\alpha)(u')^2 + 2\alpha = 0, \\ \rho\tilde{u}'' + \frac{1}{\Lambda^2}(u-\tilde{u}) = 0. \end{cases} \quad (17)$$

It is easy to show that the homogeneous solution on the trivial branch remains the same as in the case of the rigid substrate, that is $\alpha_h(\bar{\epsilon}) = \frac{\bar{\epsilon}^2}{2+\bar{\epsilon}^2}$, whereas the effective elastic energy along the homogeneous branch now reads $\tilde{\Psi}^{\text{hom}}(\bar{\epsilon}) = \frac{\bar{\epsilon}^2}{2+\bar{\epsilon}^2} + \frac{\rho}{2}\bar{\epsilon}^2$.
We once more seek the linear incremental stability of an equilibrium configuration $y_t := (u, \alpha, v)_t$ satisfying that the first order condition $\delta\tilde{\Psi}(y_t)(y - y_t) = 0$ for all admissible state perturbations $z := y - y_t \in \tilde{X}_0$ by examining the positivity of the second variation, namely

$$ (y - y_t)^T \delta^2 \tilde{\Psi}(y_t)(y - y_t) > 0, \quad \forall y - y_t \in \tilde{X}_0, $$

The second variation is given by the following bilinear form

$$ \delta^2 \tilde{\Psi}(y_t)(y, y) = \int_0^1 \left[ (1 - \alpha_t)^2 v'^2 - 4(1 - \alpha_t)u'_t v'\beta + (2 + u_t'^2)\beta^2 + \ell^2\beta'^2 + \frac{1}{\Lambda^2}(v^2 + \tilde{v}^2) + \rho\tilde{v}'^2 \right] dx. \quad (18) $$

We again proceed to extract information on the onset of instability expanding in Fourier series the fields $v, \beta$ and $\tilde{v}$ appearing in (18) such that

$$ v(x) = \sum_{n=1}^{\infty} a_n \sin(n\pi x + \phi_n), \quad \beta(x) = \sum_{n=1}^{\infty} b_n \cos(n\pi x + \psi_n), \quad \tilde{v}(x) = \sum_{n=1}^{\infty} c_n \sin(n\pi x + \theta_n). $$

Similar to the previous section we claim that, whenever the first order term vanishes, the system is stable only if $\delta^2 \tilde{\Psi} > 0$ for all sufficiently smooth admissible test fields $(v, \beta, \tilde{v})$ in the vector space $\tilde{X}_0$. The stability condition takes the form

$$ [a_n \quad b_n \quad c_n] \begin{pmatrix} a(\alpha)(n\pi)^2 + \Lambda^{-2} & \frac{\partial a}{\partial \alpha} \bar{\epsilon}(n\pi) & -\Lambda^{-2} \\ \frac{\partial a}{\partial \alpha} \bar{\epsilon}(n\pi) & \frac{\partial a^2}{\partial \alpha^2} \bar{\epsilon}^2 + \frac{\partial w^2}{\partial \alpha^2} + \ell^2(n\pi)^2 & 0 \\ -\Lambda^{-2} & 0 & \Lambda^{-2} + \rho(n\pi)^2 \end{pmatrix} \begin{bmatrix} a_n \\ b_n \\ c_n \end{bmatrix} > 0. \quad (19) $$

By substituting the homogeneous solution $\alpha_h$ into (19), we compute $\det \mathcal{H}$, as depicted in Fig. 1(b). The calculations are performed for the parameter values $\ell = 0.16$, $\Lambda = 0.34$ and $\rho = 0.5$. In Fig. 1(b), we again observe the formation of closed loops and a re-entry behavior of the affine configuration is discernible. This is marked by the emergence of two critical strains denoted as $\bar{\epsilon}^*$ and $\bar{\epsilon}^{**}$, representing the lower and upper stability limits for the homogeneous state. These critical points are highlighted by red and green dots in Fig. 1(b). The critical wave number $n_c$ for the lower limit $n_c(\bar{\epsilon}^*)$ differs from the critical wave number for the upper limit $n_c(\bar{\epsilon}^{**})$. The overall behavior of the system remains the same while the first critical wave number $n_c = 2$ is now smaller suggesting a different number of cracks will appear during the loading history.

### Nonlinear stability with irreversibility constraints
Irreversibility is introduced in the variational formulation of the evolution problem as a pointwise inequality constraint. Intuitively, irreversibility plays two distinct roles along an evolution: it acts (i) as a local constraint which prevents the damage field at a given location to decrease between two subsequent loads values (both during monotonic and non-monotonic load programs), and (ii) as a global restriction of the space of admissible variations in such a way that negative perturbations of the current damage state are no longer allowed. This changes the topological structure of the set of admissible perturbations, from a linear vector space to a convex cone. To enforce irreversibility we consider only non-decreasing damage evolutions that are sufficiently smooth with respect to the loading parameter, and seek maps $t \mapsto y_t = (u_t, \alpha_t)$ such that $\dot{\alpha}_t \ge 0$, satisfying the minimality condition (6). Because the current state can only be compared to those of with equal or higher damage, the space of admissible perturbations is a convex cone strictly contained in the vector space of admissible unconstrained perturbations, namely $K_0^+ \subset X_0$. Indeed, for $(v, \beta) \in K_0^+$, then $-(v, \beta) \in X_0$ but $-(v, \beta) \notin K_0^+$. This restriction has a profound impact on the variational characterisation of local minima and bears consequences on both for the first order (equilibrium) conditions and the second order (stability) problem which become unilateral conditions. Equilibrium states $y_t = (u_t, \alpha_t)$ of the irreversible system, are hence governed by the following (first order) necessary optimality conditions taking the form of a variational inequality

$$ \delta\Psi(y_t)(y - y_t) \ge 0, \quad \forall y - y_t \in K_0^+, \quad (20) $$

which has to hold for all admissible competitor states $y$ such that $v - u_t \in H_0^1(0, 1)$ and $\beta - \alpha_t \in H^1(0, 1)$ with $\beta \ge \alpha_t$. By testing the elasticity problem for fixed damage, and the damage problem for a given displacement field, we obtain

$$ \delta\Psi(y_t)(v - u_t, 0) = 0, \quad \delta\Psi(y_t)(0, \beta - \alpha_t) \ge 0, \quad \forall y - y_t \in K_0^+. \quad (21) $$
The last two relations are, respectively, the weak form of the mechanical equilibrium conditions and the evolution law for the damage field. The former allows to compute the elastic equilibrium displacement $u_t$ which can thus be eliminated from energy, while the latter governs the evolution of the damage field. Upon elimination of the kinematic field, the variational inequality (21).2 takes a particularly expressive form when written as a complementarity problem. This highlights the mechanical nature of the damage criterion as a threshold law. Namely, the strict convexity of the elastic model for given damage implies that (21).1 has, for given $\alpha$, a unique time-parametrized solution $u_t(\alpha)$. Substituting in (21).2 and accounting for the irreversibility constraint we are led to seek a map $t \mapsto \alpha_t$ such that

$$ \dot{\alpha}_t \ge 0 \quad -\phi_t(\alpha_t) \le 0 \quad \phi_t(\alpha_t) \dot{\alpha}_t = 0. \qquad (22) $$

Here, $\phi_t$ is the scalar function associated with the variation of elastic energy density at the equilibrium, defined by $\delta\Psi(y_t)(0, \beta) = \langle -\phi_t(\alpha_t), \beta \rangle$. Here, the subscript $t$ indicates that this quantity is computed for the equilibrium $u_t$. Consequently, $\phi_t(0)$ is the variation of the energy density at equilibrium for the undamaged structure, and all equilibrium solutions $u_t$ such that $-\phi_t(0) > 0$ belong to the interior of the damage yield surface for the sound structure. The equality $-\phi_t(0) = 0$, conversely, indicates that the damage criterion has been attained by the sound structure, or equivalently, that the state $(u_t, 0)$ has reached, from the interior, the boundary of the (damage-dependent) elastic domain. Explicitly, the function $\phi_t$ depends on the brittle material model and is defined, for both the stiff and compliant substrate models, as

$$ \phi(\alpha) := \frac{1}{2} \mathbf{a}'(\alpha) \epsilon^2 + \mathbf{w}'(\alpha), \qquad (23) $$

where $\epsilon$ is the elastic strain and $\mathbf{a}'(\alpha)$ and $\mathbf{w}'(\alpha)$ are the derivatives of the softening and dissipation energy densities with respect to the damage variable. The inequality (21).2 identifies the domains of admissible strains (and, by duality, of stresses) for homogeneous solutions as $\mathcal{R}(\alpha) := \{\epsilon \in R_{sym}^{2 \times 2} : \mathbf{E}_{2d} \epsilon^2 \le -\frac{2\mathbf{w}'(\alpha)}{\mathbf{a}'(\alpha)}\}$, and $\mathcal{R}^*(\alpha) := \{\sigma \in R_{sym}^{2 \times 2} : \frac{\sigma^2}{\mathbf{E}_{2d}} \le \frac{2\mathbf{w}'(\alpha)}{\mathbf{s}'(\alpha)}\}$, respectively. As a first order optimality condition, (21).2 states that the local elastic energy release is either smaller than or equal to the (marginal) cost of damage, whereas the complementarity condition $\phi_t(\alpha_t)\dot{\alpha}_t = 0$ ensures that the damage field evolves only if the energy release rate is critical. The three conditions above, established as necessary first order condition for constrained optimality, encode the pointwise non-negativity of the damage rate, the boundedness of the elastic domain, its dependence upon damage, and the complementarity between the attainment of the damage criterion and the conditions for the evolution of the internal order parameter.

In the current one-dimensional setup with homogeneous initial conditions $y_0 = (0, 0)$ and compatible kinematic boundary conditions, the existence of a homogeneous solution implies that the damage criterion is attained everywhere throughout the bar at the same load. This greatly simplifies the analysis of the energetic properties of the system. Using the elastic solution $u_t = 2\bar{\epsilon}_t(x - 1/2)$ the inequality in (21).2 yields the following algebraic inequality

$$ 0 \ge -\phi_t(\alpha_t) = -\frac{1}{2} \bar{\epsilon}_t^2 \mathbf{a}'(\alpha_t) - \mathbf{w}'(\alpha_t). \qquad (24) $$

The equality in (24) is an algebraic equation for $\alpha_t$ which identifies the evolution of the homogeneous damage response, as a function of the given load level $t$. The investigation of the stability properties requires considering second order energy variations with respect to all admissible perturbations that render null the first order term in the energy expansion (2). In the general case, this requires distinguishing between the regions where the damage criterion is attained, and thus damage can evolve, from the (complementary) domain where damage cannot evolve (there, the second relation in (22) is satisfied with a strict inequality). In our setup, the existence of nontrivial homogeneous solutions simplifies the analysis because the damage criterion is attained everywhere, damage can increase throughout the whole domain, and the function space of admissible perturbations is defined on the fixed domain $(0, 1)$.

Assume now that a state $y_t$ is known as a function of $t$ such that it solves (20) and is sufficienty smooth so that the (right) derivative with respect to $t$ is well-defined. As $t$ varies, $y_t$ describes a (smooth) curve in the phase space identified by its right tangent vector $\dot{y}_t =: \lim_{\tau \to 0^+} \frac{y_{t+\tau} - y_t}{\tau}$, the rate of evolution. A fundamental question is to discern whether $y_t$ is an isolated equilibrium state tracing a unique evolution path, or conversely if it lays at the intersection of multiple equilibrium curves. To this end, differentiating (20) with respect to $t$ we obtain a boundary value problem relating the rate of evolution $\dot{y}_t$ to the current state $y_t$, supposing the latter known, namely

$$ \delta^2 \Psi(y_t)(\dot{y}_t, \zeta - \dot{y}_t) + \delta \dot{\Psi}(y_t)(\zeta - \dot{y}_t) \ge 0, \quad \forall \zeta \in X_0, \qquad (25) $$
where $\delta\tilde{\Psi}$ is the time-derivative of the linear form corresponding to the first order energy variation. By construction, the homogeneous rate $\dot{y}_h$ is a solution of (25), the question is whether another solution exists. The uniqueness is thus ensured by the positive definiteness of the quadratic form in $X_0$. Thus, the non-bifurcation condition for the homogeneous evolution reads

$$ \delta^2\Psi(y^{\text{hom}})(\zeta, \zeta) > 0, \quad \forall \zeta \in X_0, \qquad (26) $$

which formally coincides with the classical linear stability problem (13) in the reversible case, yet has a different mechanical interpretation in terms of evolution rates and uniqueness of equilibrium branches. Remark that, in the general case in which the damage criterion is not attained everywhere, the space of admissible perturbations for the (second order) bifurcation problem is $X'_0 := H_0^1(0,1) \times \{\beta \in H^1(0,1) : \delta\Psi(y_t)(0,\beta) = 0\}$ for the stiff substrate model, and $\tilde{X}'_0 := X'_0 \times H_0^1(0,1)$ for the compliant substrate model. The first load at which the bifurcation inequality (26) fails, namely $t_b := \inf_t\{\delta^2\Psi(y_h)(\zeta, \zeta) = 0, \forall \zeta \in X_0\}$ corresponds to the first bifurcation load, namely, the load for which there exist (multiple) equilibrium curves intersecting the homogeneous branch. As a consequence, for $t \ge t_b$ the possibility exists of bifurcating away from homogeneous branch. The study of the bifurcation problem is functional to infer a partial response on the stability of the state. Indeed, if the current equilibrium branch is unique then, necessarily, the current state is stable. The converse is not true, however, as the existence multiple possible of bifurcation paths is not a sufficient condition exclude the stability of the current state. This holds true for the irreversible case, due to the conceptual difference between the bifurcation and the stability problems.

The stability of the homogeneous solution in the irreversible case, according to our energetic viewpoint, is governed by the positivity of $\delta^2\Psi(y^{\text{hom}})$ on the constrained space of admissible state perturbations $K_0^+$. Denoting $t_s$ the load at which the homogeneous solution loses stability by analogy to the bifurcation load, namely $t_s := \inf_t\{\delta^2\Psi(y^{\text{hom}})(\zeta, \zeta) = 0, \forall \zeta \in K_0^+\}$, the set inclusion $K_0^+ \subset X_0$ implies that necessarily $t_b \le t_s$. (equality occurs when the first bifurcation mode is nonnegative)

This indicates a qualitative conceptual distinction between the bifurcation and the stability thresholds, in the irreversible case. As a consequence, a system can persist along a critical non-unique equilibrium branch, yet be stable. A sufficient condition for the stability of the homogeneous state $y^{\text{hom}}$ is given by the strict positivity of the Hessian form, on the constrained space of admissible perturbations, namely (for the stiff substrate model)

$$ \delta^2\Psi(y^{\text{hom}})(y - y^{\text{hom}}, y - y^{\text{hom}}) > 0, \quad \forall y - y^{\text{hom}} \in K_0^+, \qquad (27) $$

and similarly for the compliant substrate model, by replacing $\Psi$ with $\tilde{\Psi}$ and $K_0^+$ with $\tilde{K}_0^+$. The variational inequality above is a constrained eigenvalue problem and a tool to characterize the stability of the state $y^{\text{hom}}$. Its solution yields, at load $\epsilon_t$, either a positive eigenpair $(\lambda_t, z_t^*) \in \mathbb{R}^+ \times K_0^+$ as a sufficient condition for the stability of current state, or a pair $(\lambda_t, z_t^*) \in \mathbb{R}^- \times K_0^+$ where $\lambda_t$ is the local (negative) energy curvature and the eigenmode $z_t^*$, indicating the direction of maximum energy decrease is interpreted as the instability mode, pointing the system towards an optimal direction of energy descent. From the numerical standpoint, the bifurcation eigen-problem in the vector space (25) may be regarded as an approximated version of the stability problem, in the time-discrete setting. As suggested in [7] through the notion of 'incremental-stability', the irreversibility constraint in the stability problem can be relaxed to a pointwise inequality with respect to the state at the previous time-step, denoted $y_-$.

This allows to replace the current state $y_t$ in the definition of the space of perturbations, with $y_-$. As a consequence, assuming that $\alpha_t$ is the equilibrium damage field solving first order optimality conditions at $t$ and that $\alpha_-$ is the solution at the previous load step, admissible perturbations for the second order problem (27) are all the $y - y_- \in K_0^+$. In this way, the set of perturbations is enlarged. It includes all sufficiently smooth functions $\beta$ which cancel the first order term, without restriction on the sign provided that $\alpha_t(x) + \beta(x) - \alpha_-(x) \ge 0$ for all $x \in (0,1)$. In this way, the space of perturbation allows for (small) perturbations $\beta$ which can be negative at points $x \in (0,1)$ where $\alpha_t(x) > \alpha_-(x)$. Such a space is a vector space and the eigen-problem can be solved by standard methods of linear algebra by projecting the Hessian to the set of active constraints, cf [1]. Conversely, the stability problem (27), a constrained eigenvalue problem in a convex cone, is of a different nature altogether due to the different topology of the underlying energy space. The associated discrete problem can be numerically solved by exploiting the orthogonality between the set $K_0^+$ and its dual $K^* := \{y \in H^1(0,1) : \langle x, y \rangle \le 0, \forall x \in K_0^+\}$, cf. [55, 58].
<figure><img src="image_2.png" /><figcaption>Figure 2: Equilibrium branches for phase-field thin film model of an elastic membrane on a compliant elastic foundation. Left, the energy difference $\Delta\Psi$ between the current and the homogeneous configurations, for the same load. Stability of solutions is color-coded: blue denotes stability while orange indicates instability, by numerical evaluation of the smallest eigenvalue of the stiffness matrix **H**. Equilibrium branches are parametrised by an integer *n* and rrows indicate branch switching events associated with the loss of stability; Right, damage and strain profiles of the minimum energy configurations on each branch.</figcaption></figure>

# 4 Numerical solutions

## 4.1 Identification of equilibrium branches

We now turn our attention to the solutions beyond the trivial homogeneous branch and will use numerical calculations to identify all the equilibrium branches corresponding to the inhomogeneous solutions. Our goal is to construct an *equilibrium map* that represents all stable and unstable equilibrium states as a function of the external load [57].

To find the equilibrium branches we make use of a pseudo-arclength continuation technique implemented in the software AUTO [24], see also [57]. It solves the nonlinear equations 9 and 17 in the case of rigid and compliant foundation, respectively, with the end displacement $\bar{\epsilon}_t$ treated as a continuation parameter. To discretize the boundary-value problem, it uses collocation with Lagrange polynomials, and in our simulations we had $N = 300$ mesh points with $N_c = 5$ collocation nodes and activated mesh adaptation. In order to assess the stability of equilibrium branches, we use the numerical evaluation of the smallest eigenvalue of the second variation by discretizing the integrals (14) and (18) to construct the stiffness matrices **H** and $\tilde{\mathbf{H}}$, investigating numerically the sign of the minimal eigenvalue $\lambda_t$ of the corresponding discrete quadratic form [69]. The finite element discretization of the displacement and damage field $(u, \alpha)$ with $n_u$ displacement degrees-of-freedom $\mathbf{u} = \{u_1, \dots, u_{n_u}\}^T$ and $n_\alpha$ damage degrees-of-freedom $\boldsymbol{\alpha} = \{\alpha_1, \dots, \alpha_{n_\alpha}\}^T$ is given by $u(\mathbf{x}) \approx u_{\text{FE}}(\mathbf{x}) := \sum_{i=1}^{n_u} \mathcal{N}_i^{(u)}(\mathbf{x})u_i$ and $\alpha(\mathbf{x}) \approx \alpha_{\text{FE}}(\mathbf{x}) := \sum_{i=1}^{n_\alpha} \mathcal{N}_i^{(\alpha)}(\mathbf{x})\alpha_i$ where $\mathcal{N}_i^{(u)}(\mathbf{x})$ and $\mathcal{N}_i^{(\alpha)}(\mathbf{x})$ are the finite element basis functions and $u_h$ and $\alpha_h$ are nodal values for the displacement and damage fields, respectively. We used quadratic one-dimensional finite elements with 3 nodes whose shape functions $N_i(x)$ at a node $i$ are given by $N_1(\xi) = -0.5\xi(1-\xi)$, $N_2(\xi) = -0.5\xi(1+\xi)$ and $N_3(\xi) = -(1-\xi)(1+\xi)$, where $\xi$ is the isoparametric coordinate. The discrete solution $u(x_i)$ provided at discrete nodes $x_i$ by AUTO was first interpolated
using B-spline basis functions of degree 3 [30], and then utilized to calculate the integrals (14) and 18 employing a three-point Gauss integration scheme. The fixed boundary conditions were enforced by removing the rows and columns corresponding to $x = 0$ and $x = 1$ from the stiffness matrices $\mathbf{H}$ and $\tilde{\mathbf{H}}$. The explicit form of the stiffness matrix for the stiff substrate model is given by

$$ \mathbf{H} = \begin{bmatrix} \int_{0}^{1} \left[ \frac{1}{\Lambda^2} \mathcal{N}_i \mathcal{N}_j + (1-\alpha)^2 \mathcal{N}_i' \mathcal{N}_j' \right] dx & -2 \int_{0}^{1} (1-\alpha) u' \mathcal{N}_i' \mathcal{N}_j' dx \\ -2 \int_{0}^{1} (1-\alpha) u' \mathcal{N}_i \mathcal{N}_j' dx & \int_{0}^{1} \left[ (2+u'^2) \mathcal{N}_i \mathcal{N}_j + \ell^2 \mathcal{N}_i' \mathcal{N}_j' \right] dx \end{bmatrix}, \quad (28) $$

whereas the compliant substrate model it reads

$$ \tilde{\mathbf{H}} = \begin{bmatrix} \int_{0}^{1} \left[ \frac{1}{\Lambda^2} \mathcal{N}_i \mathcal{N}_j + (1-\alpha)^2 \mathcal{N}_i' \mathcal{N}_j' \right] dx & -2 \int_{0}^{1} (1-\alpha) u' \mathcal{N}_i' \mathcal{N}_j' dx & - \int_{0}^{1} \left[ \frac{1}{\Lambda^2} \mathcal{N}_i \mathcal{N}_j \right] dx \\ -2 \int_{0}^{1} (1-\alpha) u' \mathcal{N}_i \mathcal{N}_j' dx & \int_{0}^{1} \left[ (2+u'^2) \mathcal{N}_i \mathcal{N}_j + \ell^2 \mathcal{N}_i' \mathcal{N}_j' \right] dx & 0 \\ - \int_{0}^{1} \left[ \frac{1}{\Lambda^2} \mathcal{N}_i \mathcal{N}_j \right] dx & 0 & \int_{0}^{1} \frac{1}{\Lambda^2} \mathcal{N}_i \mathcal{N}_j + r_s \mathcal{N}_i' \mathcal{N}_j' dx \end{bmatrix}. \quad (29) $$

Our objective is to establish a branch switching strategy that, as the external loading parameter monotonically increases, allows equilibrium branch transitions when the current branch ceases to be stable or to exist. This strategy must ensure the system's re-stabilization following an instability in a dissipative manner. In a quasi-static scenario, it should select a new locally stable equilibrium branch with inherently lower energy. Considering applications in structural mechanics, our approach to selecting the new equilibrium branch relies on a criterion of local energy minimization (LEM) which emulates the zero viscosity limit of overdamped viscous dynamics. This approach differs from a global energy minimizing (GEM) strategy, which may be more relevant in biomechanical applications [65]. According to LEM protocol, during quasi-static loading, the system will remain in a metastable state (a local minimum of energy) until it becomes unstable. Subsequently, during an isolated switching event, the new equilibrium branch will be chosen using a descent algorithm [60].

We display the equilibrium branches that solve the nonlinear equations 9 in Fig. 2. These branches are represented in Figure 3 by plotting the energy difference $\Delta\Psi$ between the energy of the current solution and the energy of the homogeneous solution $\Psi^{\text{hom}}(\bar{\epsilon}_t)$ at the current value of the loading parameter $\bar{\epsilon}_t$. The stability of solutions is color-coded, light blue and orange representing stable and unstable solutions, respectively. Stability is discerned through numerical evaluation of the smallest eigenvalue $\lambda_t$ of the stiffness matrix $\mathbf{H}$ (and $\tilde{\mathbf{H}}$) defined in (28) (respectively, in (29)). Eigenvalues represent the local curvatures of the energy functional, so the sign of the smallest $\lambda_t$ determines the stability of the solution, with $\lambda_t > 0$ indicating stability and $\lambda_t < 0$ indicating instability.

Under the LEM protocol, the system explores the fundamental branch where the homogeneous solution is stable until $\bar{\epsilon}_1^*$, identified using linear stability analysis, see Figure 2. At the critical load on the homogeneous branch at point $H$, a first instability determines a branch switching transition, from the trivial branch with wavenumber $n = 0$ to the nontrivial equilibrium branch with $n = 3$. The ensuing non-homogeneous configuration is linearly stable, as seen in Fig. 2. We depict the lowest energy configuration on this equilibrium branch in Fig. 2 at point $A$, consisting of two simultaneously nucleated localized cracks, one inside the domain and one on the boundary.

When the loading parameter $\bar{\epsilon}_t$ is further increased, the equilibrium configuration with $n = 3$ loses linear stability at point $A'$. In Fig. 2 we represent state transitions according to the LEM protocol by arrows. It can be observed that under this protocol, the only available transition from point $A'$ is to the branch with $n = 4$, which is locally stable within at the corresponding applied strain $\bar{\epsilon}_2$. The minimum energy configuration on the branch with $n = 4$ is depicted in Fig. 2-right at point $B$, consisting of two cracks: one inside the domain and two on the right and left boundaries.

Further increasing the load, the stability of the branch with $n = 4$ is lost at point $B'$, a single available transition leads the system to the branch with $n = 5$. While the branch with $n = 6$ appears accessible due to its lower energy compared to the current state, it is unstable at the current value of $\bar{\epsilon}_3$. The lowest energy configurations on the branches with $n = 5$ and $n = 6$ are illustrated in Fig. 2(right). Increasing the load along the branch with $n = 5$ which loses linear stability at point $C'$, a branch switching event occurs towards the branch with $n = 7$, the corresponding minimum energy configuration on this branch is depicted in Fig. 2(right). Finally, this branch reconnects to the homogeneous branch with $n = 0$ at point $F$. Remark that, in the case of a stiff substrate, the LEM protocol identifies a unique evolution path, with the system having a single bifurcation option, or available branch, at each instability point.

For the compliant substrate model, the equilibrium branches that solve the nonlinear equations (Eq. 17) are shown in Fig. 3. We again plot the energy difference $\Delta\tilde{\Psi}$ between the energy of the current solution and the energy of the homogeneous solution $\Psi^{\text{hom}}(\bar{\epsilon})$ at the current value of the loading
<figure><img src="image_3.png" /><figcaption>Figure 3: Equilibrium branches for phase-field thin film model of an elastic membrane on a compliant elastic foundation. Left, the energy difference $\Delta\tilde{\Psi}$ between the current and the homogeneous configurations, for the same load. Stability of solutions is color-coded: blue denotes stability while orange indicates instability, by numerical evaluation of the smallest eigenvalue of the stiffness matrix $\tilde{H}$. Equilibrium branches are parametrised by an integer $n$; Right, damage and strain profiles of the minimum energy configurations on each branch.</figcaption></figure>
parameter $\bar{\epsilon}$. The stability of solutions is color coded (light blue indicates stable states, orange identifies unstable states), the stability marker being given by the numerical evaluation of the smallest eigenvalue of the stiffness matrix $\bar{\mathbf{H}}$, defined in (29).

According to the LEM protocol with a loading history starting in the unloaded and sound configuration, the initial transition from the trivial solution to the only available branch with $n = 2$ occurs at $\bar{\epsilon}_1$, marked as point $H$. As illustrated in Fig. 3(right), the damage and strain profiles of the lowest energy configuration at point $A$ is characterised by two boundary cracks. As the loading increases, the system persists on the $n = 2$ branch until point $A'$ which marks an instability. The system can now access two equilibrium branches: $n = 3$ and $n = 4$ shown in gray and black arrows. Fig. 3(right) shows the typical damage and strain profiles on these branches. Branches $n = 3$ and $n = 4$ feature one bulk crack plus one boundary crack, and one bulk crack plus two boundary cracks, respectively.

The choice of the subsequent branch transition at point $A'$ will dictate the ensuing crack growth. On branches with $n = 3$ and $n = 4$, at the instability points $B'$ and $C'$, the system will once again encounter two available equilibrium branches. From $n = 3$ branch at point $B'$, the system can transition to either the $n = 4$ or $n = 5$ branch. Opting for the $n = 4$ branch, the subsequent branch selection occurs at point $C'$, offering branches with $n = 6$ or $n = 5$. Notably, the $n = 5$ branch smoothly reconnects with the trivial solution at point $F$. However, the $n = 6$ branch experiences another instability at point $D'$ before rejoining the trivial solution.

For the current choice of material parameters, the overall behavior of the compliant substrate model reveals additional complexity, with more branch switching events and the non-uniqueness of the global response. In the stiff substrate model, the LEM strategy identifies a unique evolution path, while in the compliant substrate model the system can follow eight different trajectories in phase-space, given by the various available path-bifurcation choices. This additional richness provides a good case to test numerical optimization methods, which we will discuss in the following.

## 4.2 Equilibrium branch selection in overdamped viscous dynamic

We consider now the equilibrium solutions reachable through the line-search based quasi-Newton algorithms such as conjugate gradient or the BFGS optimization, which effectively mimic the zero viscosity limit of overdamped viscous dynamics [68]. Quasi-Newton methods serve as alternatives to Newton's method for locating roots or local extrema of functions. Particularly advantageous in scenarios where computing the Hessian at each iteration is impractical or computationally expensive, these methods circumvent the need for explicit computation of energy derivatives. Instead, they rely on evaluating the function value and its gradient and updating the Hessian by analyzing successive gradient vectors.

Quasi-Newton methods are highly suitable for solving the phase-field equation of fracture, particularly when compared to standard Newton method-based monolithic solvers. Such solvers, which simultaneously solve the equations for both damage and displacement variables, often falter when confronted with nonconvex energy functionals. For example, as demonstrated in [82], the Newton method-based monolithic algorithm does not consistently handle brittle fracture scenarios involving abrupt crack propagation. Recently, quasi-Newton methods, particularly the BFGS variant, have been employed to effectively solve the system of coupled governing equations in a monolithic fashion within the phase-field method of fracture. These methods have demonstrated success in various engineering applications, as evidenced by [40, 83, 65, 47].

Our primary objective is not to provide a comprehensive assessment of quasi-Newton methods on a global scale. Instead, our focus lies in scrutinizing their behavior and performance specifically concerning equilibrium branch selections within our simplified framework, where all branches are readily identified. By narrowing our scope to this specific aspect, we aim to gain insights into the effectiveness and reliability of quasi-Newton methods in reaching stable states of our system.

Recall that quasi-Newton algorithms only evaluate the function value and its gradient to reach equilibrium configurations. This implies that in our case, we need to evaluate integrals (3) and (4), along with their first variations given by (7) and (16), in the stiff and compliant models, respectively. We discretized the integrals (7) and (16) to construct the residuals vectors using finite elements to obtain

$$ \mathbf{R} = \int_{0}^{1} [(1 - \alpha)^2 u' \mathcal{N}'_{i} + \frac{1}{\Lambda^2}(u - \tilde{u})\mathcal{N}_{i} + \frac{1}{2}(u'^{2}\mathbf{a}'(\alpha) + \mathbf{w}'(\alpha))\mathcal{N}_{i} + \ell^2 \alpha' \mathcal{N}'_{i}] dx, \quad (30) $$

in the stiff model and

$$ \tilde{\mathbf{R}} = \int_{0}^{1} [(1 - \alpha)^2 u' \mathcal{N}'_{i} + \frac{1}{\Lambda^2}(u - \tilde{u})\mathcal{N}_{i} + \frac{1}{2}(u'^{2}\mathbf{a}'(\alpha) + \mathbf{w}'(\alpha))\mathcal{N}_{i} + \ell^2 \alpha' \mathcal{N}'_{i} + \rho \tilde{u}' \mathcal{N}_{i}] dx. \quad (31) $$
<figure><img src="image_4.png" /><figcaption>Figure 4: Stiff substrate. Quasi-static loading simulations with L-BFGS: (a) the energy difference $\Delta\Psi$, between the quasi-Newton solutions and the homogeneous solutions are superimposed onto equilibrium branches (light blue: stable; orange: unstable). Bottom, the smallest eigenvalue $\lambda_t$ of the second variation $\Psi''$ as a function of the loading parameter $\bar{\epsilon}$. In red, indicated the unstable region $\lambda_t < 0$, where the equilibrium solution returned by the solver does not satisfy the evolution law; Right, damage profiles referring to the transition endstates.</figcaption></figure>

in the compliant substrate model.

Among iterative methods for large-scale unconstrained optimization, particularly when dealing with possibly dense Hessian matrices, quasi-Newton methods often emerge as preferable alternatives to the widely-used Newton-Raphson (NR) algorithm. The NR algorithm, conventionally utilized for solving linear equations to determine the correction $\Delta\mathbf{X}^{(k)}$ from the current estimate $\mathbf{X}^{(k)} = (\mathbf{u}^{(k)}, \alpha^{(k)})$ at iteration $k$, is expressed in our context as:

$$K_{ij}\Delta X_j^{(k)} + R_i = 0, \qquad (32)$$

where the discrete stiffness matrix $\mathbf{H}$ and bulk forces $\mathbf{R}$ are computed with the initial guess $\mathbf{X}^{(k)}$. Subsequently, the guess is updated as $\mathbf{X}^{(k+1)} = \mathbf{X}^{(k)} + \Delta\mathbf{X}^{(k)}$ after solving Equation (32) using LU factorization [69]. It's evident that the NR algorithm fails if the discrete stiffness matrix $\mathbf{H}$ isn't invertible.

On the other hand, quasi-Newton methods are well-established (see standard textbooks, e.g., [1, 56]), and generate a sequence $\{\mathbf{X}^{(k)}\}$ according to the following scheme:

$$\mathbf{X}^{(k+1)} = \mathbf{X}^{(k)} + h^{(k)}\mathbf{p}^{(k)}, \quad k = 0, 1, \dots \qquad (33)$$

with

$$\mathbf{p}^{(k)} = -(\mathbf{B}^{(k)})^{-1}\mathbf{R}, \qquad (34)$$

where $(\mathbf{B}^{(k)})^{-1}$ approximates the inverse of the Hessian matrix $\mathbf{H}$ and $h^{(k)}$ represents a step length. Particularly, instead of computing $(\mathbf{B}^{(k)})^{-1}$ at each iteration $k$, these methods update $(\mathbf{B}^{(k)})^{-1}$ in a straightforward manner to obtain the new approximation $(\mathbf{B}^{(k+1)})^{-1}$ for the subsequent iteration. Additionally, rather than storing full dense $N \times N$ approximations, they only retain a few vectors of length $N = n_{\alpha} + n_{u}$, enabling implicit representation of the approximations. Moreover, the choice of the step length $h^{(k)}$ is carried out through a line search to minimize a function $f(h) = f(\mathbf{X}^{(k)} + h\mathbf{p}^{(k)})$ in order to find an acceptable step size $h^{(k)}$ such that $h^{(k)} \in \arg\min_h f$.

Among quasi-Newton schemes, the L-BFGS method is widely regarded as one of the most efficient and well-suited for large-scale problems due to its limited and user-controlled storage requirements. This method relies on constructing an approximation of the inverse Hessian matrix, leveraging curvature
<figure><img src="image_5.png" /><figcaption>Figure 5: Compliant substrate. Quasi-static loading simulations with L-BFGS: (a) the energy difference $\Delta\tilde{\Psi}$, between the quasi-Newton solutions and the homogeneous solutions are superimposed onto equilibrium branches (light blue: stable; orange: unstable). Bottom, the smallest eigenvalue $\lambda_t$ of the second variation $\tilde{\Psi}''$ as a function of the loading parameter $\bar{\epsilon}$. In red, indicated the unstable region $\lambda_t < 0$ for the equilibrium solution; Right, damage profiles referring to the transition endstates.</figcaption></figure>

information solely from recent iterations. Note also that the update formula for the approximative in successive minimization steps depends on the adapted algorithm. These algorithms are extensively used in the literature, and we refer the reader to the references for a more detailed description [50, 84, 1, 56, 71, 44, 21]. However, quasi-Newton methods do present certain drawbacks, notably slow convergence for ill-conditioned problems, particularly when the eigenvalues of the Hessian matrix are widely dispersed [71].

In our numerical experiments, we employ the BFGS and L-BFGS solvers from the Python SciPy library [79] and the Alglib library [11]. These solvers utilize the residual vectors (see 30 and 31) at each finite element node, alongside the values of integrals (3) and (4). The following results were obtained using the Alglib library [11]; however, our findings remain consistent across the other libraries mentioned.

In Figure 4, we present the results of our numerical experiments, overlaying solutions obtained with the quasi-Newton algorithm onto the equilibrium map determined using the pseudo-arclength continuation method described in the previous section. Fig. 4-bottom shows the smallest eigenvalue $\lambda_t$ whose positivity informs on the stability of the solution. All the branch switching events are delayed in the sense they do not take place when the smallest eigenvalue of the second variation vanishes but only beyond the loads corresponding to the eigenvalue's sign transition. Notably, the expected first branch switching event—from the homogeneous solution to the branch with $n = 3$ as outlined in the LEM protocol—did not occur at point *H* as anticipated. During a monotonic loading, the system remains on the trivial branch, as shown in Fig. 2(a), until it transitions to the branch with $n = 4$ at a value of the load much higher than the critical loading parameter $\bar{\epsilon}_1$. By monitoring the smallest eigenvalue of the second energy variation at the solution, we observe that indeed the solutions computed via the quasi-Newton solver are unstable beyond $\bar{\epsilon}_1$, which is consistent with the linear stability results for the trivial branch. A closer inspection of the solution field reveals that the quasi-Newton solutions are not homogeneous but instead exhibit a small perturbation akin to the instability mode calculated analytically (i.e., the eigenvector associated with the smallest eigenvalue). For instance, the quasi-Newton solution before the first branch switching event is shown in inset (state *H*) in Fig. 4-right, where we see that the equilibrium damage shows a slight oscillation reminiscent of the eigenmode $\alpha_n(x) \sim \cos(n\pi x)$, with $n = 3$. The state transition captured by the quasi-Newton method is the one towards the $n = 4$ branch, with one bulk crack and two boundary cracks, see inset (A) in Fig. 4-right. If the transition had taken place at $\bar{\epsilon}_1$ as
<figure><img src="image_6.png" /><figcaption>Figure 6: Profiles of damage fields and perturbations showing the effect of the kick algorithm in state transitions (from unstable to stable). Branch switching events correspond to load values indicated by $\bar{\epsilon}_i$ with $i \in \mathbb{N}$ in Fig. 4 and 5. The stiff substrate model is in the left column, the compliant substrate model on the right and the light (orange) line indicates the unstable damage field $\alpha^*$, the eigenvector associated with the smallest eigenvalue, $p$, is displayed in green highlighting its positive and negative values. The dashed line represents the perturbed damage field $\alpha^* + p$ used as initial guess for the state transition, the solid black line displays the solution $\alpha_t$ returned after first-order convergence.</figcaption></figure>
anticipated by our LEM strategy, the system would have reached a configuration with two bulk cracks and one boundary crack, as shown in Fig. 2-right (branch $n = 5$).

With continued loading, the damage increases and the system evolves along the $n = 4$ branch beyond the loading value at which stability is lost, as shown in Fig. 4-right. A dissipative branch switching event takes place from the current branch ($n = 4$) at point $B'$ to the branch with $n = 6$; see Fig. 4(right) and the inset $C$ for the corresponding damage profile. The delay in bifurcation results in a branch selection event different from the one in our LEM protocol. The delays at bifurcation for the first and second branch switching events leads to a completely different evolutionary path compared to the LEM protocol until it finally reconnects with the homogeneous branch at point $C'$.

For the compliant substrate model, the solutions of our numerical experiments obtained from the quasi-Newton algorithm (black thick line) are overlaid in Fig. 5(left) onto the equilibrium map, while their stability is depicted in Fig. 5-bottom. Notably, we observe a delay, albeit less pronounced compared to the stiff substrate model, for the bifurcation from the trivial branch which occurs at a value higher than the analytically determined critical load $\bar{\epsilon}_1$. Interestingly, the quasi-Newton method switches to the branch with $n = 2$, which exhibits lower energy compared to the branch with $n = 3$. This occurs despite the delay in bifurcation because the branch $n = 3$ remains accessible for the current value $\bar{\epsilon}$. It's worth noting that the LEM protocol also anticipated a first branch selection event leading to the branch with $n = 2$ with two boundary cracks. We again observe that, even before the first branch-changing event, the quasi-Newton solution deviates from the homogeneous branch, and the damage profile has is perturbed by an oscillatory term of the form $\cos(n\pi x)$, reminiscent of the eigenmode $\alpha_n$ with $n = 2$. These perturbations of the equilibrium profiles are negligible from the global energetic standpoint, as it can be inferred in Figure 5(a), by the superposition between the energy of the computed evolution and the exact total energy of the homogeneous solution, yet have a role in the selection of the next branch as will be evidenced in the following.

The system moves along the branch with $n = 2$ as the load increases, until the stability transition at load $\bar{\epsilon}_2$, indicated by point $A'$. Energy minimization brings the system to the branch $n = 4$ although the branch with $n = 3$ is also accessible and has lower energy. Figure 5(right), displays the corresponding damage profiles before and after the branch switching events at points $B'$ and $C$. The final transition before the reconnection to the homogeneous branch takes place at point $C'$, which is also an unstable configuration. The transition brings the system to the branch $n = 6$ as seen in Fig. 5(right), the corresponding stable damage profile at point D is shown in Fig. 5(right).

In summary, while the results of the quasi-Newton minimization simulations exhibit overall similarities in both models, the trajectory taken by the system in the compliant model more closely adheres with the LEM protocol outlined above. This alignment can be attributed to the specific energy landscape of the compliant model. In case of indeterminacy of the trajectory, such as when multiple stable solutions exist at a loss of stability, how to identify an evolution direction?

## 4.3 Hybrid kick algorithm for branch switching

Our investigation shows that the solutions generated by the quasi-Newton algorithms deviate from the LEM protocol outlined in the preceding sections. We observe delayed bifurcations and the persistence of unstable solutions, both affecting branch selection events along the evolution. This is indeed related to the fact that the energy landscape is already flat when the determinant Hessian gets close to zero. Testing the energy expansion (2) along a the direction of a bifurcation mode, that is setting $y - y_t = hp_n$, where $p_n := (v_n, \alpha_n)(x)$ is the $n$-th eigenmode associated with the Hessian and $0 < h \in \mathbb{R}$, we obtain that, for admissible states $y$ $\delta$-close to the equilibrium $y_t$ we can write the following lower bound

$$ \Psi(y) - \Psi(y_t) \geq \frac{h^2}{2} \lambda_t \|p_0\|^2, \quad \forall y : \|y - y_t\| \leq \delta, $$

where $\lambda_t$ indicates the smallest eigenvalue and $w_0$ the associated eigenmode. The first order term $\delta\Psi(y_t)(v_n, \alpha_n)$ vanishes identically owing to the fact that bifurcation modes are admissible fields for the equilibrium condition. When the smallest eigenvalue continuously approaches zero, the energy landscape morphs from being locally flat (at first order) and locally convex in all admissible directions including, in particular, the directions associated with the eigenmodes, to losing local convexity in the one non-trivial direction associated with the eigenvalue that has changed sign. In the numerical practice, the convergence of quasi-Newton algorithms depends on the construction of an (approximate) strictly positive Hessian matrix. As a consequence, such approximation systematically rules out the ability to capture the change of sign (of the smallest eigenvalue) of the Hessian, that is its singularity, thus the onset of instability and
the instability mode. This can justify the systematic (and algorithm-dependent) delay of the bifurcation events via the quasi-Newton solver, observed in Fig. 4 and 5.

To address this challenge, we introduce a hybrid approach that explicitly takes into account the singular mode of the Hessian. In this method, we continuously monitor the smallest eigenvalue of the complete Hessian matrix for the equilibrium solutions obtained from the quasi-Newton algorithm.

When this eigenvalue significantly diminishes, indicating potential instability, first, we perform a full Newton-Raphson refinement using the solution returned from the quasi-Newton algorithm as an initial guess and obtain the fully homogeneous solution. Then we calculate the corresponding eigenvector $\mathbf{p}$, normalised such that $||\mathbf{p}|| = 1$. We then use this eigenvector to perturb the current solution $\mathbf{X}^*$. This perturbation sets the initial guess for the next minimization step in the quasi-Newton algorithm as $\tilde{\mathbf{X}}^{(0)} = \mathbf{X}^* + \eta\mathbf{p}$, where $\eta$ is the step size. This step size can be determined through a line-search algorithm by minimizing the one-dimensional energy slices along the energy descent mode given by the function

$$f(\eta) = \Psi(\mathbf{X} + \eta\mathbf{p}), \qquad (35)$$

to find an optimal value $\eta$, such that $\eta \in \arg\min_{\eta} f$. We then run the quasi-Newton step to obtain a new critical point $\mathbf{X}'$.

The outcome of this hybrid algorithm is the three anticipated branch-switching events identified for the first model previously reported in Fig. 2. For each of these three cases, the snapshots of the damage fields, the associated instability mode $\mathbf{p}$, the perturbed states $\tilde{\mathbf{X}}^{(0)}$, and the converged solutions $\mathbf{X}'$ are shown in Fig. 6 in the right column.

When the quasi-Newton algorithm is initiated using the perturbed state $\tilde{\mathbf{X}}^{(0)}$ as the initial guess, a state transition occurs which leads the system into a qualitatively different state, cf. the final converged solutions returned from the quasi-Newton minimization are shown in Figs. 6(a-c). For the stiff substrate model, the instability mode which corresponds to the bifurcation from the trivial branch has the form $\cos(n\pi x)$ with $n = 3$, as identified through the linear stability analysis (see Fig. 6(a)). The perturbed quasi-Newton minimization returns a state with two boundary cracks and one interior crack, as shown in Fig. 6(a), which corresponds to the solution found using the arc-length continuation algorithm on the branch with $n = 3$ (see Fig. 2(b), inset A). The hybrid algorithm is also successful in capturing the next two anticipated branch switching events: from the branch with $n = 3$ to the branch with $n = 4$ (see Fig. 6(b)) and from branch with $n = 4$ to the branch with $n = 5$ (see Fig. 6(c)). The hybrid algorithm is able to capture the unique path dictated by the LEM protocol shown in Fig. 2.

When the quasi-Newton algorithm is initiated using the perturbed state $\tilde{\mathbf{X}}^{(0)}$ as the initial guess, a state transition occurs which leads the system into a qualitatively different state, cf. the final converged solutions returned from the quasi-Newton minimization are shown in Figs. 6. For the stiff substrate model, the instability mode which corresponds to the bifurcation from the trivial branch has the form $\cos(n\pi x)$ with $n = 3$, as identified through the linear stability analysis (see Fig. 6(a)). The perturbed quasi-Newton minimization returns a state with two boundary cracks and one interior crack, as shown in Fig. 6(a), which corresponds to the solution found using the arc-length continuation algorithm on the branch with $n = 3$ (see Fig. 2(b), inset A). The hybrid algorithm is also successful in capturing the next two anticipated branch switching events: from the branch with $n = 3$ to the branch with $n = 4$ (see Fig. 6(b)) and from branch with $n = 4$ to the branch with $n = 5$ (see Fig. 6(c)). The hybrid algorithm is able to capture the unique path dictated by the LEM protocol shown in Fig. 2. We also show the evolution of the smallest eigenvalue during all loading process in Fig. 7-left where it is clear that the solution (and the evolution as a whole) is stable.

Let's now focus on the compliant substrate model. In this case, according to the LEM protocol, the path through which the system may evolve under quasi-static loading is not unique due to the multiplicity of (stable) solutions at critical loads. Note that this is not the case for the first model where a single equilibrium branch is accessible to the system for all given critical loads. Despite the quasi-Newton algorithm was able to reproduce one of the possible paths predicted by the LEM protocol for the compliant substrate model, the branch switching events are delayed with respect to the anticipated bifurcations.

The system's response under the hybrid algorithm is depicted in Fig. 6(d-g). The first column displays the perturbed damage profiles, while the second column shows the outcomes of the energy minimization process starting from this perturbed state. A stable state with two nucleated boundary cracks (branch $n = 2$) is achieved following the destabilization of the trivial branch, as illustrated in Fig. 6(d). Subsequently, an interior crack nucleates in the middle of the system (branch $n = 4$), as shown in Fig. 6(e). The third event involves a transition from branch $n = 3$ to branch $n = 5$. The event depicted in Fig. 6(f) is the transition to branch $n = 6$. The last event before reconnecting to the trivial
<figure><img src="image_7.png" /><figcaption>Figure 7: Quasi-static loading simulations with L-BFGS: the energy difference $\Delta\Psi$ (left) and $\Delta\tilde{\Psi}$ (right), between the quasi-Newton solutions and the homogeneous solutions are superimposed onto equilibrium branches. The bottom plot shows the smallest eigenvalue of the second variation as a function of the loading parameter $\bar{\epsilon}$. The positivity of the smallest eigenvalue indicates that the solution (and the evolution as a whole) is stable.</figcaption></figure>

branch is the transition from branch $n = 6$ to branch $n = 7$, the only branch available to the system at this load, see in Fig. 6(f). To summarize, the branches selected by the hybrid algorithm follow the sequence $0 \to 2 \to 3 \to 5 \to 6 \to 0$, which differs from the sequence $0 \to 2 \to 3 \to 6 \to 0$ followed by the quasi-Newton algorithm.

Our numerical findings suggest that equilibria close to the Hessian degeneracy (i.e. with small but positive eigenvalues) are robust in the sense that they are insensitive to perturbations. This is highlighted by the fact that running the quasi-Newton solver on a perturbed state close to the Hessian sign transition yields, after convergence, the same equilibrium state prior to perturbation. On the other hand, when the smallest eigenvalue of the system's Hessian is small and negative, the step size $\eta$ obtained by minimization of (35) is sufficient for the quasi-Newton algorithm to effectively escape the flat region of the energy landscape.

## 4.4 Irreversible evolutions and stability

To highlight the role of irreversibility constraints on the determination of the evolution and its consequences on the space of perturbations, we discuss the full stability problem inequality (27) in our simplified one-dimensional setting through an illustrative numerical computation.

To solve the second-order cone-constrained inequality (27), we employ a numerical method based on the orthogonal decomposition of the Hilbert space $X_0$ according to two mutually polar cones $K_0^+$ and $K^*$, cf. [55]. Given an element $z$ in $X_0$, it can be shown that there exists a unique decomposition into two orthogonal components, $x \in K_0^+$ and $y \in K^*$, where $x$ and $y$ are the closest points in $K_0^+$ and $K^*$ to $z$, respectively.

This decomposition allows one to project the problem into the cone and ensure that the eigen-solution satisfies the constraints imposed by irreversibility. This approach is particularly useful in mechanics and physics when dealing with unilateral constraints or problems where the solution space is naturally bounded by physical considerations (e.g., non-negative stress, plastic deformations, etc.)

More in particular, we implement a simple iterative Scaling-and-Projection algorithm [58] which depends upon one numerical parameter, a scaling factor $\eta > 0$. Given a convex cone $K$ and an initial guess $z_0$ (not necessarily in the convex cone), the algorithm operates by first projecting the vector in the cone, $x^{(k=0)} = P_K(z_0)$. After the projection step, the current estimate for the eigenvalue can be computed using the Rayleigh quotient

$$ \lambda^{(k)} = \frac{x^{(k)}^T H x^{(k)}}{\|x^{(k)}\|}, \qquad (36) $$

where $H$ is the Hessian operator. Then, we compute the residual vector $y^{(k)} = Hx^{(k)} - \lambda^{(k)}x^{(k)}$ and obtain the next iterate $x^{(k+1)} = v^{(k)} / \|v^{(k)}\|$ where $v^{(k)} = P_K(x^{(k)} + \eta y^{(k)})$. The algorithm is repeated
until convergence is achieved. Note that, in the cone-constrained case, the residual vector $y^{(k)}$ need not be zero at convergence.

The irreversible evolution is computed using the same parameters as in the experiment in Figure 1(a) for the rigid substrate model. In Fig. 8 we plot the profiles of the fields at the loading step corresponding with the bifurcation point of the reversible system, namely $H$ in Fig. 2.

The homogeneous damage field $\alpha_h$ and the profiles of the damage eigenfunctions associated with the minimal eigenvalue for the second order bifurcation problem inequality (25) (left panel) coincide with the damage solution and eigenmode ($n = 3$) computed by Fourier series, cf. Fig 6(a). On the other hand, the eigenmode for the full nonlinear stability problem inequality (27) is plotted in the right panel. Note that the profile of the (non-negative) instability modes do not reduce to a mere truncation (i.e. a projection on the cone $K_0^+$) of the bifurcation modes, and this is due - in general - to the necessary regularity of the solutions requiring continuity of first derivatives.

In Fig. 9, we plot the inf-eigenvalue of the Hessian operators for the bifurcation problem (with blue circles) as well as for the (irreversible) stability problem (with a thick orange line). Note that the bifurcation spectrum is singular at $\bar{\epsilon}_t = 0$ following the fact that, in our model, the damage criterion is attained as soon as the load is non-zero, thus the space of admissible state perturbations changes suddenly from $H_0^1(0,1) \times \emptyset$ at $\bar{\epsilon}_t = 0$ to the full space $X_0$ for $\bar{\epsilon}_t > 0$ which includes all (sufficiently smooth) damage perturbations. Conversely, the set of admissible damage perturbations for the stability problem changes from the empty set for the loading interval $0 < \bar{\epsilon}_t \le \bar{\epsilon}_b$ set to the solid cone of smooth positive functions $\{\beta \in H^1((0,1)) : \beta \ge 0\}$ only beyond the bifurcation load for $\bar{\epsilon}_t \ge \bar{\epsilon}_b$. Despite the occurrence of negative eigenvalues for the bifurcation problem, the eigenvalues of the stability problem are all positive, which is a sufficient condition to ascertain the stability (and thus, the observability) of the computed evolution.

The consequences of irreversibility manifest both as a pointwise constraint, ruling out transitions between branches that would require a local decrease of damage in favour of a global decrease of energy (e.g. between $n = 3$ to $n = 4$, see profiles in the Fig. 6(left) of the transition (a)$\mapsto$(b)). On the other hand, irreversibility has a major role as a global constraint on the perturbation space, by introducing a strong nonlinearity. For the current choice of parameters, this is evident in the fact that the computed evolution is qualitatively different from the unconstrained case. Despite allowing potential bifurcations, the system navigates a purely homogeneous branch, whose inf-eigenvalues that are positive indicating a sufficient condition for stability (hence observability) of the computed trajectory.

<figure><img src="image_8.png" /><figcaption>Figure 8: Profiles of the damage field (blue) and of the inf-eigenvectors (green) at the bifurcation load $\bar{\epsilon}_b$, for the bifurcation problem (left) and the stability problem (right). At bifurcation, the damage is constant (homogeneous). The bifurcation mode is the eigenmode on the branch $n = 4$, cf. 2-right. The inf-eigenvector for the stability problem (right, green) corresponds to $\lambda_b = 2 \cdot 10^{-2}$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 9: Crack nucleation should not occur in this scenario. Lower bound of the Hessian spectra for the bifurcation problem (orange line) and the stability problem (blue circles), along the homogeneous evolution. The bifurcation spectrum indicates that the evolution path is unique until $\bar{\epsilon}_b$. Then, despite multiple solution exists, the homogeneous solution is stable. Notice that the observable $\lambda_t$ (inf-eig cone) for stability is discontinuous.</figcaption></figure>

# 5 Discussion

Capturing branch switching phenomena across stability transitions is not an automatic feature of approximate numerical methods. If they rely on approximate information about the Hessian of the energy functional, these methods do not guarantee to systematically detect transitions between critical equilibrium states, when stability is lost. Indeed, this requires a careful determination of the zero-eigenmodes that render singular the exact nonlinear Hessian, which is typically not available in general purpose first order numerical algorithms. In practice, without such information, critical loads for equilibrium transitions become algorithm-dependent and are not consistent with closed-form solutions of the exact evolution problem.

Our evidence is that certain numerical methods can introduce non-physical artifacts which should be distinguished from genuine physical phenomena. Our ongoing work aims to refine numerical techniques to provide more reliable algorithms for analysing irreversible processes in variational evolutionary problems with multiple local minima and a high number of degrees of freedom.

From our numerical experience, first order solutions to strongly nonlinear, nonconvex, and singular problems, like those of interest in the applications exhibit strong sensitivity to numerical errors, possibly leading to spurious bifurcations and artificial state transitions. On the other hand, solutions which inegrate second order information are robust and their observability can be fully characterised.

More than numerical perturbations (which can always arise,) the use of numerical methods relying only on (conjugate) gradients (in lieu of exact Hessians) is prone to introducing non-physical crack nucleation.

This is an important observation, which highlights the need for a thorough investigation of the stability of solutions. If only physical factors are considered, an energetic selection mechanism is already encapsulated in the stability statement in the evolution law. As a consequence, equilibrium solutions under increasing load should be maintained as observable if stable, assuming that no nucleation should occur otherwise.

We present two main options for discussing these considerations in view of the importance of using accurate and robust algorithms in real scenarios:

1. Ignore the Numerical Artifact: focussing solely on first order considerations and acknowledging that the observed computed nucleations may be purely numerical and should not be considered in
physical terms.

2. Highlight the Numerical Artifact: Alternatively, emphasizing that homogeneous solutions should be observable, despite the sensitivity to numerical parameters (artifacts, in the quasi-Newton approach) and the abundance of admissible solutions (in the nonconvex scenario), or

3. Otherwise.

Suggesting that state transitions in complex scenarios should be carefully interpreted, the connection between observability and stability is functional to understanding real patterns that emerge, e.g., in higher dimensions or in other physical systems.

In either case, our computations show that, unless second order analysis is performed, observed nucleations are not necessarily indicative of physical cracks but rather of an interplay between purely physical phenomena, inherent to the nature of natural processes, and numerical biases inherent to the computational methods employed. This distinction is crucial for understanding the limitations and proper application of numerical techniques as predictive tools in contexts where cracks are a real concern for structures.

Future work will include a more detailed exploration of evolutionary algorithms and their implementation for stability analysis of fracture in thin films. Some notable instances are craquelures in artistic paintings [27, 12, 13], brittle stability of the cryosphere [80, 76, 74, 54] and crack-pattern selection in metallic thin films [25]. As understood in this work, the final crack patterns depend heavily on the form of the selected unstable modes. In turn, the irreversibility constraint, primarily implemented to prevent crack closure, directly impacts the stability of solutions and mode selection at transition loads. Therefore, our findings suggest that the success of the phase-field model in predicting crack patterns in thin films relies on the use of robust numerical algorithms.

**Acknowledgements** This work was in part supported by the French National Research Agency (ANR, Project No. ANR-20-CE91-0010) and the Austrian Science Fund (FWF, Project No. I 4913-N) within the framework of the project: Nanoarchitected films for unbreakable flexible electronics (NanoFilm).

## References

[1] Sequential quadratic programming. In Jorge Nocedal and Stephen J Wright, editors, *Numerical Optimization*, pages 526–573. Springer New York, New York, NY, 1999.

[2] R Alessi and F Freddi. Failure and complex crack patterns in hybrid laminates: A phase-field approach. *Composites Part B*, 179:107256, December 2019.

[3] R Baggio, E Arbib, P Biscari, S Conti, L Truskinovsky, G Zanzotto, and O U Salman. Landau-Type theory of planar crystal plasticity. *Phys. Rev. Lett.*, 123(20):205501, November 2019.

[4] R Baggio, O U Salman, and L Truskinovsky. Homogeneous nucleation of dislocations as a pattern formation phenomenon. *Eur. J. Mech. A: Solids*, 99(104897):104897, May 2023.

[5] R Baggio, O U Salman, and L Truskinovsky. Inelastic rotations and pseudoturbulent plastic avalanches in crystals. *Phys Rev E*, 107(2-2):025004, February 2023.

[6] A A L Baldelli and J F Babadjian. A variational model for fracture and debonding of thin films under in-plane loadings. *J. Mech. Phys. Solids*, 2014.

[7] A A L Baldelli and C Maurini. Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. *J. Mech. Phys. Solids*, 152(104424), 2021.

[8] Z Bažant and L Cedolin. Stability of structures: Elastic, inelastic, fracture and damage theories. *Journal of Structural Engineering-asce*, 119:1001–1002, August 2010.

[9] K Bertoldi and M C Boyce. Wave propagation and instabilities in monolithic and periodically structured elastomeric materials undergoing large deformations. *Phys. Rev. B Condens. Matter*, 78(18):184107, November 2008.

[10] R G Bettiol and P Piccione. Instability and bifurcation. *Not. Am. Math. Soc.*, 2020.
[11] S Bochkanov and V Bystritsky. Alglib. Available from: www.alglib.net, 2013.

[12] E. Bosco, A. S. J. Suiker, and N. A. Fleck. Crack channelling mechanisms in brittle coating systems under moisture or temperature gradients. *International Journal of Fracture*, 225(1):1–30, 2020.

[13] E. Bosco, A. S. J. Suiker, and N. A. Fleck. Moisture-induced cracking in a flexural bilayer with application to historical paintings. *Theoretical and Applied Fracture Mechanics*, 112:102779, 2021.

[14] B Bourdin, G A Francfort, and J-J Marigo. Numerical experiments in revisited brittle fracture. *J. Mech. Phys. Solids*, 48(4):797–826, April 2000.

[15] A Carpio and L L Bonilla. Discrete models of dislocations and their motion in cubic crystals. *Phys. Rev. B Condens. Matter*, 71(13):134105, April 2005.

[16] MM Chattaway. The anatomy of bark. vi. peppermints, boxes, ironbarks, and other eucalypts with cracked and furrowed barks. *Australian Journal of Botany*, 3(2):170–176, 1955.

[17] Yang Chen, Dmytro Vasiukov, Lionel Gélébart, and Chung Hae Park. A FFT solver for variational phase-field modeling of brittle fracture. *Comput. Methods Appl. Mech. Eng.*, 349:167–190, June 2019.

[18] J D Clayton and J Knap. A phase field model of deformation twinning: Nonlinear theory and numerical simulations. *Physica D*, 240(9):841–858, April 2011.

[19] Christelle Combescure, Pierre Henry, and Ryan S Elliott. Post-bifurcation and stability of a finitely strained hexagonal honeycomb subjected to equi-biaxial in-plane loading. *Int. J. Solids Struct.*, 88-89:296–318, June 2016.

[20] Sergio Conti and Giovanni Zanzotto. A variational model for reconstructive phase transformations in crystals, and their relation to dislocations and plasticity. *Arch. Ration. Mech. Anal.*, 173(1):69–88, July 2004.

[21] Frank E Curtis and Xiaocun Que. A quasi-newton algorithm for nonconvex, nonsmooth optimization with global convergence guarantees. *Math. Program. Comput.*, 7(4):399–428, December 2015.

[22] Y H Dai and Y Yuan. A nonlinear conjugate gradient method with a strong global convergence property. *SIAM J. Optim.*, 10(1):177–182, January 1999.

[23] 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*, pages 75–101. Springer International Publishing, Cham, 2020.

[24] Eusebius J Doedel. AUTO: A program for the automatic bifurcation analysis of autonomous systems. *Congr. Numer.*, 30(265-284):25-93, 1981.

[25] D Faurie, F Zighem, P Godard, G Parry, T Sadat, D Thiaudière, and P-O Renault. In situ x-ray diffraction analysis of 2D crack patterning in thin films. *Acta Mater.*, 165:177–182, February 2019.

[26] Alphonse Finel, Y Le Bouar, A Gaubert, and U Salman. Phase field methods: Microstructures, mechanical properties and complexity. *C. R. Phys.*, 11(3):245–256, April 2010.

[27] L. Fuster-López, F. C. Izzo, C. K. Andersen, A. Murray, A. Vila, M. Picollo, L. Stefani, R. Jiménez, and E. Aguado-Guardiola. Picasso’s 1917 paint materials and their influence on the condition of four paintings. *SN Applied Sciences*, 2(12):2159, 2020.

[28] P Godard, A Guillot, F Zighem, D Thiaudière, D Faurie, and P O Renault. Strain ratio and thickness effects on plasticity and crack patterns of nickel thin films. *Scr. Mater.*, 213:114638, May 2022.

[29] Lucas Goehring, Rebecca Conroy, Asad Akhter, William J Clegg, and Alexander F Routh. Evolution of mud-crack patterns during repeated drying cycles. *Soft Matter*, 6(15):3562–3567, July 2010.

[30] Bjarne Grimstad and Anders Sandnes. Global optimization with spline constraints: a new branch-and-bound method based on b-splines. *J. Global Optimiz.*, 65(3):401–439, July 2016.
[31] Julien Guénolé, Wolfram G Nöhring, Aviral Vaid, Frédéric Houllé, Zhuocheng Xie, Aruna Prakash, and Erik Bitzek. Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps. *Comput. Mater. Sci.*, 175:109584, April 2020.

[32] Ali Harandi, Shahed Rezaei, Soheil Karimi Aghda, Chaowei Du, Tim Brepols, Gerhard Dehm, Jochen M Schneider, and Stefanie Reese. Numerical and experimental studies on crack nucleation and propagation in thin films. *Int. J. Mech. Sci.*, 258:108568, November 2023.

[33] R Hill. A general theory of uniqueness and stability in elastic-plastic solids. *J. Mech. Phys. Solids*, 6(3):236-249, May 1958.

[34] Tianchen Hu, Johann Guilleminot, and John E Dolbow. A phase-field model of fracture with frictionless contact and random fracture properties: Application to thin-film fracture and soil desiccation. *Comput. Methods Appl. Mech. Eng.*, 368:113106, 2020.

[35] Xin Hu, Yanzhou Ji, Lei Chen, Ricardo A Lebensohn, Long-Qing Chen, and Xiangyang Cui. Spectral phase-field model of deformation twinning and plastic deformation. *Int. J. Plast.*, 143:103019, August 2021.

[36] John W Hutchinson. The role of nonlinear substrate elasticity in the wrinkling of thin films. *Philos. Trans. A Math. Phys. Eng. Sci.*, 371(1993):20120422, June 2013.

[37] Gerard Iooss and Daniel D Joseph. *Elementary stability and bifurcation theory*. Undergraduate Texts in Mathematics. Springer, New York, NY, 2 edition, December 2012.

[38] Mahdi Javanbakht and Valery I Levitas. Phase field approach to dislocation evolution at large strains: Computational aspects. *Int. J. Solids Struct.*, 82:95-110, March 2016.

[39] Mats Kirkesæther Brun, Thomas Wick, Inga Berre, Jan Martin Nordbotten, and Florin Adrian Radu. An iterative staggered scheme for phase field brittle fracture propagation with stabilizing parameters. *Comput. Methods Appl. Mech. Eng.*, 361:112752, April 2020.

[40] Philip K Kristensen and Emilio Martínez-Pañeda. Phase field fracture modelling using quasi-newton methods and a new adaptive step scheme. *Theor. Appl. Fract. Mech.*, 107:102446, June 2020.

[41] Charlotte Kuhn, Alexander Schlüter, and Ralf Müller. On degradation functions in phase field fracture models. *Comput. Mater. Sci.*, 108:374-384, October 2015.

[42] A Kumar, B Bourdin, G A Francfort, and others. Revisiting nucleation in the phase-field approach to brittle fracture. *J. Mech. Phys. Solids*, 142(104027), 2020.

[43] Andrés A León Baldelli and Blaise Bourdin. On the asymptotic derivation of Winkler-Type energies from 3D elasticity. *J. Elast.*, 121(2):275-301, December 2015.

[44] Adrian S Lewis and Michael L Overton. Nonsmooth optimization via quasi-newton methods. *Math. Program.*, 141(1):135-163, October 2013.

[45] Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. *Math. Program.*, 45(1-3):503-528, August 1989.

[46] G R Liu and S S Quek. *The Finite Element Method: A Practical Course*. Butterworth-Heinemann, August 2013.

[47] Z Liu, J Reinoso, and M Paggi. Phase field modeling of brittle fracture in large-deformation solid shells with the efficient quasi-newton solution and global-local approach. *Comput. Methods Appl. Mech. Eng.*, 399:115410, September 2022.

[48] T Lookman, S R Shenoy, K Ø Rasmussen, A Saxena, and A R Bishop. Ferroelastic dynamics and strain compatibility. *Phys. Rev. B: Condens. Matter Mater. Phys.*, 67(2):024114, January 2003.

[49] J.-J. Marigo and G. A. Francfort. Revisiting brittle fracture as an energy minimization problem. *J. Mech. Phys. Sol.*, 35(7):1319-1342, 1998.
[50] Hermann Matthias and Gilbert Strang. The solution of nonlinear finite element equations. *Int. J. Numer. Methods Eng.*, 14(11):1613–1626, January 1979.

[51] Alexander E Mayer, Vasiliy S Krasnikov, and Victor V Pogorelko. Homogeneous nucleation of dislocations in copper: Theory and approximate description based on molecular dynamics and artificial neural networks. *Comput. Mater. Sci.*, 206:111266, April 2022.

[52] C Miehe, F Welschinger, and M Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. *Int. J. Numer. Methods Eng.*, 83(10):1273–1311, 2010.

[53] Christian Miehe, Martina Hofacker, and Fabian Welschinger. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. *Comput. Methods Appl. Mech. Eng.*, 199(45):2765–2778, November 2010.

[54] R. Millan, E. Jager, J. Mouginot, M. H. Wood, S. H. Larsen, P. Mathiot, N. C. Jourdain, and A. Bjørk. Rapid disintegration and weakening of ice shelves in north greenland. *Nature Communications*, 14(1):6914, 2023.

[55] J J Moreau. Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires. *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 255:238-240, 1962.

[56] J Nocedal and S J Wright. Quasi-Newton methods. In Jorge Nocedal and Stephen J Wright, editors, *Numerical Optimization*, pages 135–163. Springer New York, New York, NY, 2006.

[57] Subrahmanyam Pattamatta, Ryan S Elliott, and Ellad B Tadmor. Mapping the stochastic response of nanostructures. *Proc. Natl. Acad. Sci. U. S. A.*, 111(17):E1678–86, April 2014.

[58] A Pinto da Costa and A Seeger. Cone-constrained eigenvalue problems: theory and algorithms. *Comput. Optim. Appl.*, 45(1):25–57, January 2010.

[59] I Plans, A Carpio, and L L Bonilla. Homogeneous nucleation of dislocations as bifurcations in a periodized discrete elasticity model. *EPL*, 81(3):36001, December 2007.

[60] G Puglisi and L Truskinovsky. Thermodynamics of rate-independent plasticity. *J. Mech. Phys. Solids*, 53(3):655–679, March 2005.

[61] Zhao Qin, Nicola M Pugno, and Markus J Buehler. Mechanics of fragmentation of crocodile skin and other thin films. *Sci. Rep.*, 4(1):1–7, May 2014.

[62] D Riccobelli, P Ciarletta, G Vitale, C Maurini, and L Truskinovsky. Elastic instability behind brittle fracture in soft solids. November 2023.

[63] A Ruffini and A Finel. Phase-field model coupling cracks and dislocations at finite strain. *Acta Mater.*, 92:197–208, 2015.

[64] O U Salman. Modeling of spatio-temporal dynamics and patterning mechanisms of martensites by phase-field and lagrangian methods. *Université Pierre er Marie-Curie*, Ph.D. thesis, 2009.

[65] O U Salman and L Truskinovsky. De-localizing brittle fracture. *J. Mech. Phys. Solids*, 154(104517):104517, September 2021.

[66] O U Salman, G Vitale, and L Truskinovsky. Continuum theory of bending-to-stretching transition. *Phys Rev E*, 100(5-1):051001, November 2019.

[67] Oguz Umut Salman, Benson Muite, and Alphonse Finel. Origin of stabilization of macrotwin boundaries in martensites. *Eur. Phys. J. B*, 92(1):20, January 2019.

[68] O.U. Salman and L. Truskinovsky. On the critical nature of plastic flow: One and two dimensional models. *International Journal of Engineering Science*, 59:219–254, 2012.

[69] Conrad Sanderson and Ryan Curtin. Armadillo: a template-based c++ library for linear algebra. *J. Open Source Softw.*, 1(2):26, June 2016.

[70] Chuang-Shi Shen, Chao Zhang, Xiaosheng Gao, and Yulong Li. Fragmentation of shells: an analogy with the crack formation in tree bark. *Philosophical Magazine Letters*, 100(6):294–305, 06 2020.
[71] V D Simone, D Serafino, and G Toraldo. Recent advances in nonlinear optimization and equilibrium problems : a tribute to marco D'Apuzzo. MATHEMATICS NOTEBOOKS, 27, 2012.

[72] Eduard Stiefel. Methods of conjugate gradients for solving linear systems. *J. Res. Natl. Bur. Stand.*, 49:409-435, 1952.

[73] Erlend Storvik, Jakub Wiktor Both, Juan Michael Sargado, Jan Martin Nordbotten, and Florin Adrian Radu. An accelerated staggered scheme for variational phase-field models of brittle fracture. *Comput. Methods Appl. Mech. Eng.*, 381:113822, August 2021.

[74] Sainan Sun and G. Hilmar Gudmundsson. The speedup of pine island ice shelf between 2017 and 2020: reevaluating the importance of ice damage. *Journal of Glaciology*, pages 1-9, 2023.

[75] E B Tadmor, M Ortiz, and R Phillips. Quasicontinuum analysis of defects in solids. *Philos. Mag. A*, 73(6):1529-1563, June 1996.

[76] Jeff Tollefson. Giant crack in antarctic ice shelf spotlights advances in glaciology. *NATURE — NEWS*, 2017.

[77] Claus Othmar Wolfgang Trost, Alice Lassnig, Patrice Kreiml, Tanja Jörg, Velislava L Terziyska, Christian Mitterer, and Megan Jo Cordill. Enthalpy-driven self-healing in thin metallic films on flexible substrates. *Adv. Mater.*, May 2024.

[78] L Truskinovsky. Fracture as a phase transition, in “contemporary research in the mechanics and mathematics of materials,” CIMNE, Barcelona, (1996), 322–332. *Received August*, 2010.

[79] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Courapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261-272, 2020.

[80] Jérôme Weiss and Véronique Dansereau. Linking scales in sea ice mechanics. *Philos Trans A Math Phys Eng Sci*, 375(2086), Feb 2017.

[81] Jérôme Weiss, Peng Zhang, Oğuz Umut Salman, Gang Liu, and Lev Truskinovsky. Fluctuations in crystalline plasticity. *C. R. Phys.*, 22(S3):1-37, March 2021.

[82] Thomas Wick. Modified newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation. *Comput. Methods Appl. Mech. Eng.*, 325:577-611, October 2017.

[83] Jian-Ying Wu, Yuli Huang, and Vinh Phu Nguyen. On the BFGS monolithic algorithm for the unified phase field damage theory. *Comput. Methods Appl. Mech. Eng.*, 360:112704, 2020.

[84] Chengxian Xu and Jianzhong Zhang. A survey of Quasi-Newton equations and Quasi-Newton methods for optimization. *Ann. Oper. Res.*, 103(1):213-234, March 2001.

[85] J Yang, J Duan, Y J Wang, and M Q Jiang. Complexity of plastic instability in amorphous solids: Insights from spatiotemporal evolution of vibrational modes. *Eur. Phys. J. E Soft Matter*, 43(9):56, September 2020.

[86] P Zhang, O U Salman, J Weiss, and L Truskinovsky. Variety of scaling behaviors in nanocrystalline plasticity. *Phys Rev E*, 102(2-1):023006, August 2020.

[87] Camilla Zolesi and Corrado Maurini. Stability and crack nucleation in variational phase-field models of fracture: Effects of length-scales and stress multi-axiality. *Journal of the Mechanics and Physics of Solids*, 192:105802, 2024.