Proving the asymptotic normality of th Non-linear Least Squares estimator Our model is $Y=X(\beta_0)+u$, where $u\sim IID(0,\sigma_0^2I)$, and $X(\beta)$ is a non-linear function of the beta.

When trying to minimize the $SSR(\beta)$ we get the following FOC:

$\nabla X(\beta)^T(Y-X(\beta))=0$, where $\nabla X(\beta)$ is the gradient.

Well, the FOC is equivalent to $n^{(-1/2)}(\nabla X(\beta)^T(X(\beta_0)+u-X(\beta))=0$.

If we apply a taylor expansion of the first order to each component $X_t(\beta)$ of $X(\beta)$, we obtain $X_t(\beta)=X_t(\beta_0)+\nabla X(\bar\beta_{(t)})^T(\beta-\beta_0)$, where $\bar\beta_{(t)}$ is a point in the line segment that joins $\beta$ and $\beta_0$. This point may be different for each taylor expansion we do, and that's why it's indexed by $t$.

Inserting the taylor expansion in the FOC:
$n^{(-1/2)}(\nabla X(\beta)^T(u-\nabla \bar X^T(\beta-\beta_0))=0$, where $\nabla \bar X$ is the matrix with $\nabla X(\bar\beta_{(i)})$ as each i-th column.

Are all of the above calculations correct? I ask this because in this book, the authors state in page 225 that we should obtain a term with second derivatives of $X(\beta)$... I do not understand why this is.

Any help would be appreciated
It appears that what you do is not a Taylor expansion, but an application of the mean-value theorem (as one should). If it was a Taylor expansion, apart from the remainder, the gradient would have to be evaluated at $\beta_0$. With the mean-value theorem, there is no remainder, and you evaluate the gradient at some $\bar \beta$ that always lies between $\beta$ and $\hat \beta$.  

As regards the issue of second derivatives/Hessian, officially speaking they/it only "temporarily" appear in the derivation of asymptotic normality of the non-linear least -squares estimator, but they vanish asymptotically (while in the Maximum likelihood estimator the Hessian stays there).  

Apart from that, we want to minimize the sum of squared residuals so, indexing the observations by $i$, and using much simpler notation (you will have to adjust it to vector-matrix notation), you set out to minimize $\sum_i[u_i(\beta)]^2 = \sum_i[y_i-h(\mathbf x_i,\beta)]^2$ w.r. to the vector $\beta$, in order to obtain a $\hat \beta$. So your FOC is (suppressing the regressors and passing the $i$ index to the function $h$)

$$\hat \beta : \sum_i\frac {\partial }{\partial \beta}[y_i-h_i(\beta)]^2 = 0 \implies \sum_i2[y_i-h_i(\hat \beta)] \frac {\partial h_i (\hat \beta)}{\partial \beta} =0, $$

Ignore "$2$" and apply the mean value theorem to the whole expression  to get 

$$\sum_i[y_i-h_i(\beta)] \frac {\partial h_i (\hat \beta)}{\partial \beta} =  \sum_i[y_i-h_i(\beta_0)] \frac {\partial h_i (\beta_0)}{\partial \beta} \\ 
+ (\hat \beta -\beta) \sum_i\left [-\frac {\partial h_i (\bar \beta)}{\partial \beta}\frac {\partial h_i (\bar \beta)}{\partial \beta}+ [y_i-h_i(\bar \beta)]\frac {\partial^2 h_i (\bar \beta)}{\partial \beta^2}\right] =0$$

Note that $y_i-h_i(\beta_0) = u_i$, the true error, divide by $1/n$ and then multiply by $\sqrt n$ as you are allowed to and re-arrange to get 

$$\sqrt n (\hat \beta -\beta) = -\left (\frac 1n\sum_i\left [-\frac {\partial h_i (\bar \beta)}{\partial \beta}\frac {\partial h_i (\bar \beta)}{\partial \beta}\right] + \frac 1n\sum_i\left [[y_i-h_i(\bar \beta)]\frac {\partial^2 h_i (\bar \beta)}{\partial \beta^2}\right]\right)^{-1} \cdot \left(\frac 1{\sqrt n} \sum_iu_i \frac {\partial h_i (\beta_0)}{\partial \beta} \right) $$

Now, we consider asymptotic normality, given that consistency holds,  $\hat \beta \xrightarrow{p} \beta_0$. Since $\bar \beta$ is sandwiched between $\hat \beta$ and $\beta_0$, it follows that $\bar \beta \xrightarrow{p} \beta_0$ also. This means that

$$\frac 1n\sum_i\left [[y_i-h_i(\bar \beta)]\frac {\partial^2 h_i (\bar \beta)}{\partial \beta^2}\right] \xrightarrow{p} \frac 1n\sum_i\left [[y_i-h_i(\beta_0)]\frac {\partial^2 h_i (\beta_0)}{\partial \beta^2}\right] \\= \frac 1n\sum_i\left [E(u_i) E\frac {\partial^2 h_i (\beta_0)}{\partial \beta^2}\right] =0$$

because $E(u_i) =0$.

So this term vanishes asymptotically and we are left with (cancelling also the negative sings)

$$\sqrt n (\hat \beta -\beta) \xrightarrow{d} \left (\text {plim}\frac 1n\sum_i\left [\frac {\partial h_i (\beta_0)}{\partial \beta}\frac {\partial h_i (\beta_0)}{\partial \beta}\right] \right)^{-1} \cdot \left(\frac 1{\sqrt n} \sum_iu_i \frac {\partial h_i (\beta_0)}{\partial \beta} \right) $$

You have to assume that the first sum converges to something positive definite, and the second converges in distribution to a normal random variable, and you do indeed make these assumptions (or deeper ones that lead to them). I don't have the specific book you mention, but you can compare the above to Davidson &amp; McKinnon "Econometric Theory and Methods"(2004) ch. 6, around eq. $(6.30)$.