function $\delta(\bm{r}'-\bm{r})$ will be denoted as $\boldsymbol{\Delta}(\bm{r}',\bm{r})$ (with a bold font in order to distinguish from $\Delta=V-S$).

In the present scheme, all operators should be transformed in the same way as the Dirac Hamiltonian. Therefore, the delta function in the FW formalism is transformed as $$\begin{aligned}
\boldsymbol{\Delta}'=e^{i\tilde{\Lambda}}\boldsymbol{\Delta}e^{-i\tilde{\Lambda}}.\end{aligned}$$ The expansion of $\boldsymbol{\Delta}'$ reads $$\begin{aligned}
\boldsymbol{\Delta}'
  =\,&\boldsymbol{\Delta}+i[\tilde{\Lambda},\boldsymbol{\Delta}]+\frac{i^2}{2!}[\tilde{\Lambda},[\tilde{\Lambda},\boldsymbol{\Delta}]]\nonumber\\
   &+\cdots+\frac{i^n}{n!}[\underbrace{\tilde{\Lambda},
   [\tilde{\Lambda},\cdots,[\tilde{\Lambda}}_n,\boldsymbol{\Delta}]\cdots]]+\cdots\end{aligned}$$ with its $n$-th component $$\begin{aligned}
\boldsymbol{\Delta}'_n
=\frac{i^n}{n!}[\underbrace{\tilde{\Lambda},[\tilde{\Lambda},\cdots,[\tilde{\Lambda}}_n,\boldsymbol{\Delta}]\cdots]].\end{aligned}$$ As a result, $$\begin{aligned}
\boldsymbol{\Delta}'=\,&\boldsymbol{\Delta}-\frac{1}{8{M^\ast}^2}\boldsymbol{\Delta}^{(2)}-\frac{1}{8M^\ast}\mathcal{B}\boldsymbol{\Delta}^{(1)}+\frac{\beta}{2M^\ast}\boldsymbol{\Delta}^{(1)}\nonumber\\
&-\frac{\beta}{48{M^\ast}^3}\boldsymbol{\Delta}^{(3)}.\end{aligned}$$ Here the relevant notations are defined as

$$\begin{aligned}
\boldsymbol{\Delta}^{(1)} &\equiv [\tilde{O},\boldsymbol{\Delta}],\\
\boldsymbol{\Delta}^{(2)} &\equiv [\tilde{O},\boldsymbol{\Delta}^{(1)}],\\
\boldsymbol{\Delta}^{(3)} &\equiv [\tilde{O},\boldsymbol{\Delta}^{(2)}].\end{aligned}$$

Following the recursion relation in the reconstituted FW transformation, the result up to the $1/{M^\ast}^3$ order reads $$\begin{aligned}
\boldsymbol{\Delta}_{\rm rFW}
=\,&\boldsymbol{\Delta}-\frac{1}{8{M^\ast}^2}\boldsymbol{\Delta}^{(2)}-\frac{1}{8M^\ast}\mathcal{B}\boldsymbol{\Delta}^{(1)}-\frac{\beta}{4{M^\ast}^3}[\tilde{O},\tilde{\varepsilon}]\boldsymbol{\Delta}^{(1)}\nonumber\\
\,&+\frac{\beta}{2M^\ast}\boldsymbol{\Delta}^{(1)}-\frac{\beta}{2{M^\ast}^3}\boldsymbol{\Delta}^{(1)}\tilde{O}^2-\frac{\beta}{2{M^\ast}^3}\boldsymbol{\Delta}^{(2)}\tilde{O}\nonumber\\
&-\frac{9\beta}{48{M^\ast}^3}\boldsymbol{\Delta}^{(3)}.\end{aligned}$$

The delta function can be expanded in terms of spherical harmonics, $$\begin{aligned}
\delta(\bm{r}'-\bm{r})=\frac{1}{r^2}\delta(r'-r)\sum_{L=0}^{\infty}\bm{Y}_L(\hat{\bm{r}}')\cdot\bm{Y}_L(\hat{\bm{r}}).\end{aligned}$$ When the spherical symmetry is taken into account, it reads $$\begin{aligned}
\delta(\bm{r}'-\bm{r})=\frac{1}{4\pi r^2}\delta(r'-r).\end{aligned}$$ Thus, $$\begin{aligned}
\boldsymbol{\Delta}(r',r)=
\frac{1}{4\pi r^2}
\left( \begin{array}{cc}
\delta(r'-r) & 0\\
0 & \delta(r'-r)
\end{array}
\right ).\end{aligned}$$

In addition, one has

$$\begin{aligned}
\boldsymbol{\Delta}^{(1)}(r',r)=\,&
\frac{1}{4\pi r^2}
\left( \begin{array}{cc}
0 & -\delta'\\
\delta' & 0
\end{array}
\right ),\\
\boldsymbol{\Delta}^{(2)}(r',r)=\,&
\frac{1}{4\pi r^2}
\left( \begin{array}{cc}
-\delta''+\frac{2\kappa}{r'}\delta' & 0\\
0 & -\delta''-\frac{2\kappa}{r'}\delta'
\end{array}
\right ),\\
\boldsymbol{\Delta}^{(3)}(r',r)=\,&
\frac{1}{4\pi r^2}
\left( \begin{array}{cc}
0& \delta'''-\frac{2\kappa}{r'}\delta''+ \frac{2\kappa}{r'^2}\delta'\\
-\delta'''-\frac{2\kappa}{r'}\delta''+ \frac{2\kappa}{r'^2}\delta' & 0
\end{array}
\right ).\end{aligned}$$

