Why don't people use Hamilton's equations for a relativistic free particle? A relativistic free particle has the Hamiltonian in general:
$$ \mathcal{H} = \sqrt{{\bf p}^2c^2+m^2c^4}.$$
I read somewhere that says, it is possible to go further and say that the EoM are Hamilton's equations. But it is not done as there is "less interest" in such a discussion.
Is there something deeper to this? Like another formalism is ''better''.
(My guess is a more trivial one though. That is, it is not useful because the equations just get very cluttered and ugly)
 A: I) Here we will assume that OP is talking about a relativistic point particle with zero spin in a $d$-dimensional Minkowski spacetime with metric $\eta_{\mu\nu}$ of sign convention $(−,+,\ldots,+)$. Also we put $c=1$ for simplicity.
Note that the relativistic point particle has world-line reparametrization invariance, which is a gauge symmetry/redundancy in the formulation. We are (to a large extent) free to parametrize the world-line of the point particle in any way we wish. Let us call the world-line parameter for $\tau$ (which does not have to be the proper time). This gauge freedom can be encoded in an einbein field $e=e(\tau)>0$. The resulting Hamiltonian Lagrangian is$^1$
$$  L_H~:=~ p_{\mu} \dot{x}^{\mu} - \underbrace{\frac{e}{2}(p^2+m^2)}_{\text{Hamiltonian}}, \tag{1} $$
cf. e.g. this Phys.SE post.
Here dot means differentiation wrt. $\tau$. The square of the momentum vector is
$$\begin{align} p^2~:=~& \eta^{\mu\nu} p_{\mu} p_{\nu}\cr ~=~&-(p^0)^2+{\bf p}^2\cr
~=~&-2p^+p^- + {\bf p}_{\perp}^2, \end{align}\tag{2}$$
where we have used light-cone coordinates in the last expression.
II) Static gauge $x^0=\tau$.  If we integrate out $p^0$ and $e$, we get OP's square root model
$$\begin{align} \left. L_H\right|_{x^0=\tau}
\quad\stackrel{p^0}{\longrightarrow}&\quad
{\bf p}\cdot \dot{\bf x}- \underbrace{\left(\frac{1}{2e} + \frac{e}{2}({\bf p}^2+m^2)\right)}_{\text{Hamiltonian}}\cr\cr
\quad\stackrel{e}{\longrightarrow}&\quad
{\bf p}\cdot \dot{\bf x}  - \underbrace{\sqrt{{\bf p}^2+m^2}}_{\text{Hamiltonian}} .\end{align}\tag{3} $$
For sufficiently short$^2$ times $\Delta \tau=\tau_f-\tau_i$, the path integral becomes$^3$
$$\begin{align}& \langle {\bf x}_f,\tau_f \mid {\bf x}_i,\tau_i\rangle\cr 
~=~&i\hbar\Delta\tau\int_{\mathbb{R_+}} \!\frac{\mathrm{d}e}{2} \int_{\mathbb{R}^d} \!\frac{\mathrm{d}^dp}{(2\pi\hbar)^d}  \exp\left[\frac{i}{\hbar}\left( p_{\mu} \Delta x^{\mu} 
-\underbrace{\frac{e}{2}(p^2+m^2)}_{\text{Hamiltonian}}\Delta\tau\right)\right]\cr
~=~& \int_{\mathbb{R}^{d-1}} \!\frac{\mathrm{d}^{d-1}{\bf p}}{(2\pi\hbar)^{d-1}}
i\hbar\Delta\tau\int_{\mathbb{R_+}} \!\frac{\mathrm{d}e}{2}
~\underbrace{\frac{1}{\sqrt{2\pi\hbar ie\Delta\tau}}}_{\text{Gauss. } p^0\text{-int.}} \cr
&\exp\left[\frac{i}{\hbar}\left( {\bf p}\cdot \Delta {\bf x} -\underbrace{\left( \frac{1}{2e} + \frac{e}{2}({\bf p}^2+m^2)\right)}_{\text{Hamiltonian}}\Delta\tau\right) \right]\cr
~\stackrel{(6)}{=}~& \int_{\mathbb{R}^{d-1}} \!\frac{\mathrm{d}^{d-1}{\bf p}}{(2\pi\hbar)^{d-1}} \frac{\hbar}{2\sqrt{{\bf p}^2+m^2}} \exp\left[\frac{i}{\hbar}\left( {\bf p}\cdot \Delta {\bf x} - \Delta \tau \underbrace{\sqrt{{\bf p}^2+m^2}}_{\text{Hamiltonian}}\right)\right]\cr
~=~&i\hbar\Delta\tau\int_{\mathbb{R_+}} \!\frac{\mathrm{d}e}{2}
~\underbrace{\frac{1}{(2\pi\hbar ie\Delta\tau)^{d/2}}}_{\text{Gauss. } p\text{-int.}}
\exp\left[\frac{i}{2\hbar}\left(  \frac{(\Delta x)^2}{e\Delta\tau} - m^2e\Delta\tau\right) \right]\cr
~\stackrel{(6)}{=}~&\frac{1}{(2\pi)^{d/2}}\Big(\frac{m/\hbar}{ \sqrt{(\Delta x)^2}}\Big)^{\frac{d}{2}-1}K_{\frac{d}{2}-1}\Big(\frac{m}{\hbar}\sqrt{(\Delta x)^2}\Big)
,\end{align} \tag{4} $$
which also happens to be the standard scalar propagator $\langle\Omega|T[\phi (x_f)\phi (x_i)]|\Omega\rangle $ in QFT/2nd quantization, cf. e.g. Refs. 1-3. From a 2nd quantized perspective, the $e$-integration in eq. (4) is a Schwinger parametrization of the Fourier transformed propagator $$\frac{i}{\hbar}\langle\Omega|T[\widetilde{\phi} (p_f)\widetilde{\phi} (p_i)]|\Omega\rangle~=~\frac{\hbar^2}{p_f^2+m^2-i\epsilon}(2\pi\hbar)^d\delta^d(p_f\!+\!p_i). \tag{5} $$
As is well-known, eq. (4) is Lorentz covariant and falls off exponentially outside the light-cone. In eq. (4) we have used the integrals
$$\begin{align}
\int_{\mathbb{R}_+} \!\frac{\mathrm{d}e}{e^{1+\nu}}\exp\left[-ae-\frac{b}{e}\right]  ~=~&2\left(\frac{a}{b}\right)^{\nu/2} K_{\nu}\left(2\sqrt{ab}\right),\cr  
\int_{\mathbb{R}_+} \!\frac{\mathrm{d}e}{e^{1-\nu}}\exp\left[-ae-\frac{b}{e}\right]  ~=~&2\left(\frac{b}{a}\right)^{\nu/2} K_{\nu}\left(2\sqrt{ab}\right),\cr
\int_{\mathbb{R}_+} \!\frac{\mathrm{d}e}{\sqrt{e}}\exp\left[-ae-\frac{b}{e}\right]  ~=~&\sqrt{\frac{\pi}{a}} \exp\left[-2\sqrt{ab}\right],\cr {\rm Re}(a), {\rm Re}(b)~>~&0.\end{align}\tag{6} $$
III) Light-cone gauge $x^+=\tau$. If we integrate out $p^-$ and $e$, we get
$$\begin{align} \left. L_H\right|_{x^+=\tau}
\quad\stackrel{p^-,~e}{\longrightarrow}\quad & -p^+\cdot \dot{x}^-
+{\bf p}_{\perp}\cdot \dot{\bf x}_{\perp}\cr &  - \underbrace{\frac{{\bf p}_{\perp}^2+m^2}{2p^+}}_{\text{Hamiltonian}} .\end{align}\tag{7} $$
IV) We stress that the Euler-Lagrange (EL) equations for either of the Hamiltonian Lagrangians (1), (3), and (7) lead to Hamilton's equations. The point is now that physical quantities should not depend on the choice of gauge-fixing. We are free to use the most convenient gauge choice. Each formulation (1), (3), and (7) are valid, and have their pros and cons. The static gauge choice (3) is disfavored because of the square root.
References:

