The following is simply a comment, grown large.
Since $x^{1/x} = \exp\left(\frac{\log(x)}{x} \right)$, Bell polynomial can be used to expresses the $n$-th order derivative using Faà di Bruno's formula:
$$
\frac{\mathrm{d}^n}{\mathrm{d}x^n} f(g(x)) = \sum_{k=1}^n f^{(k)}(g(x)) B_{n,k}\left(g^\prime(x),g^{\prime\prime}(x),\ldots, g^{(n)}(x)\right)
$$
Taking $f = \exp$ and $g(x) = \frac{\log(x)}{x}$ and using $f^\prime = f$, as well as a known result:
$$
\frac{\mathrm{d}^n}{\mathrm{d}x^n} \frac{\log(x)}{x} = (-1)^n\frac{n!}{x^{n+1}} \left(\log(x) - H_n\right)
$$
we have:
$$
\frac{\mathrm{d}^n}{\mathrm{d}x^n} x^{1/x} = x^{1/x} \sum_{k=1}^n B_{n,k}\left(-\frac{\log(x)-1}{x},\ldots, (-1)^n\frac{n!}{x^{n+1}} \left(\log(x) - H_n\right) \right)
$$
Using well known homogeneity properties of the Bell polynomial:
$$\begin{eqnarray}
B_{n,k}\left( \lambda y_1, \lambda^2 y_2, \ldots, \lambda^n y_n\right) &=& \lambda^n B_{n,k}\left( y_1, y_2, \ldots, y_n\right) \\
B_{n,k}\left( \lambda y_1, \lambda y_2, \ldots, \lambda y_n\right) &=& \lambda^k B_{n,k}\left( y_1, y_2, \ldots, y_n\right)
\end{eqnarray}
$$
the expression can be further simplified:
$$
\frac{\mathrm{d}^n}{\mathrm{d}x^n} x^{1/x} = (-1)^n x^{1/x-n} \sum_{k=1}^n \frac{1}{x^k} B_{n,k}\left(\log(x)-1,\ldots, n! \left(\log(x) - H_n\right) \right)
$$
Bell polynomials are multivariate polynomials with positive coefficients. Notice that this implies that all the real roots are bounded from above by $\exp(H_n)$.
Here are few low order Bell polynomials:
$$
\begin{array}{c|ccccc}
n\backslash k & 1 & 2 & 3 & 4 & 5 \\ \hline
1 & y_1 & \text{} & \text{} & \text{} & \text{} \\
2 & y_2 & y_1^2 & \text{} & \text{} & \text{} \\
3 & y_3 & 3 y_1 y_2 & y_1^3 & \text{} & \text{} \\
4 & y_4 & 3 y_2^2+4 y_1 y_3 & 6 y_1^2 y_2 & y_1^4 & \text{} \\
5 & y_5 & 10 y_2 y_3+5 y_1 y_4 & 10 y_3 y_1^2+15 y_2^2 y_1 & 10 y_1^3 y_2 & y_1^5
\end{array}
$$
Notice that this implies that derivatives of $x^{1/x}$ evaluated at $x=1$ are all integers. In fact this is known as sequence A008405.
Also, here are few zeros computed explicitly in Mathematica:
$$
\begin{array}{c|llllll}
n & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ \hline
1 & \mathrm{e} & \text{} & \text{} & \text{} & \text{} & \text{} \\
2 & 0.581933 & 4.36777 & \text{} & \text{} & \text{} & \text{} \\
3 & 0.342086 & 0.826389 & 6.00406 & \text{} & \text{} & \text{} \\
4 & 0.246206 & 0.453114 & 1.05675 & 7.63532 & \text{} & \text{} \\
5 & 0.19403 & 0.312965 & 0.555416 & 1.28152 & 9.26413 & \text{} \\
6 & 0.16101 & 0.239864 & 0.373541 & 0.654014 & 1.5034 & 10.8915 \\
\end{array}
$$
The following was the command used:
zeros = N[
Table[Exp[
t] /. {ToRules[
Reduce[0 ==
BellY[Table[{1, n! Exp[-t] (t - HarmonicNumber[n])}, {n, 1,
m}]], t, Reals]]}, {m, 1, 6}]];
Notice that the largest root is quite close to the upper bound $\exp(H_n)$, here are few values, side-by-side, the upper row being actual larges zeros for $1 \leqslant n \leqslant 12$ and the lower row being approximate values of $\exp(H_n)$:
