$\tilde{\varepsilon}$ satisfy that $[\tilde{\varepsilon}, \beta] = 0$ and $\{\tilde{O}, \beta\}$ = 0, respectively. It is important to notice that $[1/{M^\ast}, \tilde{O}] \neq 0$ here. The Hamiltonian  then reads $$\begin{aligned}
H&=\beta (M+S)+\boldsymbol{\alpha} \cdot \mathbf{p}+V\nonumber\\
 &=\beta M^\ast+\tilde{O}+\tilde{\varepsilon}.\end{aligned}$$

In order to simplify later calculations, some notations are made as follows,

$$\begin{aligned}
\mathcal{A} &\equiv [\tilde{O},M^\ast],\\
\mathcal{B} &\equiv [\tilde{O},\frac{1}{M^\ast}],\end{aligned}$$

and it is not difficult to derive that $$\begin{aligned}
\frac{1}{M^\ast}\mathcal{A}+M^\ast\mathcal{B}=0.\end{aligned}$$

Similarly to the operations in the conventional FW transforamtion [@Pryce1948Proc.RoyalSoc.A195.62--81; @Foldy1950Phys.Rev.78.29--36; @Tani1951; @Foldy1952Phys.Rev.87.688--693; @Greiner1990.277--290], the Dirac Hamiltonian is transformed into $$\begin{aligned}
H'=\,&e^{i\tilde{\Lambda}}He^{-i\tilde{\Lambda}}\nonumber\\
   =\,&H+i[\tilde{\Lambda},H]+\frac{i^2}{2!}[\tilde{\Lambda},[\tilde{\Lambda},H]]+\cdots\nonumber\\
      &+\frac{i^n}{n!}[\underbrace{\tilde{\Lambda},[\tilde{\Lambda},\cdots,[\tilde{\Lambda}}_n,H]\cdots]]+\cdots\nonumber\\
   =\,&H'_0+H'_1+\cdots+H'_n+\cdots\end{aligned}$$ with its $n$-th component $$\begin{aligned}
H'_n
=\frac{i^n}{n!}[\underbrace{\tilde{\Lambda},[\tilde{\Lambda},\cdots,[\tilde{\Lambda}}_n,H]\cdots]].\end{aligned}$$

Keeping all the terms up to the $1/{M^\ast}^3$ order, the unitary transformed Hamiltonian reads $$\begin{aligned}
\label{firstH}
H'
=\,&\beta M^\ast +\tilde{\varepsilon}+\frac{1}{2}\beta \tilde{O} \frac{1}{M^\ast}\tilde{O}+\frac{1}{8}\beta[\tilde{O},\mathcal{B}] -\frac{1}{8{M^\ast}^2}[\tilde{O},[\tilde{O},\tilde{\varepsilon}]]
\nonumber\\
&-\frac{1}{8M^\ast}\mathcal{B}[\tilde{O},\tilde{\varepsilon}]
-\frac{1}{8}\beta \tilde{O}^2\frac{1}{{M^\ast}^3}\tilde{O}^2
+\frac{\beta}{2M^\ast}[\tilde{O},\tilde{\varepsilon}]\nonumber\\
&-\frac{1}{24}(4\tilde{O}^2\frac{1}{{M^\ast}^2}\tilde{O}+4\tilde{O}\frac{1}{{M^\ast}^2}\tilde{O}^2 +3\tilde{O}\frac{1}{{M^\ast}}[\tilde{O},\mathcal{B}]\nonumber\\
&+
3[\tilde{O},\mathcal{B}]\frac{1}{{M^\ast}}\tilde{O}) -\frac{1}{48{M^\ast}^3}\beta[\tilde{O},[\tilde{O},[\tilde{O},\tilde{\varepsilon}]]].\end{aligned}$$

In this unitary transformed Hamiltonian (\[firstH\]), the off-diagonal parts are not zero but raised by one order from $\tilde{O}$ to $\frac{1}{2M^\ast}\beta[\tilde{O},\tilde{\varepsilon}]$. In order to make the off-diagonal parts higher than a given order (e.g., the $1/{M^\ast}^3$ order), one should repeat the FW transformation until the accuracy is achieved [@Greiner1990.277--290; @Guo2019Chin.Phys.C43.114105]. For that, the operators $\tilde{O}'$ and $\tilde{\varepsilon}'$ are redefined according to Eq. (\[firstH\]), and one has

$$\begin{aligned}
\tilde{\varepsilon}' = \,&\tilde{\varepsilon}+\frac{1}{2}\beta \tilde{O} 