*

*M.E. Peskin & D.V. Schroeder, An Intro to QFT; eq. (2.50).


*M.D. Schwartz, QFT and the Standard Model; eq. (6.25).


*O. Corradini & C. Schubert, Spinning Particles in QM & QFT, arXiv:1512.08694; subsection 1.5.1, eqs. (1.160-162)


*T. Padmanabhan, QFT: The Why, What and How, 2016; subsections 1.3.1 + 1.4.4.
--
$^1$ Strictly speaking, there are also Faddeev-Popov ghost terms and gauge-fixing terms, which we have ignored for simplicity. These action terms are consistently generated in the BFV formulation, cf. e.g. my Phys.SE post here. The normalization factor in eq. (4) can be derived via Gaussian integration in the BFV formulation over the 2 bosonic variables $x^0$, $B$; and the 4 fermionic variables $\bar{C}$, $P$, $C$, $\bar{P}$.
$^2$ Here we just consider a single time slice for simplicity. The full path integral is the continuum limit of multiple time slice discretizations with insertion of corresponding completeness relations. It turns out that the result (4) for the free theory does not depend on the number of time slice discretizations.
$^3$ Here we use the Feynman $i\epsilon$-prescription ${\rm Re}(i\Delta\tau)>0$. The Gaussian integration over $p^0_E=i p^0_M$ becomes damped after a Wick-rotation $\tau_E=i\tau_M$, $x^0_E=ix^0_M$ to Euclidean signature.
A: Well, not really. We COULD write hamiltonian as square root - if we know, what is a square root of an operator. Of course we have simple approximation:
$$\sqrt{1+x}=1+\frac x2-\frac{x^2}{8}+O(x^3)$$
Using this we could write your hamiltonian as:
$$\mathcal H=mc^2\sqrt{1+\frac{p^2}{m^2c^2}}=mc^2+\frac{p^2}{2m}+O(p^4).$$
The problem is that this form of hamiltonian allows us for superluminal particle transport - evolving particle with this hamiltonian gives nonzero probability in large distances - very problematic for hamiltonian grown from relativistic theory.
The answer is to find other solution, which after getting square gives us the hamiltonian before. That's how you derive the Dirac equation - you assume some matrices in equation and you hope that you can get a solution:
$$\mathcal H_D=\alpha^ip_ic+\beta mc^2$$
$$\mathcal H_D^2=p^2c^2+m^2c^4$$
The following procedures are basics of every relativistic quantum theory script, for example: http://inspirehep.net/record/459479/files/Forshaw.pdf , page 9.
