In our time integration scheme, to satisfy energy balance, this equation needs to be evaluated at $t_{n+\alpha_{f}}$ , thus $$\begin{equation} \left(\dot{\Psi}_{r}\right)_{n+\alpha_{f}}=J_{n+\alpha_{f}}\boldsymbol{\sigma}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}\,.\label{eq:edy-mp-energy} \end{equation}$$ However, the solution for $\boldsymbol{\sigma}_{n+\alpha_{f}}\equiv\boldsymbol{\sigma}\left(\mathbf{F}_{n+\alpha_{f}}\right)$ obtained from the momentum balance may not necessarity satisfy this equation. Thus, to satisfy energy balance over consecutive time steps, we want to evaluate an effective stress $\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}$ such that $$\begin{equation} J_{n+\alpha_{f}}\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}=\frac{\left(\Psi_{r}\right)_{n+1}-\left(\Psi_{r}\right)_{n}}{\Delta t}\,.\label{eq:edy-energy-discretized} \end{equation}$$ To find a solution for $\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}$ , we follow the procedure of Gonzalez [34] and let $$\begin{equation} \tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}=\boldsymbol{\sigma}_{n+\alpha_{f}}+f\mathbf{D}_{n+\alpha_{f}}\,,\label{eq:edy-eff-stress-model} \end{equation}$$ where $f$ is some scalar function to be determined. Substituting this relation, (6.2.6.1-6), into the previous equation, (6.2.6.1-5), produces $$\begin{equation} f=\left(\frac{\left(\Psi_{r}\right)_{n+1}-\left(\Psi_{r}\right)_{n}}{J_{n+\alpha_{f}}^{s}\Delta t}-\boldsymbol{\sigma}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}\right)\frac{1}{\mathbf{D}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}}\,.\label{eq:edy-eff-stress-f} \end{equation}$$ Hence, the equation for an effective stress needed to satisfy energy balance between consecutive time steps is $$\begin{equation} \boxed{\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}=\boldsymbol{\sigma}_{n+\alpha_{f}}+\left(\frac{\left(\Psi_{r}\right)_{n+1}-\left(\Psi_{r}\right)_{n}}{J_{n+\alpha_{f}}\Delta t}-\boldsymbol{\sigma}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}\right)\frac{\mathbf{D}_{n+\alpha_{f}}}{\mathbf{D}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}}}\,.\label{eq:edy-eff-stress-final} \end{equation}$$ In the limit when $\mathbf{D}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}=0$ , we use $\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}=\boldsymbol{\sigma}_{n+\alpha_{f}}$ . Recall that this scheme produces conservation of linear and angular momentum and total energy only with $\rho_{\infty}=1$ , or equivalently, $\alpha_{f}=\alpha_{m}=\frac{1}{2}$ , $\beta=\frac{1}{2}$ and $\gamma=1$ . Therefore, this effective stress calculation is only applied when the user employs $\rho_{\infty}=1$ .