statistical mechanics - Evaluating low-temperature dependence of the BCS gap function

How does one go about evaluating the behavior of the BCS gap $\Delta = \Delta(T)$ for $T \to 0^+$ under the weak coupling approximation $\Delta/\hbar\omega_D \ll 1$?

In Fetter & Walecka, Quantum theory of Many-Particle Systems, Prob. 13.9 it is said that the starting point is

$$\tag{1} \ln\frac{\Delta_0}{\Delta} = 2\int_0^{\hbar\omega_D}{\frac{\mathrm d\xi}{\sqrt{\xi^2+\Delta^2}}\frac{1}{e^{\beta\sqrt{\xi^2+\Delta^2}}+1}},$$ which I have no problem deriving from the theory, but I can't find a way to actually evaluate this integral even under the approximations $\hbar\omega_D \to \infty$, $\Delta \approx \Delta_0$ (in the RHS) and $\beta\Delta \to \infty$. [Of course $\beta = (k_BT)^{-1}$ and $\Delta_0 = \Delta(T = 0)$.] I have tried several approaches, used different Taylor-expansions and changes of variables, but I am simply stuck.

For the record, the expected behavior is supposed to be

$$\tag{2} \Delta(T) \sim \Delta_0\left(1 - \sqrt{\frac{2\pi}{\beta\Delta_0}}e^{-\beta\Delta_0}\right) .$$

EDIT: Just leaving it here for the posterity. I found a more complete way to tackle this integral; specifically, under the WC approximation one has

$$\int_0^{+\infty}{\frac{\mathrm d x}{\sqrt{x^2 + 1}}\frac{e^{-\beta\Delta\sqrt{x^2 + 1}}}{1 + e^{-\beta\Delta\sqrt{x^2 + 1}}}} = \int_1^{+\infty}{\frac{\mathrm d y}{\sqrt{y^2 - 1}}\frac{e^{-\beta\Delta y}}{1 + e^{-\beta\Delta y}}} =$$ $$= \sum_{k=1}^{+\infty}\int_1^{+\infty}{\frac{\mathrm d y}{\sqrt{y^2 - 1}} (-1)^{k+1}e^{-k\beta\Delta y}} = \sum_{k=1}^{+\infty}(-1)^{k+1}\int_0^{+\infty}{\mathrm d t\; e^{-k\beta\Delta \cosh{t}}} = \sum_{k=1}^{+\infty}(-1)^{k+1} K_0(k\beta\Delta),$$

$K_0$ being the 0-order modified Bessel function of the second kind, whose asymptotic behavior is known and may be used to solve the problem in a relatively clean way (and even find the corrections at higher orders, which are $\in O(e^{-\beta\Delta}(\beta\Delta)^{-k - 1/2})$. Cf. Abrikosov, Gorkov, Dzyaloshinski, Methods of Quantum Field Theory in Statistical Phyisics, 1963. Pagg. 303-304.

Answer

Hints:

1. 
   
   

   Define difference $\delta:=\Delta-\Delta_0$. Deduce from $|\delta|\ll |\Delta_0|$ that the lhs. of eq. (1) is

   $$\tag{A}\text{lhs}~\approx~ -\frac{\delta}{\Delta_0}.$$
   
   


2. 
   

   Substitute $\xi=x\Delta$ in the integral on the rhs. of eq. (1). Deduce using $\hbar \omega_D \gg \Delta$ that the rhs. is

   $$\tag{B} \text{rhs}~\approx~ \int_{\mathbb{R}} \! \frac{dx}{\sqrt{1+x^2}} \frac{1}{e^{\beta\Delta \sqrt{1+x^2}}+1}.$$
   
   


3. 
   

   Deduce from $\beta\Delta\gg 1$ that we can simplify the rhs. further to a Gaussian integral

   $$\tag{C} \text{rhs}~\approx~ \int_{\mathbb{R}} \! dx~ e^{-\beta\Delta (1+\frac{1}{2}x^2)}~=~\sqrt{\frac{2\pi}{\beta\Delta}}e^{-\beta\Delta} .$$ Such arguments are closely related to the method of steepest descent.
   
   


4. 
   

   Deduce eq. (2).