First I want to define with the Stirling numbers of the first kind $\left[ \begin{array}{c} n \\ k \end{array} \right]$ a special generalization of the Riemann Zeta function :
$$\zeta_n(m):=\sum\limits_{k=1}^\infty \frac{1}{k^m}\left(\frac{n!}{(k-1)!}\left[\begin{array}{c} k \\ n+1 \end{array} \right]\right)$$
and
$$\eta_n(m):=\sum\limits_{k=1}^\infty \frac{(-1)^{k-1}}{k^m}\left(\frac{n!}{(k-1)!}\left[\begin{array}{c} k \\ n+1 \end{array} \right]\right)$$
which are convergent for the integer values $\enspace m\geq 2$ .
For $\enspace n=0\enspace$ we have $\enspace\zeta_0(m)=\zeta(m)\enspace$ and $\enspace\eta_0(m)=\eta(m)\enspace$ .
Note: Obviously (because of the other results) these series can be expressed by sums of the polylogarithm function and modifications of that.
Please also see here, part Expansion by harmonic numbers, with $\enspace\displaystyle w(n,m):=\frac{m!}{(n-1)!}\left[ \begin{array}{c} n \\ {m+1} \end{array} \right]\enspace$ and it's recursion formula.
Secondly, an extension of an integral as a series, $n\in\mathbb{N}_0$ and $z\in\mathbb{R}\setminus \{2\mathbb{N}\}$ and $nz>-1$:
$ \displaystyle \int\limits_0^\pi x^n \left(2\sin\frac{x}{2}\right)^z dx=i^{-z} \int\limits_0^\pi x^n e^{i\frac{xz}{2}}(1- e^{-ix})^z dx= e^{-i\frac{\pi z}{2}} \int\limits_0^\pi x^n \sum\limits_{k=0}^\infty\binom{z}{k}(-1)^k e^{-ix(\frac{z}{2}-k)} dx$
$\displaystyle =\int\limits_0^\pi x^n e^{i(x-\pi)\frac{z}{2}} dx+ \sum\limits_{v=0}^n \frac{(-1)^v\pi^{n-v} n!}{i^{v+1}(n-v)!} \sum\limits_{k=1}^\infty \binom{z}{k}\frac{1}{(\frac{z}{2}-k)^{v+1}} $
$\displaystyle \hspace{3.5cm} -i^{n-1}n!e^{-i\frac{\pi z}{2}} \sum\limits_{k=1}^\infty \binom{z}{k}\frac{ (-1)^k}{(\frac{z}{2}-k)^{n+1}}$
using the main branch of the logarithm and therefore $\displaystyle i=e^{i\frac{\pi}{2}}$ .
The Stirling numbers of the first kind are usually defined by $\enspace \displaystyle \sum\limits_{k=0}^n \left[ \begin{array}{c} n \\ k \end{array} \right] x^k := x(x+1)…(x+n-1) $ .
Because of $\enspace \displaystyle (\sum\limits_{v=0}^\infty x^v \frac{d^k}{dz^k}\binom{z}{v}) |_{z=0} =\frac{d^k}{dz^k}(1+x)^z |_{z=0} =(\ln(1+x))^k=k!\sum\limits_{v=k}^\infty (-1)^{v-k} \left[\begin{array}{c} v \\ k \end{array} \right] \frac{x^v}{v!}$
we get $\enspace \displaystyle \binom{z}{k}|_{z=0}=0^k\enspace$ , $\enspace \displaystyle \frac{d}{dz} \binom{z}{k} |_{z=0} = (-1)^{k-1} \left[\begin{array}{c} k \\ 1 \end{array} \right] \frac{1}{k!}= \frac{(-1)^{k-1}}{k} \enspace$ , $\enspace \displaystyle \frac{d^2}{dz^2} \binom{z}{k} |_{z=0} = (-1)^{k-2} \left[\begin{array}{c} k \\ 2 \end{array} \right] \frac{2!}{k!}= \frac{(-1)^k 2}{k}\sum\limits_{j=1}^{k-1}\frac{1}{j} \enspace$ and $\enspace \displaystyle \frac{d^3}{dz^3} \binom{z}{k} |_{z=0} = (-1)^{k-3} \left[\begin{array}{c} k \\ 3 \end{array} \right] \frac{3!}{k!}= \frac{(-1)^{k-1} 3}{k}( (\sum\limits_{j=1}^{k-1}\frac{1}{j})^2 - \sum\limits_{j=1}^{k-1}\frac{1}{j^2} ) $ .
For $(n;k):=(3;3)$ follows
$\displaystyle \int\limits_0^\pi x^3 \left(\ln\left(2\sin\frac{x}{2} \right)\right)^3 dx =$
$\hspace{2cm}\displaystyle =\frac{9\pi^2}{2}\left(\zeta(5)+3\eta(5)-4\eta_1(4)+2\eta_2(3)\right) $
$\hspace{2.5cm}\displaystyle - 90\left(\zeta(7)+\eta(7)\right) +72\left(\zeta_1(6)+\eta_1(6)\right) - 18\left(\zeta_2(5)+\eta_2(5)\right) $
Note:
For the calculations I have used $\enspace\displaystyle\int\limits_0^\pi x^n e^{iax}dx = \frac{(-1)^{n+1} n!}{(ia)^{n+1}}+e^{i\pi a}\sum\limits_{v=0}^n\frac{(-1)^v \pi^{n-v}n!}{(ia)^{v+1}(n-v)!}$
with $\enspace\displaystyle a=-(\frac{z}{2}-k)$ .
And it was necessary to calculate $\enspace\displaystyle\frac{d^m}{dz^m} \binom{z}{k}\frac{1}{(\frac{z}{2}-k)^{v+1}}|_{z=0}\enspace$ and $\enspace\displaystyle\frac{d^m}{dz^m} e^{-i\frac{\pi z}{2}}\binom{z}{k}\frac{1}{(\frac{z}{2}-k)^{n+1}}|_{z=0}\enspace$ for $\enspace m\in\{0,1,2,3\}$ .
Since $\int \exp(b x^2+c x) dx= \frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{2 b x+c}{2 \sqrt{b}}\right)}{2 \sqrt{b}}$ we integrate by parts and we have:
\begin{eqnarray}
I&=& -\frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{c}{2 \sqrt{b}}\right)}{2 \sqrt{b}}- \frac{2 a}{\sqrt{\pi}} \int\limits_0^\infty \frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{2 b x+c}{2 \sqrt{b}}\right)}{2 \sqrt{b}} \exp(-a^2 x^2) dx\\
&=& -\frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \text{erfi}\left(\frac{c}{2 \sqrt{b}}\right)}{2 \sqrt{b}}+\frac{e^{-\frac{c^2}{4 b}}\sqrt{\pi}}{\imath \sqrt{b}}2 T(\epsilon, \imath \frac{\sqrt{b}}{a}, \imath \frac{c \sqrt{2}}{2 \sqrt{b}} ) \\
&=& -\frac{\sqrt{\pi } e^{-\frac{c^2}{4 b}} \left(4 i T\left(\frac{i c}{\sqrt{2} \sqrt{b} \sqrt{1-\frac{b}{a^2}}},\frac{i \sqrt{b}}{a}\right)-\text{erfi}\left(\frac{c}{2 \sqrt{b}
\sqrt{1-\frac{b}{a^2}}}\right)+\text{erfi}\left(\frac{c}{2 \sqrt{b}}\right)\right)}{2 \sqrt{b}}
\end{eqnarray}
where in the second line we took a small number $0 < \epsilon << 1$ and we used the definition of the generalized Owen's T function Generalized Owen's T function and in the last line we simplified the result. Here $T(\cdot,\cdot)$ is the Owen's T function. The result is valid for $0 < b < a^2$.
In[697]:= {a, c} = RandomReal[{0, 3}, 2, WorkingPrecision -> 50];
b = RandomReal[{0, a^2}, WorkingPrecision -> 50]; eps = 10^(-9);
NIntegrate[Erfc[a x] Exp[b x^2 + c x], {x, 0, Infinity}]
-((E^(-(c^2/(4 b))) Sqrt[\[Pi]] Erfi[c/(2 Sqrt[b])])/(2 Sqrt[b])) +
2/Sqrt[Pi] a NIntegrate[(
E^(-(c^2/(4 b))) Sqrt[\[Pi]] Erfi[(c + 2 b x)/(2 Sqrt[b])])/(
2 Sqrt[b]) Exp[-a^2 x^2], {x, 0, Infinity}]
-((E^(-(c^2/(4 b))) Sqrt[\[Pi]] Erfi[c/(2 Sqrt[b])])/(2 Sqrt[b])) + (
E^(-(c^2/(4 b))) Sqrt[ Pi])/(I Sqrt[b])
NIntegrate[
Erf[I (c + 2 b /(Sqrt[2] a) x)/(2 Sqrt[b])] Exp[-1/2 x^2]/Sqrt[
2 Pi], {x, 0, Infinity}]
-((E^(-(c^2/(4 b))) Sqrt[\[Pi]] Erfi[c/(2 Sqrt[b])])/(2 Sqrt[b])) + (
E^(-(c^2/(4 b))) Sqrt[ Pi])/(I Sqrt[b])
2 T[eps, I b /(a Sqrt[b]), I (c Sqrt[2])/(2 Sqrt[b])]
-((E^(-(c^2/(4 b)))
Sqrt[\[Pi]] (Erfi[c/(2 Sqrt[b])] -
Erfi[c/(2 Sqrt[b] Sqrt[1 - b/a^2])] +
4 I OwenT[(I c)/(Sqrt[2] Sqrt[b] Sqrt[1 - b/a^2]), (I Sqrt[b])/
a]))/(2 Sqrt[b]))
Out[699]= 0.789518
Out[700]= 0.789518
Out[701]= 0.789518 + 0. I
Out[702]= 0.789518482586510679235860756093903252836337770103 +
0.*10^-59 I
Out[703]= 0.789518482636559413924687222227564264392988055015
Best Answer
Is there a specific reason you need an asymptotic expansion?
If you only care about showing that the integral tends to $dD/dt$ in the limit, then you just need to formulate the integral as an approximate identity. Let $H(s)$ denote the unit step function, i.e. $$ H(s) = \begin{cases} 1 & s\geq 0\\ 0 & s < 0 \end{cases} $$ and define $$ f(s) = \exp(-s)H(s). $$ Note that for all $\theta>0$, $H(s/\theta) = H(s)$. Setting $g(s) = \frac{dD}{ds}(s)$, your integral is equal to $$ \frac{1}{\theta}\int_{-\infty}^\infty f\left(\frac{t-s}{\theta}\right)g(s)~ds. $$ Writing $f_\theta(s) = (1/\theta)f(s/\theta)$, we can therefore rewrite the integral as the convolution $(f_\theta * g)(t)$. It then remains to show that $f_\theta$ is an approximate identity (or mollifier, or good kernel, whatever your preferred language is). This is simple as $f$ is an $L^1$ function, $\int_\mathbb{R} f(s)~ds = 1$, and $f_\theta$ are just the 1-dimensional rescalings of $f$. Therefore $(f_\theta * g)(t) \to g(t)$ pointwise almost everywhere as $\theta\to\infty$ for any reasonably nice function $g$. The convergence can be upgraded to stronger than pointwise a.e. depending on how nice $g$ is.