functions $\mathrm{e}^{\Phi(r)}$, $\xi(r)$, and $\phi(r)$ uniquely determines another two functions.

For the sake of brevity, we will explore a simple, fully analytic, one-parameter family of solutions with $$\label{xi}
\xi= \sqrt{r^2+2ar+5a^2}-a,\quad \mathrm{e}^{\Phi}= \frac{r+a}{\sqrt{r^2+2ar+5a^2}}.$$ By direct integration in (\[A-f\]), we obtain $$\begin{gathered}
\label{Aex}
A=1+\frac{2a}{3r}\,-\,2\,\frac{a+3m}{15a} \vphantom{\int\limits_{r}^{\infty}}\\
\times\left\{\!\!\frac{\sqrt{r^2+2ar+5a^2}}{r}\! \left(\!\!1+\frac{r}{a}-\frac{r^2}{a^2}\!\right)+ \frac{r^2}{a^2}\!\right\},\end{gathered}$$ where $a$ is the parameter of ’intensity’ of the scalar field.

In studying bound orbits, we are interested only in the metric functions and can, therefore, use an arbitrary unit of length. On the other hand, the solution (\[A-f\])–(\[V\]) is invariant under the scale transformations $$r\rightarrow r/\lambda,\quad m\rightarrow m/\lambda,\quad V\rightarrow \lambda^2V,\quad\;\;\lambda>0.\nonumber$$ It means that by applying $\lambda=m$ in this transformation, we can take, as it is usually done in general relativity, the mass of a SFNS as the current unit of length. Thus, without loss of generality, we suppose everywhere below that $m=1$.

Using (\[R1\])–(\[R3\]) we find that the Kretchmann invariant diverges at the centre, namely, $$\label{}
K= \frac{\,5 \big(3-\sqrt{5}\,\big) \big(4a-3\sqrt{5}-3\big)^{\!2}} {6\;\!\!r^6\vphantom{\tilde{A}}}\; + \; O\big(r^{-5}\big)\,,\quad r\rightarrow0.
\nonumber$$ In accordance with (\[3m&lt;xi0\]), the condition $a>3/(\sqrt{5}-1)$ determines the subfamily of SFNSs.

The scalar field can be obtained by solving the problem $\phi'=\sqrt{\Phi'/r\,},\,\;\phi(\infty)=0$. The result is $$\label{phi-F}
\phi=\sqrt[4]{4/5}\, \big(F(\lambda_0,k)-F(\lambda(r),k)\big),$$ where $$k=\sqrt{1/2+1/\sqrt{5}}, \quad
\lambda_0= \arcsin\! \left(\frac{2\sqrt[4]{20}}{\,\sqrt{5}+2\,}\right),
\nonumber$$ $$\lambda(r)=\arcsin\! 