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 \tag{$*$} $$

(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. For our single equation in two variables, 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 - \text{number of solutions to $(*)$ over $\mathbb{F}_p$}$$

Question: When is it true that the numbers $\{a_p\}$ are "modular", in the  sense that there exists a modular form 

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

such that $b_p = a_p$ for almost all primes $p$? (The "almost all" is to avoid problems with bad primes. Note that $f$ is still uniquely determined by the above requirement.)

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.)