Nonnegativity of coefficients of a modular form defined in terms of the Jacobi thetanull functions Question
Let
\begin{align*}
\theta_2(q) & = \sum_{n=-\infty}^{\infty} q^{(n+1/2)^2}
\\
\theta_3(q) & = \sum_{n=-\infty}^{\infty} q^{n^2}
\\
\theta_4(q) & = \sum_{n=-\infty}^{\infty} (-1)^n q^{n^2}
\end{align*}
be the Jacobi thetanull functions. Define
$$
f(q) = \theta_3^{12} \theta_4^8  + \theta_3^8\theta_4^{12} + \theta_3^{12} \theta_2^8  + \theta_3^8\theta_2^{12}.
$$
This function has the $q$-series expansion
$$
f(q) = 2 + 240 q^2 + 10240 q^3 + 134640 q^4 + 1007616 q^5 + 5215680 q^6 + 
 20828160 q^7 + 69131760 q^8 + \ldots.
$$
Conjecture: The coefficients in the $q$-series expansion of $f(q)$ are all nonnegative.
Can someone help me prove this?
Motivation
Maryna Viazovska gave a miraculous solution to the sphere packing problem in 8 dimensions. One component of her proof is a proof of two inequalities that are satisfied by explicit modular forms. Her proof of those inequalities is heavily numerics-based and not as simple or conceptual as one might wish for. (To quote Henry Cohn who mentioned this issue in his paper A conceptual breakthrough in sphere packing: "Unfortunately, no simple proof of these inequalities is known at present, but one can verify them by reducing the problem to a finite calculation".) So, it seems desirable to find a more conceptual proof of these inequalities, or at least a proof that can be read and understood by a human without having to rely on a large amount of computer calculations.
I've now found what appears to be a relatively simple proof of Viazovska's inequalities that avoids the numerical work her original proof requires, but one part of my argument assumes the above conjecture. If the conjecture is proved, we'll have an unconditional proof.
 A: Not an answer but a possible approach. Set $T=\theta_3$ and
$F_2(q)=\sum_{m\ge0}\sigma(2m+1)q^{2m+1}$, so that $T^4$ and $F_2$ form a basis of $M_2(\Gamma_0(4))$. It is easy to see that $f\in M_{10}(\Gamma_0(4))$, and although this has dimension $6$ we have simply $f/2=T^{20}-40T^{16}F_2+640T^{12}F_2^2$.
Now $T^4-8F_2$ is a well-known Eisenstein series of weight 2 with positive coefficients, so the same is true for $T^{20}-8T^{16}F_2$. The first three
coefficients of the remaining part $32(20T^{12}F_2^2-T^{16}F_2)$ are negative,
but all the others seem to be positive. Maybe this is easier to prove.
Alternatively, one can write $f$ as a linear combination of $E_{10}$, $E_{10}(2\tau)$, and $E_{10}(4\tau)$ plus a cusp form; the Eisenstein part must have positive coefficients and trivially proved since the coefficients are explicit, and the cuspidal part coefficients can be bounded by the Ramanujan--Deligne bound and shown to be smaller than the Eisenstein part. I have not written the details, but this should be possible.
A: Hmm, after a bit of googling I found that my conjecture is in fact proved in the 2018 Harvard honors thesis "Modular Magic" by Aaron Slipper. (The lucky google search that yielded this was "2607840 theta series" -- 2607840 being the coefficient of $q^6$ in $f(q)/2$ -- but it only occurred to me to try this after an evening spent googling and reading related literature.)
Funnily enough, Slipper proved the fact that I need as part of his exposition of Viazovska's work on sphere packing, and needed it for reasons that seem a bit related to my own reasons for needing it, but not exactly to prove Viazovska's inequalities in the simple fashion I discovered.
Slipper's argument, found in page 76 of his thesis (as a component of his proof of Proposition 4.4.6 from page 75), is quite elegant. In my notation, he found that the function $f$ can be represented as
$$ 
f(q) = \frac{2 \Theta(q) + \Theta(-q)}{3}
$$
where $\Theta(q)$ is the theta function associated with a certain 20-dimensional lattice, called the DualExtremal(20,2)a lattice. Since the Fourier coefficients of the theta series of a lattice are nonnegative by definition, we immediately get the claim about the nonnegativity of the coefficients of $f(q)$.
Now, this is already very nice, but I'm wondering if this argument can be reworked to get a direct proof of the nonnegativity from simple properties of the Jacobi thetanull functions. In particular, I worked out that the function $\Theta(q)$ can be represented as
$$
\Theta = 2X^5 - 5 X^4 Y + 5 X^3 Y^2 + 5 X^2 Y^3 - 5 X Y^4 + 2Y^5,
$$
where $X=\theta_3^4$, $Y=\theta_2^4$. So maybe from here it's not too hard to show that $\Theta$ has nonnegative Fourier coefficients. I'll think more about this tomorrow. Edit: yes, I was right about that -- see below.

