## 2 Methods
### Numerical Method  
The analysis of adhesive interaction at the HDI requires the use of numerical methods for Equation (11). To improve the computational efficiency, the inexact Newton-Bi-conjugate stabilized method is applied, and the numerical procedure is given in Figure 4. The key techniques involved are introduced in the following. (1) Determination of the equilibrium distance \(z_{01}\):  
At the equilibrium distance, the net interaction force of two surfaces vanishes, which implies that \(P_{\rm vdW}+P_{\rm rep}=0\) at \(h=z_{01}\). Therefore, \(z_{01}\) satisfies  
\[-\frac{kT}{8\pi z_{01}^{3}}\sum_{n=0}^{\infty}\;\;,^{\prime}\int_{r_{n}}^{ \infty}dr\;r^{2}\cdot\left[\frac{\overline{\Delta}_{31}^{E}\overline{\Delta}_ {32}^{E}e^{-r}}{1-\overline{\Delta}_{31}^{E}\overline{\Delta}_{32}^{E}e^{-r} }+\frac{\overline{\Delta}_{31}^{F}\overline{\Delta}_{32}^{E}e^{-r}}{1- \overline{\Delta}_{31}^{F}\overline{\Delta}_{32}^{E}e^{-r}}\right]+\frac{B}{4 5\pi z_{01}^{9}}=0 \tag{17}\]  
The relaxation iterative method is used to solve the above equation. (2) Determination of the vdW pressure \(P_{\rm vdW}\):  
According to Equation (5), the vdW pressure is a two-dimensional summation operation, which is time-consuming in the numerical procedure of the Newton method. To reduce computational burdens, the vdW pressure at different distances is calculated at first and then used through interpolation during the iterative solving of Equation (11). Figure 4: Flowchart for the computational procedure. (3) Determination of the initial estimate of distance \(h\):  
For the Newton method, the initial value of \(h\) can influence the convergence of the iterative procedure. When the nominal flying height \(z_{r}\) is larger compared with the TP height, the adhesive interaction at the HDI is small and the surface deformation can be neglected. Thus, \(h=l\text{sin}\theta+z_{r}-z_{TP}\) can be taken as the initial guess of \(h\). However, with the decrease in the nominal flying height \(z_{r}\), the adhesive interaction becomes remarkable and the initial estimate of \(h\) should take the surface deformation into account. In addition, jump-in and jump-out instabilities [29] may occur at the HDI. Near the jump-in and jump-out points, small changes in the nominal flying height may lead to the abrupt variation of the interface interaction between the head and the disk, which will influence the convergence of the numerical simulation. To solve this problem, the nominal flying height is substituted by the distance \(h(0,\,0)\). Denoting \(h_{\text{c}}=h(0,\,0)\), the following equation is derived from Equation (11). \[h_{c}=h(0,0)=l\sin\theta+z_{r}-z_{TP}(0,0)+u(0,0) \tag{18}\]  
By combining Equation (11) with Equation (18), we have  
\[h(x,y)=h_{c}-[z_{TP}(x,y)-z_{TP}(0,0)]+[u(x,y)-u(0,0)] \tag{19}\]  
Through the above substitution, the distance \(h\) becomes a function of the central distance \(h_{\text{c}}\). To obtain the initial value of \(h\) at a preset nominal flying height \(z_{rest}\), the simulation procedure starts from a large central distance \(h_{\text{c}}\) and solves Equation (19) using an initial value \(h=h_{\text{c}}-[z_{TP}(x,y)-z_{TP}(0,0)]\). After achieving the convergent solution, \(z_{r}\) can be inversely calculated by Equation (18). If \(z_{r}<z_{rest}\), then the central distance is decreased by a preset interval \(\Delta h_{\text{c}}\), and the existing solution at \(h_{\text{c}}\) is used as the initial estimate at \(h_{\text{c}}-\Delta h_{\text{c}}\). The above process is repeated until the condition \(z_{r}<z_{rest}\) does not satisfy. Using the obtained solution as the initial guess at \(z_{r}=z_{rest}\), and solving Equation (11), the desirable convergent solution can be acquired. The detailed flowchart for the computational procedure is presented in Figure 4, and the element of residuals \(R_{r},R_{r}^{\prime}\) and that of the Jacobian Matrix \(J,J^{\prime}\) in discrete form can be written as  
\[