\frac{1}{M^\ast}\tilde{O}-\frac{1}{8{M^\ast}^2}[\tilde{O},[\tilde{O},\tilde{\varepsilon}]] +\frac{1}{8}\beta[\tilde{O},\mathcal{B}]\nonumber\\
&-\beta \tilde{O}^2\frac{1}{8{M^\ast}^3}\tilde{O}^2
-\frac{1}{8M^\ast}\mathcal{B}[\tilde{O},\tilde{\varepsilon}],\\
\tilde{O}' = \,&\frac{\beta}{2M^\ast}[\tilde{O},\tilde{\varepsilon}]-\frac{1}{24}\Big(4\tilde{O}^2\frac{1}{{M^\ast}^2}\tilde{O} +4\tilde{O}\frac{1}{{M^\ast}^2}\tilde{O}^2
          \nonumber\\
&+3\tilde{O}\frac{1}{{M^\ast}}[\tilde{O},\mathcal{B}]+
3[\tilde{O},\mathcal{B}]\frac{1}{{M^\ast}}\tilde{O}\Big)\nonumber\\
&-\frac{1}{48{M^\ast}^3}\beta[\tilde{O},[\tilde{O},[\tilde{O},\tilde{\varepsilon}]]],\\
    \tilde{\Lambda}' = \,&-\frac{i\beta}{4}\left(\frac{1}{M^\ast}\tilde{O}'+\tilde{O}'\frac{1}{M^\ast}\right).\end{aligned}$$

The corresponding reconstituted FW transformation reads $$\begin{aligned}
H''= e^{i\tilde{\Lambda}'} H' e^{-i\tilde{\Lambda}'}.\end{aligned}$$

With the recursions [@Guo2019Chin.Phys.C43.114105], the Hamiltonian of the reconstituted FW transformation eventually reads $$\begin{aligned}
\label{FWH}
\mathcal{H}_{\rm rFW}=\,&
\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}]\nonumber\\
&-\frac{1}{8{M^\ast}^2}[\tilde{O},[\tilde{O},\tilde{\varepsilon}]]
-\frac{1}{8M^\ast}\mathcal{B}[\tilde{O},\tilde{\varepsilon}]\nonumber\\
&-\frac{1}{8}\beta \tilde{O}^2\frac{1}{{M^\ast}^3}\tilde{O}^2
-\frac{1}{8{M^\ast}^3}\beta[\tilde{O},\tilde{\varepsilon}][\tilde{O},\tilde{\varepsilon}],\end{aligned}$$ and its off-diagonal parts have been raised up to at least the $1/{M^\ast}^4$ order.

For the systems with the spherical symmetry, the corresponding radial single-nucleon Dirac equation reads [@Meng2006Prog.Part.Nucl.Phys.57.470--563; @Liang2015Phys.Rep.570.1--84] $$\label{eq:Dirac1}
\left(
\begin{array}{cc}
\Sigma(r)+M & -\frac{\textrm{d}}{\textrm{d}r}+\frac{\kappa}{r}\\
\frac{\textrm{d}}{\textrm{d}r}+\frac{\kappa}{r} & \Delta(r)-M
\end{array}
\right )
\left(
\begin{array}{c}
G(r) \\ F(r)
\end{array}
\right )
=E
\left(
\begin{array}{c}
G(r) \\ F(r)
\end{array}
\right ),$$ with the normalization condition $\int \mathrm{d}r\,[G^2(r) + F^2(r)] = 1$, where $\kappa$ is a good quantum number defined as $\kappa=\mp~(j+{1}/{2})$ for $j=l\pm{1}/{2}$, and $\Sigma(r) = V(r) + S(r)$ and $\Delta(r) = V(r) - S(r)$ are the sum of and the difference between the vector and scalar potentials, respectively. The single-particle energies $E = \varepsilon +M $ include the mass of nucleon. The initial conditions read $$\begin{aligned}
\label{initialnew}
\tilde{\varepsilon}={\left( \begin{array}{cc}
V(r) & 0\\
0 & V(r)
\end{array}
\right )}&,\quad
\tilde{O}={\left( \begin{array}{cc}
0 &-\frac{\textrm{d}}{\textrm{d}r}+\frac{\kappa}{r}\\
\frac{\textrm{d}}{\textrm{d}r}+\frac{\kappa}{r}& 0
\end{array}
