I want to evalute the leading order term of the following integral as a series of $1/N$ and $\epsilon$,
$\int_{-\epsilon}^\epsilon dt (\sqrt{1-(\rho+t)^2}-\sqrt{1-\rho^2})e^{-N t^2}$, where $\rho=1-\epsilon$, $\epsilon\rightarrow 0$, $N\rightarrow \infty$.
I think the difficulties is related to the singularity of derivation of function $f(\rho)=\sqrt{1-\rho^2}$ at $\epsilon\rightarrow0$. Namely, we cannot just expand $f(\rho)$ using Taylor series to the second order. All higher order terms will contribute significantly to the result. It also seems unclear how to apply Waston's lemma since the integral is from $-\epsilon$ to $\epsilon$. Is there a way to solve this problem?
While it is true that the integrand is singular at one of the endpoints of integration, this is not the biggest problem about the question you're asking. In fact, it is the requirements for a series in $\epsilon$ and $1/N$ simultaneously to exist that seem irreconcilable. First of all, note that the two obvious pieces the integral can be broken into exist separately, and we can evaluate the second one explicitly:
$$\sqrt{1-\rho^2}\int_{-\epsilon}^\epsilon dt e^{-Nt^2}=\sqrt{\epsilon(2-\epsilon)}\sqrt{\frac{\pi}{N}}~\text{erf}(\epsilon\sqrt{N})$$
There is no way to expand this integral in the series form that you desire using traditional expansions (convergent or asymptotic). No matter what, one cannot obtain series with ascending powers of $\epsilon$ and descending powers $N$. The first part of the integral is a bit more complex, but it can be rewritten in the form:
$$\int_{-\epsilon}^\epsilon dt \sqrt{1-(\rho+t)^2}e^{-N t^2}=\int_0^{2\epsilon}\sqrt{t(2-t)}e^{-N(t-\epsilon)^2}$$
Now expand the term $\sqrt{2-t}=\sqrt{2}\sum_{n}c_n (t/2)^n$, for some combinatorial coefficients $c_n:=\Gamma(3/2)/(\Gamma(3/2-n)\Gamma(n+1))$ and note that
$$\int_{0}^{2\epsilon}t^{n+1/2}e^{-N(t-\epsilon)^2}dt=\epsilon^{n+3/2}\int_{-1}^1(1-x)^{n+1/2}e^{-N\epsilon^2x^2}dx:=\epsilon^{n+3/2}f_n(\epsilon\sqrt{N})$$
and therefore
$$\int_{-\epsilon}^\epsilon dt \left(\sqrt{1-(\rho+t)^2}-\sqrt{1-\rho^2}\right)e^{-N t^2}=\sum_n \frac{c_n}{2^{n-1/2}}\epsilon^{n+3/2}\left(f_n(\epsilon\sqrt{N})-\sqrt{\pi}\frac{\text{erf}(\epsilon\sqrt{N})}{\epsilon\sqrt{N}}\right)$$
and now we have run into the exact same situation as above. The appearance of functions of the variable $\epsilon\sqrt{N}$ preclude us from obtaining the desired series. Another way to understand the issue at hand is to look at the first equation in three cases: $\epsilon^{-1}=N, \epsilon^{-2}=N, \epsilon^{-4}=N$. Clearly, all these cases satisfy the conditions of the question, however they all result in wildly different first order large $N$ asymptotics for the integral in the first equation. Thus a unique asymptotic expansion in two variables for the integral in question cannot possibly exist.