Here is another simple (symmetrization) trick, which may be better known, though. Let $Y$ be an independent copy of $X$. Note that $\ln x$ is increasing in $x>0$. So, twice the difference between the left- and right-hand sides of the inequality $(*)$ in question is $$(**)\qquad E(Y-X)(\ln Y-\ln X)X^a Y^a\ge0,$$ since $(Y-X)(\ln Y-\ln X)X^a Y^a\ge0$. Inequality $(**)$, and hence $(*)$, are strict iff the random variable $X$ is non-degenerate; that is, iff $P(X=b)<1$ for all real $b$. We also see that the assumptions $X>e$ and $0<a<1$ can be relaxed to $X>0$ and $a\in\mathbb{R}$.
Addendum: I would like to put the above solution into the following perspective.
The Chebyshev integral association inequality states that $$(1)\qquad\int_I fg\,d\mu\ge \int_I f\,d\mu\,\int_I g\,d\mu,$$ where $\mu$ is a probability measure on an interval $I\subseteq\mathbb R$ and functions $f$ and $g$ are both increasing (or both decreasing) on $I$. A formally more general version of this inequality is $$(2)\qquad\int_I fg\,d\nu\,\int_I d\nu\ge \int_I f\,d\nu\,\int_I g\,d\nu $$ for any (say) finite measure $\nu$ on $I$. However, $(2)$ follows immediately from $(1)$ by taking $\mu=\nu/\int_I d\nu$.
One can view (the nonstrict version of) the inequality $(*)$ in question as an instance of $(2)$, with $I=(0,\infty)$, $f(x)\equiv x$, $g(x)\equiv\ln x$, and $\nu(dx)=x^a\,P(X\in dx)$. On the other hand, of course $(2)$ is proved by the same symmetrization argument as the one used directly in the above solution.
Since you said you were also interested in "similar inequalities but with more than one variable", you may want to look at the FKG inequality, which extends the Chebyshev integral association inequality to increasing functions on a lattice in place of the one-dimensional interval $I$, and possibly at further extensions and variations of FKG; see e.g. [FKG-Wikipedia]. However, then you need to impose certain structural restrictions on the measure $\nu$.