The "surface" of the SDF is given as the set $S = f^{-1}[0] =\{p\in\mathbb{R}^n\,:\, f(p) =0 \}$. Thus you can generate a number of points $p_1,\ldots,p_m\in\mathbb{R}^n$ and set them as initial guesses to a Newton iteration:
\begin{align}
p^{k+1}_i &= p^k_i - \theta^k_i (J_f(p^k_i))^+f(p^k_i) \\
&= p_i^k -\theta^k_i(J_f(p_i^k)^TJ_f(p_i^k))^{-1}J_f(p^k_i)^Tf(p^k_i)\\
&= p^k_i + \underbrace{\frac{f(p_i^k)\theta^k_i}{\|\nabla f(p^k_i)\|_2}}_{\alpha^k_i}\underbrace{\left(-\frac{\nabla f(p^k_i)}{\|\nabla f(p^k_i)\|_2}\right)}_{d^k_i}\\
&= p^k_i +\alpha^k_i d^k_i. 
\end{align}

For an actual sdf $\|\nabla f(p^k_i)\|_2 = 1$, and if you set the damping factor $\theta^k_i=1$, then you get $$p^{k+1}_i = p^k_i - f(p^k_i)\nabla f(p^k_i).$$

Note that $f(p^k_i)$ is the signed distance to the surface. So you are doing sphere tracing but with the direction from the negated gradient.

Another option, if you have a point $p$ alread on the surface, is to construct an orthonormal basis for the space orthogonal to $\nabla f(p)$. The latter would be the tangent space. And then take steps along each axis of this space and subsequently project onto the surface. The size of the steps could use information from the Hessian in order to figure out the curvature in each direction. This is more in line with the idea of tracing out the surface starting from a point, while the former is based on arbitrary points traveling down to the surface along the gradient field induced by the sdf. On that note, sometimes it may be beneficial to take $\theta_k$ larger and sometimes smaller in order to boost convergence speed. 

One could equivalently look at this problem as $$\min_{p}\frac{1}{2}\|f(p)\|^2_2.$$
The gradient of the above is precisely $f(p)\nabla f(p)$, so you could just as well interpret this as a gradient descent for the above energy.

If you don't have an sdf but a function that underestimates the distance then I would use $$p^{k+1}_i = p^k_i - \theta^k_i g(p^k_i) \frac{\nabla g(p^k_i)}{\|\nabla g(p^k_i)\|_2}.$$ The convergence would be slower depending on how much $g$ underestimates the real distance. To accelerate both methods you can use something like a heavy-ball method, or also look at (nonlinear) conjugate gradients to handle more unpleasant sdfs. If your sdf is non-differentiable then maybe subgradient descent may be useful.