nt.number theory - Odd powers of the theta function as eigenforms

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:

Proposition. If $k$ is an odd integer greater than $7$, then $Theta^k$ is not an eigenform for $T_{3^2}$.

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 "square-partitions" 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 square-partition 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 ldotscdot binom{k-alpha_1-ldots-alpha_{r-1}}{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 square-partitions 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$):

Theorem (Shimura). Let $p$ be an odd prime number, and $f in G_k(N, chi)$ (an integral modular form of half-integer weight $k/2$, level $N$, and character $chi$).
Put: $$f=sum_{n=0}^{infty} a_n q^n$$
$$T_{p^2}f=sum_{n=0}^{infty} b_n q^n$$
Then:
$$b_n = a_{p^2n} + chi (p) bigg(frac{(-1)^lambda n}{p}bigg) p^{lambda -1}a_n + chi(p^2) p^{k-2}a_{n/p^2}$$
where $lambda = frac{k-1}{2}$, and $a_{n/p^2}=0$ if $p^2 not| n$.

For $Theta^k$ we have $N = 4$ and $chi = 1$. Using the above with $p=3$:
$$b_1 = a_9 - (-3)^{lambda-1}a_1$$
$$b_2 = a_{18} + (-3)^{lambda-1}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 $kgt 7$.

The idea is simple: writing $b_1a_2-b_2a_1$ out we see a polynomial of $lambda$ of degree 19, and $(-3)^{lambda-1}$ times a quadratic polynomial in $lambda$. As $lambda rightarrow infty$ (so $krightarrow 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.

---

def square_partitions(n, min=1):
    """ Returns a list of partitions of n into squares, the least of which can be min^2. """
    if n == 0:
        return [[(-1,1)]]
    else:
        l = []
        for i in [min..sqrt(n)]:
            for p in square_partitions(n-i^2, i):
                if p[0][0] == i:
                    p[0] = (i, p[0][1]+1)
                else:
                    p = [(i,1)] + p
                l += [p]
        return l

def theta_coeff(k, n):
    """ Computes the nth coefficient of $theta^k$. """
    f = 0
    for p in square_partitions(n):
        g = 1
        j = k
        for a in p:
            if a[0] != -1:
                g *= 2^a[1] * binomial(j, a[1])
                j -= a[1]
        if (not k in ZZ) or j >= 0:
            f += g
    return f

R.<l,s> = QQ['l','s']
# s = (-3)^(l-1)

k = 2*l+1
a1, a2, a9, a18 = theta_coeff(k,1), theta_coeff(k,2), theta_coeff(k,9), theta_coeff(k,18)
b1, b2 = a9 - s*a1, a18 + s*a2

g = b1*a2-b2*a1
lc = LCM([denominator(c) for c in g.coefficients()])
g *= lc

# Check the quadratic factor increases in absolute value when $xge7$
x = QQ['x'].0
print factor(g.coefficient(s)(l=x)) # (-6252318072000) * x * (x + 1/2)^2

# Lower bound where exponential factor dominates. Can do better since we're not using the quadratic factor.
m1 = sum(abs(c) for c in g(s=0).coefficients())
print m1*116^19 < 3^115 # True

# Check remaining cases manually.
print [2*i+1 for i in [0..120] if g(i, (-3)^(i-1)) == 0] # [1, 3, 5, 7]