Can a bivariate polynomial evaluated modulo 9 produce a valid Sudoku grid?

Question

Can a bivariate polynomial with integer coefficients $f\colon \mathbb{Z}\times\mathbb{Z}\to\mathbb{Z}$ generate a valid $9\times 9$ Sudoku grid, when evaluated modulo 9?

A consequence is that for every $a=1,2,\ldots,9$, both $x\mapsto f\left(x,a\right)$ and $y\mapsto f\left(a,y\right)$ are permutation polynomials $\pmod{9}$. However, this does not capture the requirement of no repeating digits in a subsquare.

---

Random search found two Latin squares (in both cases, there are repeating digits in a subsquare):

The polynomial $8\cdot y^3 + 5\cdot y + 7\cdot x^3 + 3\cdot x^2 + 7\cdot x$ generates a Latin square when evaluated modulo 9:

(%i41) f(x, y) := 8 * y^3 + 5 * y + 7 * x^3 + 3 * x^2 + 7 * x;
                                 3            3      2
(%o41)             f(x, y) := 8 y  + 5 y + 7 x  + 3 x  + 7 x
(%i42) f(x,y);
                           3            3      2
(%o42)                  8 y  + 5 y + 7 x  + 3 x  + 7 x
(%i43) for y: 1 thru 9 do print(makelist(mod(f(x,y), 9), x, 1, 9));
[3, 5, 7, 6, 8, 1, 0, 2, 4] 
[1, 3, 5, 4, 6, 8, 7, 0, 2] 
[5, 7, 0, 8, 1, 3, 2, 4, 6] 
[0, 2, 4, 3, 5, 7, 6, 8, 1] 
[7, 0, 2, 1, 3, 5, 4, 6, 8] 
[2, 4, 6, 5, 7, 0, 8, 1, 3] 
[6, 8, 1, 0, 2, 4, 3, 5, 7] 
[4, 6, 8, 7, 0, 2, 1, 3, 5] 
[8, 1, 3, 2, 4, 6, 5, 7, 0] 
(%o43)                               done

Likewise $6\cdot y^3+6\cdot x\cdot y^2+6\cdot y^2+3\cdot x^2\cdot y+3\cdot x\cdot y+5\cdot y+5\cdot x^3+2\cdot x$:

(%i45) f(x,y);
              3        2      2      2                      3
(%o45)     6 y  + 6 x y  + 6 y  + 3 x  y + 3 x y + 5 y + 5 x  + 2 x
(%i46) for y: 1 thru 9 do print(makelist(mod(f(x,y), 9), x, 1, 9));
[0, 1, 5, 6, 7, 2, 3, 4, 8] 
[8, 3, 7, 5, 0, 4, 2, 6, 1] 
[4, 5, 3, 1, 2, 0, 7, 8, 6] 
[6, 7, 2, 3, 4, 8, 0, 1, 5] 
[5, 0, 4, 2, 6, 1, 8, 3, 7] 
[1, 2, 0, 7, 8, 6, 4, 5, 3] 
[3, 4, 8, 0, 1, 5, 6, 7, 2] 
[2, 6, 1, 8, 3, 7, 5, 0, 4] 
[7, 8, 6, 4, 5, 3, 1, 2, 0] 
(%o46)                               done

Answer 1

There are no such polynomials.

Note that the polynomial $h(x)=(x-1)(x-2)(x-3)$ is always zero modulo $3$ for integral $x$, and (of course) similarly for $h(y)=(y-1)(y-2)(y-3)$. One can certainly write any integral polynomial $P(x,y)$ as

$$F(x,y)=A(x,y)+h(x)R(x,y)+h(y) S(x,y),$$

where $A(x,y)$ has degree at most $2$ in $x$ and $y$, and where $R$ and $S$ are integral, just by polynomial division. Now note that in the first $3 \times 3$ square with $x=1,2,3$ and $y=1,2,3$ we have $h(x)=h(y)=0$, and so

$$P(x,y) \equiv A(x,y) \bmod 9, A(x,y) = a+bx+cy+dxy+ex^2+fy^2$$

in this range. By the Sudoku property, the values in this square must be all the residues mod $9$. But additionally, we have:

Lemma: In the first $3 \times 3$ square, the values mod $3$ in each column and each row are $0,1,2 \bmod 3$ precisely once in some order.

Proof: WLOG, suppose there was a duplication in some row, with $x$ fixed. Since $P(x,y)$ is integral, $P(x,y+3i) \equiv P(x,y) \bmod 3$. That would imply there would be six values in the same row which are congruent mod $3$, and thus they can not all be distinct mod $9$.

It's now simple enough to loop through the values of $a,b,c,d,e,f \bmod 9$ to determine when the first square has distinct entries mod $9$ for all the $9$ values, and also distinct values mod $3$ in every row and column. There are some simplifications of course, the constant $a$ makes no difference, and similarly we can also scale the entire equation by a unit mod $9$ and this also makes no difference. The conclusion of the Lemma already implies that $$A(x,y)=a+bx+cy \mod 3$$ with $bc \not\equiv 0 \bmod 3$. One can prove this by a calculation, but (added) more theoretically as follows:

Lemma: A quadratic $q(x) = r+sx+tx^2 \bmod 3$ takes on $3$-distinct values for $x=-1,0,1 \bmod 3$ if and only if $t=0$ and $s \ne 0$.

Proof: If $t=0$ the result is clear, so it suffices to prove $t=0$. If the values are distinct then they are $0,1,2 \bmod 3$ in some order. Thus

$$0 \equiv 0+1+2 \equiv q(-1)+q(0)+q(1) = (r-s+t) + (r) + (r+s+t) = 2t \bmod 3,$$

which implies that $t=0$.

Now returning to $A(x,y)=a+bx+cy+dxy+ex^2+fy^2$, we know specializations with $x$ or $y$ fixed must give polynomials which take on distinct values mod $3$ by our lemma. Thus:

1. If you specialize to $y=0$ you get $q(x)=a+bx+ex^2$ and then $e \equiv 0 \bmod 3$ and $b \not\equiv 0 \bmod 3$;
2. If you specialize to $x=0$ you get $q(y)=a+cy+fy^2$ and then get $f \equiv 0 \bmod 3$ and $c \not\equiv 0 \bmod 3$.
3. If you take $y=1$ or $y=-1$ you get $q(x)=a+bx+c+dx \bmod 3$ or $q(x)=a+bx+c-dx \bmod 3$, and you deduce that $b,b+d,b-d \bmod 3$ are all non-zero. But if $d \not\equiv 0 \bmod 3$ then these are all distinct and hence not all non-zero, so $d \equiv 0 \bmod 3$.

Hence we have: $$A(x,y)=a+bx+cy+3(dxy + e x^2 + f y^2)$$

which (with $a=0$ and scaling $b$ to $1$ since it is a unit) leaves only $6 \cdot 9^3 = 4374$ cases to check. The conclusion of this computation however is that there are no such solutions. Hence the answer to the question is negative.