I would like to try and evaluate the following gamma function inspired sum.
$$\Xi(z)=\sum_{t=1}^{\infty}\frac{t^z}{e^t}$$
According to my computations, for large $z$,
$$\Xi(z)\approx\Gamma (z+1)$$
and perhaps even
$$\Xi(z) \sim \Gamma (z+1)$$
Does a closed form exist for this sum?
Your sum is in the form of the Polylogarithm. In fact, it is equal to $\operatorname{Li}_{-z}(1/e).$ When $z$ is a negative integer, this sum is easily computed in closed form using the identity $\displaystyle \frac{d}{dx} \operatorname{Li}_n (x) = x\operatorname{Li}_{n-1} (x).$
By applying the Abel-Plana summation to the Polylogarithm series, we get
$$\operatorname{Li}_s(z) = {z\over2} + {\Gamma(1 \!-\! s, -\ln z) \over (-\ln z)^{1-s}} + 2z \int_0^\infty \frac{\sin(s\arctan t \,- \,t\ln z)} {(1+t^2)^{s/2} \,(e^{2\pi t}-1)} \,\mathrm{d}t $$
and so $$\operatorname{Li}_{-z}(1/e) = \frac{1}{2e} + \Gamma(z+1,1) + \frac{2}{e} \int^{\infty}_0 \frac{ (1+t^2)^{z/2} \sin(-z \tan^{-1}t+t)}{e^{2\pi t} -1} dt .$$
where $\Gamma(s,x)$ is the incomplete gamma function. Your asymptotic would be explained if you could show why the remaining integral is comparatively small.
It appears that we can obtain a simple proof of the proposed asymptotics using harmonic sum techniques when $\Re(z) = \sigma > 0$.
To see this, introduce the sum $$S(x; z) = \sum_{n\ge 1} \frac{(nx)^z}{\exp(nx)}$$ so that $$\Xi(z) = S(1; z).$$
The sum term is harmonic and may be evaluated by inverting its Mellin transform.
Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have $$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad g(x) = \frac{x^z}{\exp(x)}.$$ We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty \frac{x^z}{\exp(x)} x^{s-1} dx = \int_0^\infty e^{-x} x^{s+z-1} dx.$$
No additional computation is necessary as this integral is the defining integral of the gamma function $$\Gamma(s) = \int_0^\infty e^{-x} x^{s-1} dx$$ so that $$g^*(s) = \Gamma(s+z)$$ with fundamental strip $\langle -\sigma, +\infty \rangle.$
It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x; z)$ is given by
$$Q(s) = \zeta(s) \Gamma(s+z) \quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(s).$$ The Mellin inversion integral here is $$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$ which we evaluate by shifting it to the left for an expansion about zero.
The line for the integral is chosen in the intersection of the fundamental strip of $g^*(s)$ and the half-plane of convergence of $\zeta(s)$, which is $\Re(s)>1.$
Collecting the residues from the poles we first obtain that $$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = \frac{1}{x} \Gamma(1+z).$$
The remaining poles are from the gamma function. Keeping in mind that $\Re(z) = \sigma > 0$ we find with $q\ge 0$ $$\mathrm{Res}\left(Q(s)/x^s; s=-z-q\right) = x^{z+q} \frac{(-1)^q}{q!} \zeta(-z-q).$$
We conclude that $$\Xi(z) = S(1; z) \sim \Gamma(1+z) + \sum_{q\ge 0} \frac{(-1)^q}{q!} \zeta(-z-q).$$
For $z=m$ with $m$ a positive integer this simplifies to $$\Xi(m) \sim m! - \sum_{q\ge 0} \frac{(-1)^q}{q!} \frac{B_{m+q+1}}{m+q+1}.$$
We may answer the original question in the affirmative, it is true that $$\Xi(z) \sim \Gamma(1+z).$$
This MSE link points to another interesting Mellin transform computation.
Addendum. The convergence properties of the original sum as proposed above are quite interesting, there is an initial segment where it appears to diverge before it then begins to converge. This observation just in case someone decides to study the numerics of this problem.
$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle #1 \right\rangle} \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace} \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left( #1 \right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ \begin{align} \color{#00f}{\large\Xi\pars{z}}& = \sum_{t = 1}^{\infty}{t^{z} \over \expo{t}} =\sum_{t = 1}^{\infty}{\pars{\expo{-1}}^{t} \over t^{-z}} = \color{#00f}{\large{\rm L_{i}}_{-z}\pars{1 \over e}} \end{align} where $\ds{{\rm L_{i}}_{s}\pars{z}}$ is the Polylogarithm Function.
According to the 'link' given above, we can get the following relations: $$ \lim_{\Re\pars{z} \to -\infty}\Xi\pars{z} = {1 \over \expo{}}\,,\qquad\qquad \Xi\pars{z} \sim \Gamma\pars{1 + z}\quad\mbox{when}\quad\Re\pars{z} \gg 1 $$
