Is it "wellknown" which odd powers of the theta function are eigenforms for the halfintegral weight Hecke operators? If so, what is a good reference? Is there a slick algorithm for proving computationally that a specific weight k/2 noncusp form is an eigenform?

The case of even powers was (asked and) solved by P. T. Bateman in "Problem E 2051", Amer. Math. Monthly 76 (1969), the solution being that the powers 2, 4, 8 are the only even powers of $\Theta$ that are eigenforms. A similar idea shows that indeed the only odd powers of theta that are eigenforms are 1, 3, 5 and 7:
Proof. The idea is to compute explicit expressions for the first two fourier coefficients of $\Theta^k$ and $T_{3^2} \Theta^k$, depending only on $k$. For a fixed $n$, in order to compute the $n$th coefficient of $\Theta^k$ we consider "squarepartitions" of $n$  all the ways to write $n$ as a sum of squares of positive numbers, regardless of order. Each such partition contributes to all large enough powers of $\Theta$, specifically, as soon as the power is greater than the number of summands in the partition. We can calculate exactly how much a squarepartition contributes to a coefficient. Namely, the partition: $$n = \alpha_1 \cdot n_1^2 + \alpha_2 \cdot n_2^2 + \ldots + \alpha_r \cdot n_r^2$$ where the $n_i$ are distinct, constributes to the $n$th coefficient of $\Theta^k$: $$2^{\alpha_1+\ldots+\alpha_r} \binom{k}{\alpha_1} \binom{k\alpha_1}{\alpha_2}\cdot \ldots\cdot \binom{k\alpha_1\ldots\alpha_{r1}}{\alpha_r}$$ The power of 2 accounts for different choices of sign, and each binomial coefficient accounts for the number of ways left to choose the positions of $n_i$'s. If $n$ is fixed, each such expression is a polynomial in $k$. So, we can sum over all squarepartitions of $n$ to get the $n$th coefficient of $\Theta^k$. Now we need to get the coefficients of $T_{3^2}\Theta^k$. This is well known. We use theorem 1.7 of "On Modular Forms of Half Integral Weight", G. Shimura (1973), which can also be found as proposition 13 of chapter 4 in "Introduction to Elliptic Curves and Modular Forms", N. Koblitz (1984) (I am slightly weakening the theorem to exclude $p=2$):
For $\Theta^k$ we have $N = 4$ and $\chi = 1$. Using the above with $p=3$: $$b_1 = a_9  (3)^{\lambda1}a_1$$ $$b_2 = a_{18} + (3)^{\lambda1}a_2$$ And $a_1$, $a_2$, $a_9$, $a_{18}$ are polynomials in $k$, or $\lambda$ (which is more convenient). If $\Theta^k$ is an eigenform for $T_{3^2}$, we must have $\frac{b1}{a1}=\frac{b2}{a2}$. We will show this is doesn't happen for $k\gt 7$. The idea is simple: writing $b_1a_2b_2a_1$ out we see a polynomial of $\lambda$ of degree 19, and $(3)^{\lambda1}$ times a quadratic polynomial in $\lambda$. As $\lambda \rightarrow \infty$ (so $k\rightarrow \infty$ as well), the exponential factor dominates, and the expression cannot be zero. We compute a lower bound for this, and check remaining cases manually. The last part isn't very enlightening, but since you expressed keenness for an algorithm (and someone might want to check that I'm not rambling nonsense), I will include the sage computation below.


