## 1 Introduction
In contrast, the CPM only satisfies the consistency condition without requiring the second derivative of the yield function; thus, this approach is also referred to as a semi-explicit method. Regarding the stress integration algorithms on the distortional hardening models, similar approaches, based on either the CPPM or CPM schemes, have been reported. [PERSON] et al. [55] implemented the first version of the HAH (denoted as original HAH hereafter) model into commercial ABAQUS software, using both the CPPM and CPM approaches. This implementation was also extended to the e-HAH model; however, in this case, only the CPM scheme was applied in combination with a sub-stepping numerical method [56]. Recently, [PERSON] and [PERSON] [57] also reported the implementation of the HAH model, using the Euler backward scheme, but they calculated the derivatives of the yield function using numerical finite differences. More recently, [PERSON] et al. [58] coupled the CPPM scheme with a line-search algorithm for an updated version of the HAH model [59] to improve numerical efficiency. However, the existing numerical integration schemes do not fully exploit the implicit CPPM scheme. In other words, these approaches only account for the partial number of residuals in order to attain the solution process of linearized equations under the Newton-Raphson method. In this sense, these methods can be regarded as a partial or semi-implicit integration algorithm based on the Euler backward scheme. In this study, a fully implicit stress integration algorithm based on the Euler backward scheme is proposed for the e-HAH model. To the best of our knowledge, the present work is the first to report this novel scheme. The present algorithmic study differs from previous studies [55; 56; 57; 58] in that the complete set of residuals is established from the stress tensor, microstructure deviator, and other variables controlling the distortion of the yield function. Furthermore, the calculations of the first and second derivatives of the distorting yield function are both analytically expressed. Note that the CPPM or Euler backward scheme presented in previous studies adopted only the consistency condition and flow rule for the residuals. This may lead to non-physical evolutions of the distortional hardening state variables when there is an abrupt change in the loading path. For the comparative study on the effect of the employed stress integration algorithms, four different algorithms were implemented into the FE simulation. These algorithms consist of the CPM and EBM schemes with different numbers of residuals. Again, the EBM scheme is further classified based on either analytical derivatives or finite-difference--based numerical derivatives of the yield function. Comprehensive validations are presented in this study. Firstly, one-element analyses with two different material sets are presented for the investigated stress integration algorithms. The flow stress and anisotropy in deformation were evaluated under two distinctive loading conditions, i.e., reverse and cross-loading paths. Finally, the industrial-size S-rail forming and springback simulations were performed to evaluate the numerical accuracy and computational efficiency of the proposed algorithm. Note that even this large-scale industrial problem was solved using the purely static implicit FE solver, whereas most previous studies employed a dynamic explicit method for the forming process and a static implicit method for the springback, due to the divergence problem resulted from complex contact changes during forming. The manuscript is organized as follows. In Section 2, a summary of the e-HAH model is presented. In Section 3, the numerical integration algorithms investigated in this study are introduced, alongside a summary of previous studies. Then, in Section 4, the accuracy and stability of the individual algorithms are comparatively evaluated via one-element analyses under two distinctive loading path changes.