Firstly,
$$I_1=\int\limits_0^e x^{-x}\,\text dx = e\int\limits_0^1 (ex)^{-ex}\,\text dx = e\int\limits_0^1 e^{-ex(1+\ln x)}\,\text dx
= e\int\limits_0^1\,\sum_{k=0}^\infty \dfrac{(-ex(1+\ln x))^k}{k!}\,\text dx$$
$$= \sum_{k=0}^\infty \dfrac{(-1)^k e^{k
+1}}{k!} \int\limits_0^1\, x^k(1+\ln x)^k\,\text dx
= e+\sum_{k=1}^\infty \dfrac{(-1)^k e^{k+1}}{k!} \sum_{j=0}^k\dbinom kj \int\limits_0^1\, x^k(\ln x)^j\,\text dx$$
$$= e+\sum_{k=1}^\infty \dfrac{(-1)^k e^{k+1}}{k!} \sum_{j=0}^k\dbinom kj \left(\dfrac2{k+1}\right)^{j+1} \int\limits_0^1\, \dfrac{k+1}2\,x^k \ln^j\left(x^{\large\frac{k+1}2}\right)\,\text dx$$
$$= e + \sum_{k=1}^\infty \dfrac{(-1)^k e^{k+1}}{k!} \sum_{j=0}^k\dbinom kj \left(\dfrac2{k+1}\right)^{j+1} \int\limits_0^1\, y\ln^j\left(y\right)\,\text dy$$
$$=e+\sum_{k=1}^\infty \dfrac{(-1)^k e^{k+1}}{k!} \sum_{j=0}^k\dbinom kj \left(\dfrac2{k+1}\right)^{j+1} (-1)^j \dfrac{j!}{2^{j+1}}$$
$$=e+\sum_{k=1}^\infty (-1)^k e^{k+1} \sum_{j=0}^k\dfrac{(-1)^j}{ (k-j)!(k+1)^{j+1}}$$
$$=e + \sum_{k=1}^\infty (-1)^k e^{k+1} e^{-k-1}\dfrac{(-1)^k}{(k+1)^{k+1} k!}\Gamma(k+1,-k-1),$$
$$I_1=e + \sum_{k=2}^\infty \dfrac{\Gamma(k,-k)}{k^k \Gamma(k)}$$
received an alternating series, wherein
$$e + \sum_{k=2}^{22}\dfrac{\Gamma(k,-k)}{k^k \Gamma(k)}
\le 1.96466\,14544\,475,\quad
\dfrac{\Gamma(23,-23)}{\Gamma(23)}{23^{-23}}\le 2\cdot 10^{-13}.$$
Therefore,
$$I_1\ge 1.96466\,14544\,477.$$
$\color{green}{\textbf{Edit of 21.02.23}}$
Obtained estimation can be significantly improved.
Let $x\in\big(a,(1+\delta a\big),\;\delta\ll1,\;$ then
$$x=a(1+\delta z)=a +\delta a z,\qquad \big(z\in(0,1),\;a\in(1,\infty)\big),$$
$$x^{-x}=a^{\large -a-\delta a z}\;\left(1+\delta z \right) ^{\large -a(1+\delta z)},$$
Taking in account the inequality
$$(1+t)^{\large-a(1+t)} -e^{-at} \le \dfrac a2 t^2\left(1-\dfrac{1+3a}3 t\right)\tag1$$
and identity
$$\int\limits_0^d t^k a^{-a t}\,\text dt =\dfrac{Γ(k+1)-Γ(k+1,ad\ln a)}{(a\ln a)^{k+1}},\tag2$$
one can get
$$I(a, \delta)=\int\limits_a^{a+\delta a} x^{-x}\,\text dx
= \delta a\cdot a^{-a} \int\limits_0^1 a^{-\delta az}\left(1+\delta z\right)^{\large- a(1+\delta z)}\,\text dz$$
$$= a^{2-a}\int\limits_0^\delta a^{-at}\left(1+t\right)^{\large-a(1+t)}\,\text dt
\le a^{2-a} \int\limits_0^\delta a^{\large-at}
\left(e^{-at}+\frac a2 t^2 - \frac {a(3a+1)}6 t^3\right)\,\text dt$$
$$=a^{1-a} \int\limits_0^\delta (ea)^{\large-a t}\,\text dt
+ \dfrac12a^{2-a} \int\limits_0^\delta t^2 a^{\large-a t}\,\text dt
-\dfrac {3a+1}6 a^{2-a} \int\limits_0^\delta t^3 a^{\large-a t}\,\text dt,$$
$$\color{brown}{\mathbf{I(a, \delta)=a^{-a}\left(\dfrac{1-(ea)^{-\large \delta a}}{\ln a}
+\dfrac{2-\Gamma(3,a\delta\ln a)}{2a\ln^3a}
-(3a+1)\dfrac{6-\Gamma(4,a\delta\ln a)}{6a\ln^4a}\right)}},\tag3$$
$$\int\limits_e^\infty x^{-x} dx = \lim\limits_{\delta \to 0} \sum\limits_{k=0}^\infty I(e\delta^k, \delta).\tag4$$
Comparing of this expression (for $\delta=\sqrt[\large 6000]2-1,\; n=3000$) with the numeric value of the integral is shown in the table below
$$\left[\begin{matrix}
[a,b] & I(a,\delta) & \text{numeric}\\
\left[e, e\sqrt2\right] & 0.02848\,08301\,58887 & 0.02848\,08301\,58764\\
\left[e\sqrt2, 2e\right] & 0.00227\,71750\,39440 & 0.00227\,71750\,39412\\
\left[2e, 2e\sqrt2\right] & 0.00003\,64476\,51537 & 0.00003\,64476\,51535\\
\left[2e\sqrt2,4e\right] & 0.00000\,00502\,01224 & 0.00000\,00502\,01224\\
[e,4e] & 0.03079\,45030\,51082 & 0.03079\,45030\,50937
\end{matrix}\right]$$
Besides,
$$(1+t)^{-a(1+t)}\le e^{-at},$$
$$I(a,\infty)=\int\limits_a^\infty x^{-x}\,\text dx
=a^{-a}\int_0^\infty a^{-at}\left(1+t\right)^{\large- a(1+t)}\,\text dt$$
$$\le a^{-a} \int\limits_0^\infty a^{-at} e^{-at}\text dz
= \dfrac{a^{-a}}{\ln(ae)} (ae)^{-at}\bigg|_0^\infty
= \dfrac{a^{-a}}{\ln(ae)},$$
$$I(4e,\infty) \le \dfrac{(4e)^{-4e}}{\ln(4e^2)} \le 1.592\cdot 10^{-12}.$$
Finally,
$$I\le 1.96466\,14544\,477 + 0.03079\,45030\,51082 + 1.5917\cdot 10^{-12}
< \color{brown}{\mathbf{1.99545\,59576}}.$$
Numeric calculations give
$$I\approx 1.99545\,59575.$$
$\textbf{Previous version.}$
On the other hand,
$$I_2=\int\limits_e^\infty x^{-x}\,\text dy
=e\int\limits_0^\infty (e(1+y))^{-e(1+y)}\,\text dy
=\int\limits_0^\infty \dfrac{e^{1-e(1+y)}}{(1+y)^3} (1+y)^{3-e(1+y)}\,\text dy=I_{21}+I_{22},$$
where
$$I_{21}=\int\limits_1^2 \dfrac{e^{1-ez}}{z^3} z^{3-ez}\,\text dz,\quad
I_{22}=\int\limits_2^\infty \dfrac{e^{1-ez}}{z^3}\, z^{3-ez}\,\text dz.$$
Applying inequality $z^{3-ez}\le f(z),$ where
$$f(z)=\begin{cases}
e^{-0.7177ez}(-26.937+49.92z-15.94z^2),\quad\text{if}\quad z\in[1,2)\\[8pt]
(z-0.48) e^{6-1.4868ez},\quad\text{if}\quad z\in[2,+\infty)
\end{cases}$$
(see also WA plot), one can get
$$I_{21}\leq \int\limits_1^2 \dfrac{e^{1-ez}}{z^3}\,(e^{0.7177ez}(-26.937+49.92z-15.94z^2))\,\text dz$$
$$=\int\limits_1^2 \dfrac1{z^3}\,(e^{1-1.7177ez}(-26.937+49.92z-15.94z^2))\,\text dz$$
$$=\left(\dfrac{36.6112 - 306.641z}{z^2}\,e^{-4.66919 z} - 1475.1\operatorname{ Ei}(-4.66919 z)\bigg|_1^2\right)\le 0.03105\,33417\,253,$$(see also WA checking)
$$I_{22}\leq \int\limits_2^\infty \dfrac{e^{1-ez}}{z^3}\,((z-0.48) e^{6-1.4868ez})\,\text dz
=\int\limits_2^\infty \dfrac{z-0.48}{z^3}\,e^{7-2.4868ez})\,\text dz$$
$$=\dfrac{263.192-2875.76 z}{z^2}e^{-6.75982 z} - 19439.7 Ei(-6.75982 z) \bigg|_2^\infty \le 0.00003\,70357\,268$$
(see also WA checking).
Finally,
$$I=I_1+I_{21}+I_{22}\le 1.96466\,14544\,477+0.03105\,33417\,253+0.00003\,70357\,268,$$
$$\color{brown}{\mathbf{I\le1.99575\,18318\,998 \le \pi-\ln\pi.}}$$
