Rigorously prove the period of small oscillations by directly integrating This answer proved that
$$\lim_{E\to E_0}2\int_{x_1}^{x_2}\frac{\mathrm dx}{\sqrt{2\left(E-U\!\left(x\right)\right)}}=\frac{2\pi}{\sqrt{U''\!\left(x_0\right)}},$$
where $E_0:=U\!\left(x_0\right)$ is a minimal value of the second-differentiable potential $U$.
However, I don't think this proof is rigorous enough in that it just throws away terms with order higher than $\left(x-x_0\right)^2$ in $U$.
There are some skeptical points.
We cannot just expand the integrand in Taylor series.
We cannot guarantee that the integrand is analytic at $x_0$ (we cannot even assume that $U$ has third or higher derivatives).
Also, we haven't argued that those terms do not have significant effects on the integral boundaries (i.e. the value of $x_1$ and $x_2$).
I hope there will be an $\varepsilon$-$\delta$ version of the proof.
 A: While I won't try to make a rigorous proof in this answer, I'll point out some important remarks that I believe are sufficient to address this question.

*

*As mentioned in the comments, "small oscillations" means expanding the potential to second-order. We assume the oscillations are small enough to neglect higher order terms.


*This might be the most important point. The function doesn't need to be analytic to be expanded to second order. You don't need a Taylor series. A Taylor polynomial with a convenient form of the remainder (see Wikipedia on Taylor's Theorem) is enough. This requires way weaker assumptions.


*Due to the small oscillations hypotheses, higher order terms can't make a difference on the boundary. "Small oscillations" means $x_1$ and $x_2$ are sufficiently close to $x_0$ such that those terms are negligible.


*I don't believe this expression holds exactly, but using Taylor's theorem you'll probably be able to derive an estimate on the error. Again, this doesn't require analyticity of the involved quantities, only being differentiable/continuously differentiable a couple of times.
A: If $U$ is twice-differentiable at the point $x_0$ and $U''(x_0)\neq 0$, then there exists a function $R$ such that
$$U(x) = U(x_0) + U'(x_0) (x-x_0) +\frac{1}{2} U''(x_0)(x-x_0)^2 + R(x)(x-x_0)^2$$
and $\lim_{x\rightarrow x_0} R(x) = 0$.  Given that $x_0$ is taken to be a stable equilibrium point of $U$, we have that $U'(x_0)= 0$ and $U''(x_0)>0$. We therefore have that
$$\sqrt{E-U(x)} = \sqrt{E-E_0 - \frac{1}{2} U''(x_0)(x-x_0)^2 - R(x)(x-x_0)^2}$$
where in accordance with OP's notation we define $U(x_0) \equiv E_0$.  For $E>E_0$, we may define $\epsilon \equiv E-E_0$ and write
$$\sqrt{E-U(x)} = \sqrt{\epsilon} \sqrt{1- \frac{1}{2\epsilon}\big[U''(x_0)+2R(x)\big] (x-x_0)^2}$$
From this it follows that we may write
$$\frac{\mathrm dx}{\sqrt{E-U(x)}}= \frac{dx/\sqrt{\epsilon}}{\sqrt{1-\frac{1}{2}\big[U''(x_0)+2R(x)\big]\left(\frac{x-x_0}{\sqrt{\epsilon}}\right)^2}}$$
$$= \frac{d\xi}{\sqrt{1-\frac{1}{2}\big[U''(x_0)+2R\left(\xi\sqrt{\epsilon}+x_0\right)\big] \xi^2}}$$
where we've defined $\xi \equiv (x-x_0)/\sqrt{\epsilon}$. The turning points between which we are integrating occur when $U(x)=E$, which can be rewritten implicitly as
$$\xi= \pm\frac{1}{\sqrt{\frac{1}{2}\big[U''(x_0)+2R\big(\xi\sqrt{\epsilon}+x_0\big)\big] }}$$
Let the solutions to this equation be $\xi_1$ and $\xi_2$, respectively. Our integral then becomes
$$T_{\epsilon} = \frac{2}{\sqrt{2}}\int_{\xi_1}^{\xi_2} \frac{\mathrm d\xi}{\sqrt{1-\frac{1}{2}\big[U''(x_0)+2R\left(\xi\sqrt{\epsilon}+x_0\right)\big] \xi^2}}$$
$\epsilon$ appears in the argument of $R$ as well as in the integral bounds, so taking the $\epsilon\rightarrow 0$ limit is a bit subtle - however, I claim that we obtain the correct answer by naively taking all three limits simultaneously (see below for a more rigorous justification). If we do so, we may note that $\lim_{\epsilon\rightarrow 0} R\big(\xi\sqrt{\epsilon}+x_0\big) = 0$ and $\lim_{\epsilon\rightarrow 0} \xi_{1/2} = \pm \xi_0$ where $\xi_0 \equiv \sqrt{2/U''(x_0)}$. Plugging all of this in, we obtain
$$\lim_{\epsilon\rightarrow 0} T_\epsilon =\frac{2}{\sqrt{2}}\int_{-\xi_0}^{\xi_0} \frac{1}{\sqrt{1-\frac{1}{2}U''(x_0) \xi^2}}$$
We conclude in the obvious way - define $\sigma \equiv \xi \sqrt{U''(x_0)/2}$ to obtain
$$\lim_{\epsilon\rightarrow 0}T_\epsilon = \frac{2}{\sqrt{U''(x_0)}}\int_{-1}^1 \frac{d\sigma}{\sqrt{1-\sigma^2}} = \frac{2\pi}{\sqrt{U''(x_0)}}$$

For those interested, I'll justify my naive limit-taking more rigorously. Let $\{\epsilon_n\}$ and $\{\alpha_n\}$ be two positive sequences tending to zero such that $0<\epsilon<\epsilon_n \implies |R\big(\xi\sqrt{\epsilon} + x_0\big)|<\alpha_n$ for all each $\xi\in[\xi_1,\xi_2]$.
Define $\gamma_n = \sqrt{\frac{2}{U''(x_0)+2\alpha_n}}$ and observe that $\xi_1 < -\gamma_n$ and $\gamma_n < \xi_2$. We may now define the sequence
$$t_n = \frac{2}{\sqrt{2}}\int_{-\gamma_n}^{\gamma_n} \frac{\mathrm d\xi}{\sqrt{1-\frac{1}{2}\big[U''(x_0)+2R\big(\xi\sqrt{\epsilon_n}+x_0\big)\big]\xi^2}}$$
and note that $\lim_{n\rightarrow \infty} t_n = \lim_{\epsilon\rightarrow 0}T_\epsilon$.
Note that the integrand is bounded below by $\left(1-\frac{1}{2}\big[U''(x_0)-2\alpha_n\big]\right)^{-1/2}$.  The integral of this lower bound is
$$L_n = \frac{2}{\sqrt{2}}\int_{-\gamma_n}^{\gamma_n} \frac{\mathrm d\xi}{\sqrt{1-\frac{1}{2}\big[U''(x_0)-2\alpha_n\big]\xi^2}}$$
$$= \frac{4}{\sqrt{U''(x_0)-2\alpha_n}}\sin^{-1}\left(\sqrt{\frac{U''(x_0)-2\alpha_n}{U''(x_0)+2\alpha_n}}\right)$$
Similarly, the integrand is bounded above by  $\left(1-\frac{1}{2}\big[U''(x_0)+2\alpha_n\big]\right)^{-1/2}$, and the integral by
$$U_n = \frac{2}{\sqrt{2}}\int_{-\gamma_n}^{\gamma_n} \frac{\mathrm d\xi}{\sqrt{1-\frac{1}{2}\big[U''(x_0)+2\alpha_n\big]\xi^2}} = \frac{2\pi}{\sqrt{U''(x_0)+2\alpha_n}}$$

Because $L_n < t_n < U_n$ and $\lim_{n\rightarrow \infty} L_n = \lim_{n\rightarrow \infty} U_n = \frac{2\pi}{\sqrt{U''(x_0)}}$, we have by the squeeze theorem that $\lim_{n\rightarrow \infty} t_n = \lim_{\epsilon\rightarrow 0} T_\epsilon = \frac{2\pi}{\sqrt{U''(x_0)}}$.
