## 3 Implicit Backward-Difference Time Integration  
Let \(\Delta t=t^{n+1}-t^{n}>0\), with \(t^{n}\), \(t^{n+1}\in[0-T]\), a time interval along which the constitutive law (1) has to be integrated. A fully implicit, backward-difference scheme is adopted here, whereby all the quantities are evaluated at the end \(t^{n+1}\) of the step, whereas the state of the system is assumed to be completely known at \(t^{n}\). In the following derivations, the isotropic elastic tensor \(\mathbf{D}\) is considered constant in the time step, as is commonly assumed in the literature; see, e.g., [20]. Its value will be updated based on Equation (7) only at the end of the considered time step. An implicit formulation of the algorithm considering also the evolution of \(\mathbf{D}\) within the time step would be possible, at the cost however of additional complications in the formulation of the integration scheme  
whose corresponding residuals are:  
\[\Psi_{1} =\Delta e^{vp}-\left(p_{A}^{\text{eq}}\right)^{n+1}\frac{\Delta e_{v }^{vp}}{\left(p_{p}^{\text{eq}}\right)^{n+1}}=0 \tag{23}\] \[\Psi_{2} =\Delta e_{v}^{vp}-\frac{\mu^{*}}{\tau}\,\Delta t\left[\left(p_{ p}^{\text{eq}}\right)^{n+1}\right]^{\frac{\lambda^{*}-\kappa^{*}}{p^{*}}}=0\] (24) \[\Psi_{3} =\left(p_{p}^{\text{eq}}\right)^{n+1}-p_{p0}^{\text{eq}}\,\exp \left(\frac{e_{v}^{vp}\,n+\Delta e_{v}^{p}}{\lambda^{*}-\kappa^{*}}\right)=0\] (25) \[\Psi_{4} =p^{n+1}-p^{\text{trial}}+K\,\Delta e_{v}^{vp}=0\] (26) \[\Psi_{5} =q^{n+1}-q^{\text{trial}}+3\,G\,\Delta e^{vp}=0 \tag{27}\]  
For convenience, these can be collected in a residual vector \(\mathbf{\
abla}\):  
\[\mathbf{\Psi}_{(5\times 1)}=\left[\Psi_{1}\quad\Psi_{2}\quad\Psi_{3}\quad\Psi_{4} \quad\Psi_{5}\right]^{T} \tag{28}\]  
The trial values in (26) and in (27) are given by:  
\[\sigma^{\text{trial}} =\mathbf{D}:\left(\epsilon^{n+1}-\epsilon^{vp}\,\pi\right),\quad p^{ \text{trial}}=\mathbf{m}:\sigma^{\text{trial}},\quad\mathbf{s}^{\text{trial}}=\mathbf{n}: \sigma^{\text{trial}}, \tag{29}\] \[q^{\text{trial}} =\sqrt{\frac{3}{2}\mathbf{s}^{\text{trial}}:\mathbf{s}^{\text{trial}}}\]  
It should be noted that not all residuals explicitly depend on all unknowns, as can be seen below:  
\[\Psi_{1} =\Psi_{1}\Big{(}\Delta e^{vp},\Delta e_{v}^{vp},p^{n+1},q^{n+1} \Big{)} \tag{30}\] \[\Psi_{2} =\Psi_{2}\Big{(}\Delta e_{v}^{vp},\left(p_{p}^{\text{eq}}\right) ^{n+1},p^{n+1},q^{n+1}\Big{)}\] (31) \[\Psi_{3} =\Psi_{3}\Big{(}\Delta e_{v}^{vp},\left(p_{p}^{\text{eq}}\right) ^{n+1}\Big{)}\] (32) \[\Psi_{4} =\Psi_{4}\Big{(}\Delta e_{v}^{vp},p^{n+1}\Big{)}\] (33) \[\Psi_{5} =\Psi_{5}\Big{(}\Delta e^{vp},q^{n+1}\Big{)} \tag{34}\]  
Proceeding iteratively with a first-order linearization of the resolving system, at the \((i+1)\)-th iteration, one has:  
\[\mathbf{\Psi}_{i+1}\cong\mathbf{\Psi}_{i}+\left.\frac{d\mathbf{\
abla}}{d\mathbf{X}}\right|_{ i}\delta\mathbf{X}=\mathbf{0} \tag{35}\]  
where \(\mathbf{\Psi}_{i}\) defines the residual vector at iteration \(i\) and \(\delta\mathbf{X}\) is the increment of the unknown variables between two subsequent iterations:  
\[\delta\mathbf{X}_{(5\times 1)}=\left[\delta\Delta e^{vp}\quad\delta\Delta e_{v}^{vp} \quad\delta p_{p}^{\text{eq}}\quad\delta p\quad\delta q\right]^{T} \tag{36}\]  
The residual gradient \(d\mathbf{\Psi}/d\mathbf{X}\) is a \(5\times 5\) matrix containing the derivatives of the residuals w.r.t. the unknowns. Its detailed expression is given in Appendix B. It can be shown that the gradient \(d\mathbf{\Psi}/d\mathbf{X}\) is non-symmetric. Consequently, a non-symmetric solver must be used for the local problem, if a consistent tangent matrix is desired to assure an asymptotic quadratic convergence.