Edit: here is a self-contained proof of my nonnegativity conjecture, inspired by the argument from Slipper's thesis but reworked to give a more direct proof that does not rely on facts about lattices and their theta series:

*

*Denote $X=\theta_3^4$, $Y=\theta_2^4$ as above.


*Because of the well-known identity $\theta_3^4=\theta_2^4+\theta_4^4$, $f(q)$ can be written as
\begin{align*}
f &= X^3 Y^2 + X^2 Y^3 + X^3 (X-Y)^2 + X^2 (X-Y)^3
\\ &= 2 X^5 - 5 X^4 Y + 5 X^3 Y^2
\end{align*}


*Define a function $\Theta(q)$ by
$$
\Theta(q) = 2X^5 - 5 X^4 Y + 5 X^3 Y^2 + 5 X^2 Y^3 - 5 X Y^4 + 2Y^5.
$$


*Observe that, since $X(-q)=X-Y$ and $Y(-q)=-Y$, we have
\begin{align*}
\Theta(-q) &= 2 (X-Y)^5 - 5 (X-Y)^4 (-Y) + 5 (X-Y)^3 (-Y)^2 
\\ & \qquad + 5 (X-Y)^2 (-Y)^3 - 5 (X-Y) (-Y)^4 + 2 (-Y)^5
\\ &=
2 X^5 - 5 X^4 Y + 5 X^3 Y^2 - 10 X^2 Y^3 + 10 X Y^4 - 4 Y^5.
\end{align*}


*From the last two equations it follows that
$$
f(q) = \frac{2\Theta(q) + \Theta(-q)}{3}.
$$


*Note that an alternative way to express $\Theta$ is as
$$
\Theta = \frac{(2X-Y)^5 + 30 (2X-Y)^2 Y^3 + 15 (2X-Y) Y^4 + 18 (2X-Y)^5}{16}
$$


*Note that both $Y$ and $2X-Y$ have nonnegative Fourier coefficients. For $Y$ this is immediate from the definition, and for $2X-Y$ this follows from the identity
$$
\frac{2X-Y}{2} = 1 + 24 \sum_{n=1}^\infty \sigma_{\textrm{odd}}(n) q^{2n},
\qquad (*)
$$
where $\sigma_{\textrm{odd}}(n)$ is the odd divisor function, defined by
$$
\sigma_{\textrm{odd}}(n) = \sum_{d\,|\,n,\ \  d\textrm{ odd}} d.
$$


*We expressed $\Theta$ as a linear combination with positive coefficients of monomials in $Y$ and $2X-Y$, which have nonnegative Fourier coefficients. Therefore $\Theta$ has nonnegative Fourier coefficients.


*By the way we expressed $f$ in terms of $\Theta$, it's clear that $f$ too must have nonnegative Fourier coefficients. QED.
Remarks:

*

*The identity $(*)$ is standard, see for example this OEIS page. (I haven't checked, but I think it can be proved using the identities in this MathOverflow post expressing $\theta_2^4$ and $\theta_3^4$ in terms of Eisenstein series.)


*The modular form $2X-Y$ is twice the Eisenstein series of weight 2 referred to in Henri Cohen's answer, where it was written as $T^4-8F_2$. (Thanks Henri for the suggestion to look for a representation of $f$ involving this function!)


*The function I'm denoting by $\Theta$ is the same as the function $\Theta$ from page 76 of Slipper's thesis. But I defined it in terms of the Jacobi thetanulls $\theta_2$, $\theta_3$, and he defined his $\Theta$ as the theta series of a certain (pretty esoteric as far as I can tell) 20-dimensional lattice. Proving that the two functions coincide would require additional work, but isn't necessary to make my version of the nonnegativity proof work.


*The proof above does not explain how I found the representation of $\Theta$ as a linear combination with positive coefficients of monomials in $2X-Y$ and $Y$. I did this by solving a linear program in 21 variables (using Mathematica's Minimize command), but maybe there's an easier way.
