form of the degradation function. Nonetheless having $\phi_c > 0$ allows the stress paths for subsequent unloading and reloading to coincide with the initial loading curve, which cannot be achieved otherwise. For the present work, we have chosen to set $\phi_c = 0.5$ in all of our simulations.

The coupled system described in is implemented in a classical finite element framework, with the primary unknowns being the displacement $\bm{u}$ and phase-field $\phi$. In a 2-dimensional setting, these are expressed in terms of the corresponding nodal degrees of freedom as $$\bm{u} = \sum\limits_{I=1}^m \bm{N}_I^{\bm{u}} \bm{u}_I \quad\text{and}\quad \phi = \sum\limits_{I=1}^m N_I \phi_I$$ wherein $$\bm{N}_I =
    \left[ \begin{matrix}
        N_I & 0 \\[0.5em] 0 & N_I
    \end{matrix} \right]$$ with $N_I = N_I \left( \bm{x} \right)$ denoting the shape function associated with node $I$, and $\bm{u}_I$ and $\phi_I$ the respective displacement and phase-field degrees of freedom at node $I$. The strain and phase-field gradient are given by $$\bm{\epsilon} = \sum\limits_{I=1}^m \bm{B}_I^{\bm{u}} \bm{u}_I \quad\text{and}\quad \nabla\phi = \sum\limits_{I=1}^m \bm{B}^\phi_I \phi_I$$ in which $$\bm{B}_I^{\bm{u}} =
    \left[ \begin{matrix}
        N_{I,x} & 0 \\[0.5em] 0 & N_{I,y} \\[0.5em] N_{I,y} & N_{I,x}
    \end{matrix} \right]
    \quad\text{and}\quad
    \bm{B}_I^\phi =
    \left[ \begin{matrix}
        N_{I,x} \\[0.5em] N_{I,y}
    \end{matrix} \right].$$ The former is the symmetrized gradient matrix associated 