$$X_{n+1}=\frac{1}{2}(X_n)(1+\frac{a}{(X_n)^2})$$
From this scheme, we know by Newton-Raphson that this is for $x = \sqrt{a}$. Let's say $\theta$ be the solution.
$\theta^2 = a$
To prove convergence of order p, we need to show :
$|X_{n+1} - \theta| = \lambda * |X_n - \theta|^p$ where $\lambda$ is the asymptotic error constant.
I am getting :
$$|X_{n+1} - \theta| = \frac{1}{2*X_n} |X_n - \theta|^2$$
Can we say it is of order $2$ ? $\lambda$ here is not constant. It is dependent on $X_n$.
Some general considerations: if $f(x_0)=0$ and $f'(x_0) \neq 0$ then the standard error estimate for Newton's method looks like this:
$$x-\frac{f(x)}{f'(x)}-x_0=x-\frac{f'(x_0)(x-x_0)+1/2 f''(x_0)(x-x_0)^2+o((x-x_0)^2)}{f'(x_0)+O(x-x_0)}-x_0 \\ = x-(1+O(x-x_0)) \left (x-x_0+1/2 \frac{f''(x_0)}{f'(x_0)}(x-x_0)^2+o((x-x_0)^2) \right ) - x_0 \\ = -1/2 \frac{f''(x_0)}{f'(x_0)}(x-x_0)^2+o((x-x_0)^2).$$
Some consequences: you will not have exact equality as you wrote it. What you hope for is an inequality, $|X_{n+1}-x_0| \leq \lambda |X_n-x_0|^2$ for some $\lambda$. In general this does not hold globally, it only holds once $X_n$ is in some interval around $x_0$. To understand this interval better you need to look at what that $o((x-x_0)^2)$ term actually is. Assuming you have a third continuous derivative, we can do the calculation over again using explicit Lagrange remainders:
$$x-\frac{f(x)}{f'(x)}-x_0=x-\frac{f'(x_0)(x-x_0)+1/2f''(x_0)(x-x_0)^2+1/6 f'''(\xi_x)(x-x_0)^3}{f'(x_0)+f''(\eta_x)(x-x_0)}-x_0$$
where $\xi_x$ and $\eta_x$ are numbers between $x_0$ and $x$. Simplifying a bit you wind up at:
$$(x-x_0) \left ( 1-\frac{1+1/2 f''(x_0)(x-x_0)/f'(x_0)+1/6 f'''(\xi_x)(x-x_0)^2/f'(x_0)}{1+f''(\eta_x)(x-x_0)/f'(x_0)} \right ).$$
Simplifying further by getting a common denominator and you wind up at:
$$(x-x_0)^2 \left ( \frac{f''(\eta_x)/f'(x_0)-1/2 f''(x_0)/f'(x_0)-1/6 f'''(\xi_x)(x-x_0)/f'(x_0)}{1+f''(\eta_x)(x-x_0)/f'(x_0)} \right ).$$
So in general we shouldn't expect global quadratic convergence because of that third derivative term.
When you specialize the above to the case of Newton's method for square roots, you get something which is quite algebraically obvious. I could've derived this from the start, but I wanted to go over the general considerations as well.
$$\frac{x+a/x}{2}-\sqrt{a}=\frac{1}{2x} (x-\sqrt{a})^2.$$
Thus we indeed do not have global quadratic convergence for square roots; the coefficient blows up as $x \to 0$. (This should not be surprising, because if $X_1$ is the same when $X_0=x$ and when $X_0=a/x$, for any $x$.) Thus a really robust error estimate should try to argue that $X_n$ can't get too close to zero. A good result that I think holds is $X_n \geq \alpha := \min \{ X_0,\sqrt{a} \}$. In this case it you get your desired estimate with $\lambda=\frac{1}{2\alpha}$.
For $a>0$ and $0<X_1\ne \sqrt a,$ there do not exist constants $\lambda, p $ such that $|X_{n+1}-\sqrt a|=\lambda |X_n=\sqrt a|^p$ for all $n.$ We have $$X_{n+1}-\sqrt a=\frac {X_n}{2}( 1+\frac {a}{X_n^2}) -\sqrt a=$$ $$=\frac {1}{2X_n}(X_n^2-2X_n\sqrt a +a) =\frac {1}{2X_n}(X_n-\sqrt a)^2 .$$ So if $X_n\ne \sqrt a$ then $X_{n+1}\ne \sqrt a. $ But $X_n$ converges to $\sqrt a$ so for large $n$ the ratio of $|X_{n+1}-\sqrt a|$ to $|\frac {1}{2\sqrt a} (X_n-\sqrt a)^2|$ converges to $1.$
Historical note: This special case of Newton's method (as it is called in older books) was described by Heron of Alexandria about 2000 years ago.