Your (related?) question brought me here. I'm not sure if I'm following you right, but I guess your first question asks: is there {$h_1,...,h_8$},a basis of the Cartan subalgebra, such that $${\rm tr}_{L(\Lambda_0)} q^{L(0)-1/3} x_1^{h_1} \cdots x_8^{h_8}$$ equals the right hand side for some $\alpha$? (Here $L(\Lambda_0)$ is your basic representation). I assume you are flexible on the choice of a basis of the Cartan, right? For instance, your formula obviously doesn't hold if you take (for some purposes natural) $h_i=\omega_i$ (fundamental co-weights).
Let me ignore the denominator which doesn't change.
Start from the usual Euclidean $\mathbb{Z}^8$-lattice with orthonormal basis $\epsilon_i$. Then, $D_8$ lattice, as you know, consists of all $v \in \mathbb{Z}^8$ such that the sum of coordinates of $v$ is even. To get the $E_8$-lattice you also include all $v \in (\mathbb{Z}+\frac{1}{2})^8$ such that (again) the sum of coordinates of $v$ is even.
Set
$$\theta_3(q^{1/2},x)=\sum_{n \in \mathbb{Z}} q^{n^2/2} x^{n}$$
$$\theta_4(q^{1/2},x)=\sum_{n \in \mathbb{Z}}(-1)^n q^{n^2/2} x^n $$
$$\theta_2(q^{1/2},x)= \sum_{n \in \mathbb{Z}} q^{(n+1/2)^2/2} x^{n+1/2}$$
$$\theta_1(q^{1/2},x)=\sum_{n \in \mathbb{Z}} (-1)^n q^{(n+1/2)^2/2} x^{n+1/2} $$
So you will be summing over the $D_8$-lattice and another "shifted" $D_8$.
For the $D_8$ part, you get
$$\prod_{i=1}^8 \theta_3(q^{1/2},x_i)+\prod_{i=1}^8 \theta_4(q^{1/2},x_i)$$
multiplied by $\frac{1}{2}$.The shifted part contributes with
$$\prod_{i=1}^8 \theta_2(q^{1/2},x_i)+\prod_{i=1}^8 \theta_1(q^{1/2},x_i)$$
multiplied by $\frac{1}{2}$. The resulting formula looks very similar to yours (with $h_i=\epsilon_i$).
This construction is essentially decomposition
$L(\Lambda_0)=L_{D_8}(\Lambda_0) \oplus L_{D_8}(\Lambda_8)$ into
irreducible (standard) $D^{(1)}_8$-modules (I hope $\Lambda_8$ is correct here). All this is very natural if you are familiar with the spinor construction of ${E_8}^{(1)}$ (see this monograph ).
Just in case you prefer $${\rm tr}_{L(\Lambda_0)} q^{L(0)-1/3} x_1^{\omega_1} \cdots x_8^{\omega_8},$$ simply express $\omega_i$ as a linear combination of $\epsilon_j$'s, and rearrange $x_j$'s accordingly.
I hope this helps. My apologies if I said something wrong.
n.b. $x_i$ should be $e^{2 \pi i x_i}$ to be consistent with " $\theta_1$ vanishes at $x_i=0$".
EDIT: (Too long for a comment.) Regarding your second question and the fact that up to degree $8$ in $x_i$'s there is only contribution from the power sum $p_2(x)=x_1^2+ \cdots + x_8^2$.
Actually something similar happens if you consider one-point functions on the torus so it might be related. For a holomorphic rational CFT $V$ (resp. VOA) of central charge $c$, meaning that there is a unique irreducible sector (resp. module), the vector space of one-point functions (essentially traces of degree zero operators coming from vertex operators)
$$\langle Z(v,\tau), v \in V \rangle$$
is a fairly nice space of modular forms with a mild pole at infinity. Basically after you multiply with $\eta^{c}(\tau)$ you expect to get all holomorphic modular form.
For the ($c=24$) moonshine module this was shown by Dong and Mason here , but I think it must also hold for $c=8k$ (shown maybe in this more recent paper.) .
Now back to $V=L_{E_8}(\Lambda_0)$, which is holomorphic with $c=8$. If you accept what I said, each $\eta^8(\tau) Z(v,\tau)$ is holomorphic (possibly zero) for $\Gamma(1)$ of weight $deg(v)$. So if you pick $v$ of degree $2,4,6$ you will be generating one-dimensional spaces spanned by $\frac{E_6(\tau)}{\eta^8}$, $\frac{E_8(\tau)}{\eta^8}$ and $\frac{E_{10}(\tau)}{\eta^8}$. However in degree $8$ you have a contribution from the cusp form $\Delta$, so the dimension jumps by one. If I remember correctly, up to degree $8$ you can use conformal vector to get those, but higher you also have to construct
some primary vectors coming from spherical harmonics to get a nonzero $Z(v,\tau)$, and then you apply a result of Waldspurger . So I agree with you, the Weyl invariant of degree $8$ (for $E_8$) is most likely relevant in your case.
Dong and Mason (and their collaborators/students) wrote several papers on this subject, so it might be useful to read some of their work. For instance, this paper , might be also useful for your purposes.
Best Answer
An eta-product identity for $k=3$, $s=1$ similar to that given in my comment above may be found in "Some eta-identities arising from theta series" by Günter Köhler (Math. Scand. 66, 1990, p. 146, identity (3)): $$ \frac{\eta^2(\tau)\eta^2(4\tau)}{\eta(2\tau)}=\sum_{n=1}^\infty\left(\frac n3\right)ne^\frac{2\pi in^2\tau}3, $$ which is equivalent to $$ \prod_{n=1}^\infty\frac{(1-q^n)^2(1-q^{4n})^2}{1-q^{2n}}=\sum_{n=1}^\infty\left(\frac n3\right)nq^{\frac{n^2-1}3}=1-2q+4q^5-5q^8+7q^{16}-8q^{21}+... $$ The rhs is the same as $q^{-\frac13}\sum_{n\in\mathbb Z}(3n+1)q^{\frac13{(3n+1)^2}}$, i. e., up to rescalings, $\frac\partial{\partial z}$ of your $\theta_1(z;\tau)_3$ at $z=0$.
I don't know whether such identities for other $k$ are systematically treated anywhere. Köhler also has a book "Eta Products and Theta Series Identities" but I have never looked inside. From "About this book":
(Second part) This might be a separate answer but I decided to just add it here. In "Eta-quotients and theta functions" by Robert J. Lemke Oliver (free online version available), a classification of all those theta-functions of the form $\theta_\psi(z)=\sum_n\psi(n)n^\delta q^{n^2}$ which admit an eta-product decomposition is given. Here $\psi$ is a Dirichlet character, and $\delta=0$ or $1$ according to the parity of $\psi$. There turns out to be very few of this kind - just eight in the even case and five in the odd case. The case of your $k=3$, $s=1$ is the third identity of Theorem 1.1 (2) in that paper; I am not reproducing it here as it is essentially the same as the Köhler's case above. If one allows some twists of the characters, there are a little bit more such products.