Rigorous derivation of relativistic energy-momentum relation I wish to derive the relativistic energy-momentum relation $E^2 = p^2c^2 + m^2 c^4$ following rigorous mathematical steps and without resorting to relativistic mass.
In one spatial dimension, given $p := m \gamma(u) u$ with $\gamma(u) := (1 - \frac{|u|^2}{c^2})^{-1/2}$, the energy would be given by
$$E = \int{ \frac{d}{dt}p \space dx}$$
I'm having a hard time with this this integration.

How is the relation $E^2 = p^2c^2 + m^2 c^4$ rigorously derived starting from relativistic momentum, without resorting to relativistic mass?


To give an idea of the rigour I expect in an answer, in example, an answer I'd accept for the derivation of $ E = \frac{1}{2} m v^2$ in classical mechanics would have been as follows:
We seek to integrate the differential form $F \space dx$. Parametrising $x$ by $t$, we obtain $dx = \frac{d}{dt} x \space dt$. 
The integral of interest is $\int F \space dx = m \int \frac{d^2}{dt^2}x \space dx = m \int (\frac{d^2}{dt^2}x) (\frac{d}{dt} x) dt$ after changing variables.
We recognize the integrand as $\frac{d}{dt} \left( \frac{1}{2} \left(\frac{d}{dt}x \right)^2 \right) $, and so the result $E = \frac{1}{2} m v^2$ follows from the fundamental theorem of calculus.

Again, as an example, a derivation of $E = \frac{1}{2} m v^2$ I would definitely not accept would be as follows:
$ \int F \space dx = m \int a \space dx = m \int \frac{dv}{dt} \space dx = m \int dv \frac{dx}{dt} $ = $ m \int v \space dv = \frac{1}{2}m v^2$.
Please carry out rigorous mathematical manipulations only.
 A: Since $P = Fv$ we have
$$\frac{dE}{dt} = \frac{dp}{dt} v$$
by Newton's second law. Integrating both sides with respect to $t$ gives
$$\int \frac{dE}{dt} \, dt = \int v \frac{dp}{dt} \, dt = \int v \, dp$$
by the chain rule, aka ordinary $u$-substitution. We have
$$p = \gamma m v = \frac{m v}{\sqrt{1-v^2}} \quad \Rightarrow \quad dp = \frac{m \, dv}{(1-v^2)^{3/2}}$$
where I set $c = 1$ for convenience and used the quotient rule. Integrating with initial and final velocities zero and $v_0$ gives
$$E(v_0) - E(0) = \int_0^{v_0} \frac{mv}{(1-v^2)^{3/2}} \, dv = \frac{m}{\sqrt{1 - v_0^2}} - m.$$
At this point we cannot proceed further since we don't know the constant of integration. One can show by physical arguments that $E(0) = m$. Thus
$$E(v) = \frac{m}{\sqrt{1-v^2}}$$
as desired. This isn't a hard derivation, but you're right: a lot of textbooks botch it.
A: For completeness, here's an arguably cleaner and simpler formulation of @knzhou 's answer:
We obtain
$$E = \int_{0}^{x_0} (\frac{d}{dt} p) \space dx = \int_{0}^{t_0} (\frac{d}{dt}  p) \space v \space dt = \int_{0}^{p_0} v \space dp = \int_{0}^{v_0} v \space (\frac{d}{dv} p) \space dv$$
by applying a sequence of reparametrizations $dx = v \space dt$, $dp = (\frac{d}{dt} p) \space dt$ and $dp = (\frac{d}{dv} p) \space dv$ to the integral for $E$. Since $ \frac{d}{dv} p = m \space (1 - \frac{v^2}{c^2})^{-3/2}$, it follows that
$$ E = \int_{0}^{v} \dfrac{m v}{(1-\frac{v^2}{c^2})^{3/2}} dv  = \frac{mc^2}{(1 - \frac{v^2}{c^2})^{1/2}} - mc^2.$$ 
Defining the total energy $\Sigma = E + mc^2$, since $\Sigma = \gamma m c^2$ and $p = \gamma m v$, it is easy to see by direct computation that $\Sigma^2 - c^2 p^2 = m^2 c^4$, hence
$$\Sigma^2 = m^2 c^4 + c^2 p^2 \space .$$
A: I want to elaborate a little bit inspired by the comment made by AccidentalFourierTransform and how it relates to the answer by knzouh. You start from the four-velocity
$$
u = \gamma~( c, \mathbf{v}),
$$
on which you base your definition of the four-momentum as
$$
p = m_0\gamma~( c, \mathbf{v}) = ( m_0\gamma c,m_0\gamma  \mathbf{v}) =  ( m_0\gamma c,m_0\gamma  \mathbf{v}) .
$$
Something you know about the four-velocity is, that is always squares to $u^\mu u_\mu=\gamma^2 (c^2-v^2) = c^2$ from which you immediately conclude that $p^\mu p_\mu = m_0 c^2$, thus representing a Lorentz scalar. On the other hand, if you recognize $ m_0\gamma c^2$ as the energy $E$ and the term $m_0\gamma  \mathbf{v}$ as the momentm $\mathbf{p}$ this can as well be written as $p^\mu p_\mu = E^2 /c^2 - \mathbf{p}^2$ which thus leads to desired formula
$$
E^2 /c^2 - \mathbf{p}^2 = m_0c^2 .
$$
In order for this argument to work we need to justify why we interpret $p^0$ as the energy. One possible way to proceed is the relativistic action $S$ of a free particle and define $E=- \partial_t S$ and $\vec{p} = \nabla_\mathbf{x} S$ as outlined for example here.
