There's no such family of functions for $\epsilon = 0$, and even for $\epsilon_n$ depending on $n$ decaying at a rate to be determined later. I don't know what happens if $\epsilon > 0$ is independent of $n$.
For a multiset of indices $I = \{ i_1, \dots i_n \}$ write $a_I = \prod_{i \in I} a_i$ for the corresponding monomial. We're looking for a family of functions $f_n$ which can be written
$$f_n(a_1, \dots a_n) = \sum_{i=1}^N f_{n, i}(a_{I_{n, i}})$$
for some functions $f_{n, i}$ and some indices $I_{n, i}$, where $N = O(n^d)$ for some absolute constant $d$, such that $f_n(a_1, \dots a_n) = 1$ if $\sum a_i = 0$ and $|f_n(a_1, \dots a_n)| < \epsilon_n$ otherwise.
We'll restrict our attention to the Hamming cube $(a_1, \dots a_n) \in \{ -1, 1 \}^n$ and identify a point in this cube with the subset $S(a) = \{ i : a_i = -1 \}$ of indices with coordinate $-1$. On this cube the monomials $a_I$ behave very simply. Since $a_i^2 = 1$ for $a_i \in \{ -1, 1 \}$ we can restrict our attention to sets of indices $I$ rather than multisets. Then we have
$$a_I = (-1)^{|S(a) \cap I|}$$
hence that $f_{n, i}(a_{I_{n, i}})$ takes on one of two possible values depending on the parity of $|S(a) \cap I|$, and in fact we always just have
$$f(a_I) = \frac{f(1) + f(-1)}{2} + \frac{f(1) - f(-1)}{2} a_I$$
(since this identity holds for $a_I = 1$ or $a_I = -1$). It follows that any function of the form $\sum_{i=1}^N f_{n,i}(a_{I_{n,i}})$, when restricted to the Hamming cube, is a constant plus a linear combination of $N$ functions of the form $a_I$.
The functions $a_I$ are precisely the characters of $\{ -1, 1 \}^n$ regarded as a finite abelian group under pointwise multiplication, and now standard discrete Fourier analysis (the Hadamard transform) tells us that any function $f : \{ -1, 1 \}^n \to \mathbb{R}$ on the Hamming cube can be represented uniquely as a linear combination of the $a_I$, with coefficients determined by the inner products
$$\langle f, a_I \rangle = \frac{1}{2^n} \sum_{(a_1, \dots a_n) \in \{ -1, 1 \}^n} f(a_1, \dots a_n) (-1)^{|S(a) \cap I|}.$$
Our strategy from here will be to show that there are infinitely many $n$ such that exponentially many of the Fourier coefficients of $f_n$ don't vanish. The points in the Hamming cube satisfying $\sum a_i = 0$ are exactly the ${n \choose \frac n 2}$ points with the same number of $1$s as $-1$s, and in particular there are no such points if $n$ is odd; from here on we assume that $n$ is even. We have
$$\langle f_n, a_I \rangle = \frac{1}{2^n} \left( \sum_{J \subseteq \{ 1, 2, \dots n \} : |J| = \frac n 2} (-1)^{|J \cap I|} + O(2^n \epsilon) \right).$$
If we replace $I$ by its complement the sum $\sum_{J \subseteq \{ 1, 2, \dots n \} : |J| = n/2} (-1)^{|J \cap I|}$ is multiplied by $(-1)^n$, so to analyze this sum we can assume WLOG that $|I| \le \frac n 2$. By breaking it up depending on the possible values of $|J \cap I|$ we find that
$$\sum_{J \subseteq \{ 1, 2, \dots n \} : |J| = \frac n 2} (-1)^{|J \cap I|} = \sum_{k=0}^{|I|} (-1)^k {|I| \choose k} {n-|I| \choose \frac n 2 - k}.$$
This is the coefficient of $x^{\frac n 2}$ in $(1 - x)^{|I|} (1 + x)^{n-|I|}$, which is a little messy. This product simplifies as much as possible if $|I| = \frac n 2$, where it becomes $(1 - x^2)^{\frac n 2}$, which gives
$$\sum_{k=0}^{\frac n 2} (-1)^k {\frac n 2 \choose k} {\frac n 2 \choose \frac n 2 - k} = \begin{cases} 0 \text{ if } \frac n 2 \equiv 1 \bmod 2 \\ (-1)^{\frac n 4} {\frac n 2 \choose \frac n 4} \text{ if } \frac n 2 \equiv 0 \bmod 2 \end{cases}$$
So from here we assume that $n = 4m$ is divisible by $4$. In this case we've shown that if $\epsilon_{4m} = 0$ then at least ${4m \choose 2m} \sim \frac{2^{4m}}{\sqrt{2\pi m}}$ (by e.g. Stirling's formula) of the Fourier coefficients don't vanish, which is exponentially many in $n$, and in particular more than $N = O(n^d)$ for any $d$. We still get the same conclusion as long as $2^{4m} \epsilon_{4m} < {2m \choose m} \sim \frac{2^{2m}}{\sqrt{\pi m}}$, hence as long as
$$\epsilon_n = o \left( \frac{1}{2^{ \frac n 2} \sqrt{n}} \right).$$
If we only want to rule out $N = O(1)$ we can instead look at the ${n \choose 2}$ Fourier coefficients where $|I| = 2$; these are all equal to
$${n-2 \choose \frac n 2} - 2 {n-2 \choose \frac n 2 - 1} + {n-2 \choose \frac n 2 - 2}$$
which, after some simplifications I'll save you (I'm a little surprised I don't know an elegant way to do this), is equal to
$$\frac{(2n-1)}{(\frac n 2)(\frac n 2 + 1)} {n-2 \choose \frac n 2 -1} = \Theta \left( \frac{2^n}{n^{\frac 3 2}} \right)$$
so to guarantee that at least ${n \choose 2}$ of the Fourier coefficients are nonzero, which even rules out $N = O(n)$, it suffices to take
$$\epsilon_n = o \left( n^{- \frac 3 2} \right).$$
Similar computations are available for other small values of $|I|$ and using more of these it ought to be possible to improve the bounds.
(Also, if we only want the conclusion for $\epsilon = 0$ there's a similar but much simpler argument given by restricting to the cube $\{ 0, 1 \}^n$ instead and observing that $f_n(a_1, \dots a_n) = \prod_{i=1}^n (1 - a_i)$ on this cube, showing that the monomials $a_I$ for $I$ a set (which are all we need to consider on $\{ 0, 1 \}^n$) are linearly independent, and arguing that this implies that $N \ge 2^n - 1$, for all $n$. But I don't know how to adapt this argument to the $\epsilon_n \neq 0$ case.)
