Suppose we are given some polynomial with integer coefficients, which we regard as carving out an affine variety $E$, for example:

$ 3x^2y - 12 x^3y^5 + 27y^9 - 2 = 0 \quad \quad (*) $

(We might consider a bunch of equations, we might work over projective space, but let's keep it simple for now).

We are interested in the number of points on $E$ when we reduce modulo $p$, i.e. over the finite field $\mathbb{F}_p$ as the prime $p$ varies. As a rough approximation, we would expect $p$ points in general. So for each prime $p$, we define numbers $a_p$ which measure the deviation from this,

$a_p := p - \mbox{number of solutions to (*) over $\mathbb{F}_p$}$

Question: When is it true that there exists a modular form

$f = \sum_{n=1}^\infty b_n q^n$

such that $b_p = a_p$ for all primes $p$? (Note that this would determine $f$ uniquely.)

The Modularity Theorem of Breuil, Conrad, Taylor and Diamond says that this is true when $E$ is an elliptic curve, i.e. takes the form $y^2 = 4x^3 - g_2 x - g_3$ for some integers $g_2, g_3$. In that case, $f$ is a weight 2 modular form of level $N$ where $N$ is the "conductor" of $E$.

But is it true for more general varieties?

(Note: I am aware of a generalized "Modularity Theorem" for certain Abelian varieties, but it's not clear to me that what people mean by "Modular" in that context is the same as the simple-minded notion I'm using --- that an adjusted count of points mod $p$ gives the Fourier coefficients of a modular form.)