To see how the square lattice determines how a wave looks at large distance, let me copy from my notes how response of the slightly simpler problem of how the solution of a massive scalar wave equation decays on a two-dimensional square lattice. I expect that the gapless wave equation works the same way except that the exponential decay will be replaced by a power-law decay.
Start with the Green function for the square lattice Laplacian
$$
G({\bf n}, m^2)= \int_{-\pi}^{\pi} \frac{d^2 k}{(2\pi)^2}\frac{e^{-i{\bf k}\cdot
{\bf n} }}{m^2+ 2 \sum_1^d(1-\cos k_i)}.
$$
Suppose $\bf n$ is becoming large in a particular direction specified by
a unit vector $\bf e$:
$$
{\bf n}=r{\bf e}, \qquad {\bf e}^2=1.
$$
We expect $G$ to fall off exponentially at large $|r|$
$$
G(r{ \bf e}, m^2)\approx e^{-\kappa({\bf e)}|r|}.
$$
The problem is to compute the direction depend decay rate $\kappa({\bf e})$, the inverse correlation
length in the direction $\bf e$. We will see that $\kappa$ is anisotropic
at large $m^2$ but becomes isotropic as $m^2$ becomes small.
Define
$$
f(\xi)= \int_{-\infty}^\infty dr e^{-ir\xi}G(r{\bf e}, m^2).
$$
This quantity should be analytic in a neighborhood of the real $\xi$ axis and the
asymptotic behavior at large $r$ of $G(r{\bf e}, m^2)$ will be
determined by the nearest singularity to the real $\xi$ axis. The nearest
singularity will be on the imaginary axis if $G$ does not oscillate. For example,
a singularity at $\xi=i\zeta_0$ would make $G\propto
e^{-\zeta_0r}$. This is because
$$
\int_{-\infty}^{\infty} e^{ir\xi}\frac 1{\xi-i\zeta_0} \frac {d\xi}{2\pi i}
=e^{-\zeta_0r}, \quad \hbox{ for } r>0.
$$
Now
$$
f(\xi)=
\pi\int_{-\pi}^{\pi} \frac{d^2 k}{(2\pi)^2}\frac{\delta({\bf k}\cdot{\bf e}-\xi)}{m^2+ 2 \sum_1^d(1-\cos k_i)}
\nonumber\\ = \pi \int_{-\pi}^{\pi} \frac{d^2 k}{(2\pi)^2}\frac{\delta({\bf k}\cdot{\bf e}-\xi)}
{D({\bf k})}.
$$
For small $m^2$, the singularity is expected to be caused by zeros of
the denominator $D({\bf k})$ pinching the contour of integration. This is certainly
what happens in the continuum where
$$
f(i\zeta) =\int_{-\infty}^{\infty} dk \frac 1{k^2+m^2-\zeta^2}.
$$
The integrand has poles at $k=\pm i\sqrt{m^2-\zeta^2}$ and these pinch when $\zeta=m$.
The resultant singularity
makes the propagator fall off as $e^{-m|r|}$, a result we know already.
To locate the pinch we move the $k$ integral off the real axis by setting ${\bf k} = i{\bf
K}$. (The fact the $k$ integrals have real endpoints at $k=\pm \pi$ is not relevant to finding the pinch.)
Then a pinch at $\xi=i\zeta_0$ requires
$$
D({\bf K})=0,\nonumber\\
\frac{\partial D}{\partial {\bf K}}=0, \quad {\rm on}\quad {\bf K}\cdot
{\bf e}=\zeta_0.
$$
We can impose the constraint by means of a Lagrange multiplier
$$
\frac{\partial}{\partial {\bf K}}\left\{D({\bf K})-\lambda {\bf K}\cdot
{\bf e}\right\}=0.
$$
These equations have a simple geometric interpretation. The set of
points ${\bf K}\cdot {\bf e}=\zeta_0$ is a straight line perpendicular to
${\bf e}$ and at a distance $\zeta_0$ from the origin. Th equation above says that, if there is a pinch singularity at at
$i\zeta_0$, then this line must be tangent to the curve $D({\bf
K})=0$. A more direct way to see this is to note that the points of
intersection of the contour of integration ${\bf K}\cdot{\bf e}\zeta_0$
with $D({\bf K})=0$ are the location of the poles of the inegrand in
the complex ${\bf k}$ plane. Clearly they can only pinch if they are
coincident and the line is tangent.
In our case
$$
D({\bf K})= 4+m^2-2 (\cosh K_1+\cosh K_2).
$$
For large $m^2$ the curve $D({\bf K})=0$ is essentially a square with sides at $\pm \cosh^{-1}(4+m^2)/2$
A little geometry shows that
$$
\zeta_0 \approx (\cos\theta +\sin\theta)\cosh^{-1}(4+m^2)/2\nonumber\\
\approx (\cos\theta +\sin\theta)\ln (4+m^2).
$$
Therefore, for large $m^2$
$$
G(|r|{\bf e},m^2) \approx e^{-|r|\zeta_0} = \frac
1{(4+m^2)^{|n_1|+|n_2|}}.
$$
as you conjecture.
This should perhaps not be surprising. At large mass the Green function is dominated by the shortest
lattice path between $0$ and ${\bf n}=(n_1,n_2)$. There are many of these, but they all have the same length
length $|n_1|+|n_2|$.
As $m^2$ becomes smaller the corners of the curve $D({\bf K})=0$ round off,
tending to circularity at at small $m^2$.
It is easy to show that along any of the four axes
$$
\kappa= \zeta_0=\cosh^{-1}(1+\frac 12 m^2),
$$
while at $45^\circ$ we have
$$
\kappa=\zeta_0= \sqrt 2 \cosh^{-1}(1+\frac 14 m^2).
$$
Both expressions are equal to $m$ when $m^2$ is small, so the fall-off of the wave becomes circularly symmetric.
