Any solution.

Very similar to https://oeis.org/A071053. This is a cellular automaton. Here are the rules. All lattice points can either be activated(on), or deactivated(off). A neighbor $(x,y)$ of a lattice point $(m,n)$ is a different point such that $|x-m|, |y-n| \leq 1$. To go from step $m$ to $m+1$, the following happens: If a cell is deactivated, it will be activated if it has an odd number of neighbors. If a cell is activated, it will be deactivated if it has an even number of neighbors. Only the point $(0, 0)$ is activated on step $0$. Unless this has a nice pattern, it isn't gonna be easy. Not sure how to count this... (It's not exactly the right picture due to some weird integer offset)

Attachments:

Trivially $a_n$ equals the number of non-zero terms of $f_n(x,y)=[(x+\frac{1}{x}+1)(y+\frac{1}{y}+1)-1]^n$ in $F_2$. By manipulting functions, we have $\forall k,t\in N$, $0\le i\le 2^{k-1}-1$, $a_{t\cdot 2^k+i}=a_{t}\cdot a_{i}$. This means that if there is a $0$ in $n$'s binary expression, we could cut the expression in two parts and compute them instead. (For example, $a_{(101101110011)_2}=a_{(1011)_2}\cdot a_{(1110011)_2}$.) Thus the problem reduces to finding $a_{2^k-1}$ for all $k\in N$. Since $f_{2^{k+1}-1}=[(x^{2^k}+\frac{1}{x^{2^k}}+1)(y^{2^k}+\frac{1}{y^{2^k}}+1)-1]\cdot f_{2^k-1}$ (in $F_2$ of course), the value matrix of $f_{2^{k+1}-1}$ is equal to the sum of some $8$ tranlations of that of $f_{2^{k}-1}$. A (not) simple reduction would show that $a_{2^k-1}=\frac{5\cdot 4^k+(-2)^{k+1}}{3}$ (it satisfies a second-order linear recurrence with constant coefficients). Combining this with the above claim and we get the desired result. (Unless I'm very mistaken, this is not a closed form in terms of $n$, since it might need arbitrarily many steps of calculation. Anyway, I believe it's the simplest expression- if I didn't get wrong in the induction step.)