Integrals involving fractions of exponentials I am trying to calculate the average degree of a complex network, which requires me to solve for the following integral:
$$\int \mathrm{d} x \frac{\exp{\left[-x -\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right]}}{1 + a\exp{(-x)}},$$
where $\mu$, $\sigma$, and $a$ are known constants. Any advice on how to tackle this integral would be greatly appreciated.
So far, I have tried integration by parts, as well as substitutions of the type $s = \exp{(-x)}$ or identities in terms of the hyperbolic trig functions. Unfortunately, none of these got me far.
 A: Risch algorithm "is a method for deciding whether a function has an elementary function as an indefinite integral, and if it does, for determining that indefinite integral."
Mathematica knows this algorithm.
However, it returns the integral unevaluated, even in the special case when $\mu=0,\sigma=1,a=1$, and even when using the Rubi package:

So, it is highly unlikely that this indefinite integral is an elementary function.
A: We can obtain power series in $a$ by expanding in geometric series. Since the specific treatment depends on the integration region, let's introduce lower and upper limits, $y_0 <y$,
$$
I = \int_{y_0 }^{y} dx\, \frac{\exp \left[ -x -\frac{1}{2} \left( \frac{x-\mu }{\sigma } \right)^{2} \right] }{1+a\exp (-x)}
$$
For $|a|<e^{y_0 } $, we have $|a\exp (-x)| < 1$ and we can expand in a geometric series as
\begin{eqnarray*}
I &=& \sum_{n=0}^{\infty } (-a)^n \int_{y_0 }^{y} dx\, \exp \left[ -(n+1) x -\frac{1}{2} \left( \frac{x-\mu }{\sigma } \right)^{2} \right] \\
&=& \sum_{n=0}^{\infty } (-a)^n \sqrt{\frac{\pi }{2} } \sigma \ e^{((n+1) \sigma^{2} -2\mu)(n+1)/2} \cdot \\
& & \hspace{2.5cm} \left[ \mbox{erf} \left( \frac{(n+1)\sigma^{2} -\mu +y}{\sqrt{2} \sigma }\right) - \mbox{erf} \left( \frac{(n+1)\sigma^{2} -\mu +y_0 }{\sqrt{2} \sigma }\right) \right]
\end{eqnarray*}
On the other hand, for $|a|>e^y $, we have $|\exp (x) /a| <1$, and we can alternatively write and expand
\begin{eqnarray*}
I &=& \frac{1}{a} \int_{y_0 }^{y} dx\, \frac{\exp \left[-\frac{1}{2} \left( \frac{x-\mu }{\sigma } \right)^{2} \right] }{1+\exp (x)/a} \\
&=& \frac{1}{a} \sum_{n=0}^{\infty } \left(-\frac{1}{a} \right)^{n} \int_{y_0 }^{y} dx\, \exp \left[ nx -\frac{1}{2} \left( \frac{x-\mu }{\sigma } \right)^{2} \right] \\
&=& \sum_{n=0}^{\infty } \left(-\frac{1}{a} \right)^{n+1} \sqrt{\frac{\pi }{2} } \sigma \ e^{(n \sigma^{2} +2\mu)n/2} \cdot \\
& & \hspace{2.5cm} \left[ \mbox{erf} \left( \frac{n\sigma^{2} +\mu -y}{\sqrt{2} \sigma }\right) - \mbox{erf} \left( \frac{n\sigma^{2} +\mu -y_0 }{\sqrt{2} \sigma }\right) \right]
\end{eqnarray*}
Finally, for $|a|\in \ ]e^{y_0 } ,e^y [$, we can separate into two integration intervals, $[y_0,\ln |a|\, [$ and $]\, \ln |a|,y]$, and use the applicable expansion from above in each case to obtain the combined result. In the case of negative $a$, one will have to specify how the boundary at $\ln |a|$ is approached from either side, corresponding to the pole in the original integrand; for example, one might specify a principal value prescription.
A: Decided to mention in more detail what I said in a comment to the question.
The corresponding definite integral is expressible through the Mordell integral.
From the introduction of Vector-valued higher depth quantum modular forms and higher Mordell integrals by Bringmann, Kaszian and Milas (free version in arXiv):

The Mordell integral is usually defined as a function of two variables$$h(z)=h(z;\tau):=\int_{\mathbb R}\frac{\cosh(2\pi zw)}{\cosh(\pi w)}e^{\pi i\tau w^2}dw,\tag{1.1}$$where $z\in\mathbb C$ and $\tau\in\mathbb H$, the complex upper half-plane. Integrals of this form were studied by many mathematicians including Kronecker, Lerch, Ramanujan, Riemann, Siegel, and of course Mordell, who proved that a whole family of integrals reduces to $(1.1)$.

Two questions here on MO about Mordell integrals: Asymptotic expansion of the Mordell integral, Conceptual meaning of a non-linear relation connecting $6$ Mordell integrals? It is also mentioned in my own (incomplete) answer to one of my questions here.
A: Note
$$
\int  \frac{\exp{\left[-x -\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right]}}{1 + a\exp{(-x)}}dx
=
\int \frac{\exp{\left[-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right]}}{ \exp(x)+ a}dx .
$$
Even this
$$
\int\frac{e^{-x^2/2}}{e^x+1}dx
$$
is not known to Maple, although
$$
\int_{-\infty}^{+\infty}\frac{e^{-x^2/2}}{e^x+1}dx = \frac{\sqrt{\pi}}{\sqrt{2}} .
$$
