I would like to simplify the following sum if possible: $$ e^{-x} \sum_{k=1}^\infty \frac{x^k}{k!} \sum_{j=1}^k \frac{\lambda^j}{j} $$ where $x \geq 0$ and $\lambda \in[0, 1]$.
When $\lambda = 1$ the inner sum is just the $k$th harmonic number, $H_k$, so the sum is known from the exponential generating function for the harmonic series: $$ e^{-x}\sum_{k=1}^\infty \frac{x^k}{k!} H_k = \operatorname{Ein}(x) = \Gamma(0,x) + \ln(x) + \gamma, $$ where $\gamma$ is Euler's constant and $\Gamma(\cdot)$ is the incomplete gamma function. And the case $\lambda = 0$ is, of course, easy.
But what about when $0 < \lambda < 1$? Any help appreciated.
(The sum arises in a problem involving the order statistics of a Poisson-distributed number of iid exponentially distributed random variables.)
Sounds like you want some kind of integral representation. Maybe try the following. Call the sum $f(\lambda)$. Then
$$f'(\lambda)=e^{-x}\sum_{k=1}^\infty\frac{x^k}{k!}\sum_{j=0}^{k-1}\lambda^j=\frac{e^{-x}}{(1-\lambda)}\sum_{k=1}^\infty \frac{x^k(1-\lambda^k)}{k!}=\frac{e^{-x}}{1-\lambda}\left(e^x-1-(e^{\lambda x}-1)\right)$$
(verify you can actually differentiate through the sum). So,
$$f(\lambda)=e^{-x}\int_0^\lambda\frac{e^x-e^{\lambda x}}{1-\lambda}d\lambda.$$
Which looks like an error function integral to me if you split up the numerator.