\lefteqn{\mathrm{d}\left(V{}_\mathrm{eff}^\mathrm{int}(r, R, b)\right)/\mathrm{d}r  \equiv}\nonumber \\
      & & \left[2 r R^2 \left(b (Y - 4 R Z) + 6 R^2 Z\right) \left\{-b^3 \left(r^2 (Y-4 R Z)+8 R^2 Y\right) + b^2 R\times \right.\right. \nonumber \\
      & & \times \left(r^2 (3 Y-14 R Z)+R^2 \left(6 R^2 Z - 9 R Y + 40 Y\right)\right)+b R^2 \left(R^2(-13 R^2 Z + \right.\nonumber \\
      & & \left.\left.\left.\left. + 25 R Y- 66 Y\right)- 2 r^2 (Y- 6 R Z)\right)+ 6 R^5 \left(R^2 Z-3 R Y+ 6 Y\right)\right\}\right]\times \nonumber \\
      & &\times \left[9 \left(b+(R-2) R\right) Y \left(R (Y-3 R Z)-b \left(Y-\nonumber 2 R Z\right)\right)^3\right]^{-1} = 0.\\
    \end{aligned}$$ This equation can be solved by using numerical methods and determines the loci $r_\mathrm{c(i)}(R, b)$ of the internal circular null geodesics that are stable and correspond to a local maximum of the internal effective potential. From the radius $r_\mathrm{c(i)}$ we easily obtain the impact parameter $\ell^{2}_{\mathrm{c(i)}}(R, b)$ using equation (\[EQellsb\]). The existence of the internal stable circular null geodesics is allowed only for the surface radius $R$ limited by values explicitly given by $$R_\mathrm{I\pm} = \frac{3}{2} \pm \sqrt{\frac{9}{4} - \frac{5}{3}b}$$ that has to satisfy naturally also the condition $R > R_\mathrm{min}(b)$.

In the external spacetime, the photon circular geodesics and 