"Spherical wave" in a 2D square lattice of masses and springs Assume a 2D square array of masses with mass $m$ connected by springs with constant $\kappa$. Then the equation of motion for small perturbations in the plane of the array to the mass at $(l,n)$ is
$$
\ddot{\vec{\psi}}_{l,n} = \Omega^2(\vec{\psi}_{l+1,n}-2\vec{\psi}_{l,n} + \vec{\psi}_{l-1,n})+\Omega^2(\vec{\psi}_{l,n+1}-2\vec{\psi}_{l,n} + \vec{\psi}_{l,n-1})
$$
where $\Omega^2=\kappa/m$.
How will a "spherical wave" will look like in this system? I.e. say I take the mass at $(0,0)$ and make it oscillate along the $y$ axis (i.e. along the $n$ direction). What will the disturbance look like as a function of $(l,n)$?
If I had to guess, it would be something like $e^{ik(|n|+|l|)}$, however I have failed to prove it. Also, how will it decay in space?
EDIT
Here is an attempt. In such a scenario, only the $y$ coordinate of each mass will oscillate, so we are reduced to a scalar wave. From symmetry we can focus on the region $l,n>0$. If $\psi(l,n,t)=A(l,n)e^{-i\omega t}$, then we must have $A(l+1,n)=A(l,n+1),\ \ A(l-1,n)=A(l,n-1)$. Also, $A(l,n)=A(l+n)$. Let $N=l+n$. Then plugging the solution in, we get
$$
-\omega ^2A(N)=2\Omega^2\left(A(N+1)-2A(N)+A(N-1)\right)
$$
which has a solution $A(N) = e^{ikN}$ when $\omega^2=4\Omega^2\sin^2\left(k/2\right)$.
However, this solution does not decay, which makes me skeptical.
 A: 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.
