$\DeclareMathOperator\u{\mathbf u}$$\DeclareMathOperator\e{\mathbf e}$I think I know how to answer my question. First, note that $p(x,y) = x^{-1}y^{-1} \tilde p(x,y)$, where $\tilde p(x,y) = (x-1)^2y+x(y-1)^2 \in \mathbb Z[x,y]$. More generally, any $q(x,y)\in\mathbb Z[x,x^{-1},y,y^{-1}]$ can be written as $q(x,y) = x^a y^b \tilde q(x,y)$ for some $a,b\in\mathbb Z$ such that $\tilde q(x,y)\in\mathbb Z[x,y]$. So, my question can be rephrased as follows:
Is there a $\tilde q(x,y)\in\mathbb Z[x,y]$ that is not divisible by $\tilde p(x,y)$ such that $\tilde q(x,y) \hat f(x,y)=0$ for any $f\in\mathcal F$?
I will show that the answer is no, i.e., any $\tilde q(x,y)$ that satisfies $\tilde q(x,y) \hat f(x,y)=0$ for any $f\in\mathcal F$ must be divisible by $\tilde p(x,y)$.
Choosing a lexicographic monomial order with $x\succ y$ for multivariate division, any $\tilde q(x,y)$ can be uniquely written as
$$
\tilde q(x,y) = x^2\alpha(x) + x \beta(y) + \gamma(y) + \tilde r(x,y)\tilde p(x,y)~,\tag{1}\label{1}
$$
where $\tilde r(x,y)\in\mathbb Z[x,y]$ and $\alpha(X),\beta(X),\gamma(X)\in\mathbb Z[X]$. Let $u$, $v$, and $w$ be the degrees of $\alpha(X)$, $\beta(X)$, and $\gamma(X)$ respectively, and let
$$
\alpha(X) = \sum_{i=0}^u a_i X^i~,\qquad \beta(X)=\sum_{j=0}^v b_j X^j~,\qquad \gamma(X)=\sum_{k=0}^w c_k X^k~.
$$
I will show that $a_i = b_j = c_k = 0$ if $\tilde q(x,y) \hat f(x,y)=0$ for all $f\in\mathcal F$.
Consider the function
$$
f(m,n;t) = (-1)^m t^{-m+n} \left(\frac{1+t}{1-t}\right)^{-m-n}~.
$$
It is easy to check that it is a discrete harmonic function on a square lattice for any $t\ne 0,\pm 1$ (see page 13 of this note for a related function). Now, using \eqref{1}, the equation $\tilde q(x,y) \hat f(x,y;t) = 0$ becomes
$$
x_t^2 \alpha(x_t) + x_t \beta(y_t) + \gamma(y_t) = 0~,
$$
where $x_t := t\left(\frac{1+t}{1-t}\right)$ and $y_t := t\left(\frac{1-t}{1+t}\right)$. Let $\nu := \max(v-1,w)$. Multiplying by $(1-t)^{u+2} (1+t)^\nu$ gives a polynomial in $t$ given by
$$
\sum_{i=0}^u a_i t^{i+2} (1+t)^{\nu+i+2} (1-t)^{u-i} + \sum_{j=0}^v b_j t^{j+1} (1+t)^{\nu-j+1} (1-t)^{u+j+1} + \sum_{k=0}^w c_k t^k (1+t)^{\nu-k} (1-t)^{u+k+2} = 0~.\tag{2}\label{2}
$$
Since this equation holds for any $t\ne0,\pm 1$, the polynomial in $t$ must vanish identically (even at $t=0,\pm1$).
We can write \eqref{2} as
$$
\sum_{i=0}^u a_i \alpha_i(t;\u) + \sum_{j=0}^v b_j \beta_j(t;\u) + \sum_{k=0}^w c_k \gamma_k(t;\u) = 0~.\tag{3}\label{3}
$$
where $\u:=(u,v,w)$. I will proceed by induction to show that the set $S(\u)$ of polynomials $\alpha_i(t;\u)$'s, $\beta_j(t;\u)$'s, and $\gamma_k(t;\u)$'s is linearly independent for any $\u$:
Base cases: For $\u=\mathbf 0$, $S(\mathbf 0)$ contains
$$
\alpha_0(t;\mathbf 0) = t^2(1+t)^2~, \qquad \beta_0(t;\mathbf 0) = t(1-t^2)~,\qquad \gamma_0(t;\mathbf 0) = (1-t)^2~,
$$
which are clearly linearly independent. Similarly, for $\u=(0,1,0)$, $S(0,1,0)$ contains
$$
\alpha_0(t;\u) = t^2(1+t)^2~, \qquad \beta_0(t;\u) = t(1-t^2)~, \qquad \gamma_0(t;\u) = (1-t)^2~,
\\
\beta_1(t;\u) = t^2(1-t)^2~,
$$
which are also linearly independent.
Induction step: Assume that $S(\u)$ is linearly independent.
For $\u+\e_1$, where $\e_1:=(1,0,0)$, $S(\u+\e_1)$ contains
$$
\alpha_i(t;\u+\e_1) = \begin{cases}
(1-t)\alpha_i(t;\u)~,&\text{for } i = 0,\ldots,u~,
\\
t^{u+3} (1+t)^{\nu+u+3}~,&\text{for } i=u+1~,
\end{cases}
\\
\beta_j(t;\u+\e_1) = (1-t) \beta_j(t;\u)~,\quad\text{for } j = 0,\ldots,v~,
\\
\gamma_k(t;\u+\e_1) = (1-t) \gamma_k(t;\u)~,\quad\text{for } k = 0,\ldots,w~.
$$
Except for $\alpha_{u+1}(t;\u+\e_1)$, all the other polynomials are linearly independent because they are all $(1-t)$ times polynomials that are linearly independent by hypothesis. On the other hand, $\alpha_{u+1}(1;\u+\e_1)\ne0$, so it is linearly independent of the other polynomials, which vanish at $t=1$. Therefore, $S(\u+\e_1)$ is linearly independent.
For $\u+\e_2$ and $\u+\e_3$, the above argument is a bit subtle. Instead, we can proceed as follows. Say $S(u,v_0,v_0-1)$ is a linearly independent set. Then $S(u,v_0,w)$ with $w<v_0-1$, and $S(u,v,v_0-1)$ with $v<v_0$ are also linearly independent sets because $S(u,v,v_0-1)\subset S(u,v_0,v_0-1)$ and $S(u,v_0,w)\subset S(u,v_0,v_0-1)$. So, all I need to show is that if $S(u,v_0,v_0-1)$ is linearly independent, then $S(u,v_0+1,v_0)$ is also linearly independent. This can be shown using an argument similar to the one above (it's only slightly more complicated).
Since the polynomials $\alpha_i(t;\u)$'s, $\beta_j(t;\u)$'s, and $\gamma_k(t;\u)$'s are all linearly independent for any $\u$, it follows from \eqref{3} that $a_i = b_j = c_k = 0$. Therefore, $\tilde q(x,y)$ is divisible by $\tilde p(x,y)$.
