This is a tough question, and I don’t think there’s a definitive answer yet. For some mathematical details, see the following survey articles:
https://arxiv.org/abs/1611.01685
https://arxiv.org/abs/1603.05202
Instead, I’ll focus on the big picture here. Why modular forms? I can see a couple of potential answers:
(1) Why not modular forms? Before Viazovska’s proof, numerical experiments indicated that there were remarkable special functions in 8 and 24 dimensions that would prove the optimality of $E_8$ and the Leech lattice. However, nobody had any idea how to construct them explicitly, or prove existence at all. Modular forms are by far the most important class of special functions related to lattices (in higher dimensions, since arguably trig functions and the exponential function are the most important special functions related to lattices), so lots of people expected that the magic functions for sphere packing should be connected somehow to modular forms. The proof had to wait for Viazovska to come up with a beautiful integral transform, but the fact that it used modular forms wasn’t such a great surprise. I.e., her contribution wasn’t the idea that modular forms should play a role, but rather figuring out how to use them, which was quite subtle and ingenious.
You’re right that nobody has any idea how to use modular forms to optimize the linear programming bound in other dimensions. However, it’s possible that they will continue to play a role. For example, see the example Felipe Gonçalves and I found at the end of Section 2.1 of our paper https://arxiv.org/abs/1712.04438 (which is not a sphere packing bound, but closely related). It really looks like a small perturbation of a function based on modular forms (see https://arxiv.org/abs/1903.05737), and I wouldn’t be surprised if the optimal function has a nice series expansion based in some way on modular forms. From this perspective, the remarkable thing about 8 and 24 dimensions wouldn’t be the appearance of modular forms, but rather the fact that the series collapses to a single term, with a matching sphere packing. However, this is all speculative.
(2) The other perspective is that we have very little understanding of why 8 and 24 dimensions are special in the first place. For example, why shouldn’t sphere packing in 137 dimensions also admit an exact solution via linear programming bounds? It sure doesn’t look like it does, but perhaps we just don’t know the right sphere packing to use, and some currently unknown packing might match the upper bound. That would be very surprising, since our experience is that exceptional phenomena occur in clusters. We would expect to see some sort of remarkable symmetry group, probably a finite simple group, and there aren’t any candidates acting on 137 dimensions. However, this expectation is just based on our limited experience, and mathematics can confound our expectations. So far, nobody has found even a convincing heuristic argument for why there shouldn’t be an exact solution in 137 dimensions, and that’s a major gap in our understanding. The most we can say is that it would have to differ in some important ways from 8 and 24 dimensions, which is far from an explanation of why it can’t happen.
I guess I’d summarize it like this. If you accept that lattices in 8 and 24 dimensions play a special role, then modular forms feel naturally connected. However, we’re missing a deeper explanation of the role of these special dimensions.
Let me add one more specific mathematical comment. The magic functions in 8 and 24 dimensions fit into a general picture of building radial functions that vanish on all but finitely many vector lengths in a lattice, and whose Fourier transforms vanish on all but finitely many vector lengths in the dual lattice. If you can do this in full generality, then Poisson summation lets you solve for the number of lattice vectors of each length. These are the coefficients of a modular form, namely the theta series of $E_8$ or the Leech lattice, so the conclusion is that this family of functions somehow “knows about” the theta series. In other words, you can’t expect to construct the whole family without running into modular forms in some way. This leaves a couple of possibilities: maybe the magic functions for sphere packing are simpler than most functions in this family, and could be constructed without modular forms, or maybe these functions are deeper than modular forms (and require some mysterious special functions not yet known to mathematicians). What we know now is that modular forms suffice, and in a sense are necessary because the magic functions are unique.