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.