What you are seeing on the square plate are the resonant modes of the structure. Each of these modes has a particular frequency associated with it, and is rung up when the plate is driven at that frequency. These resonant modes act like standing waves on a string: where some parts of the plate are moving a lot while other parts are standing still. The sand is bounced away from the parts which are moving a lot (the nodes) and is left in the places where the plate is not moving at all (the anti-nodes).
I'll address your specific questions in more detail in reverse order.
what importance has the source point of the vibration?
Almost none (with a few caveats). If you drive the plate at one of its nodes (a point where it will be moving a lot), then you will excite the resonant mode more than if you drive it at a point which is near an anti-node. This will only effect the amplitude though, not the shape of the mode.
One caveat is that if the structure has a lot of internal friction (damping), then the mode structure will change slightly and the height of the nodes will die out as you get far from the source point.
Are the patterns created related to the Laplace eigenfunctions of the shape?
I'm not sure what you mean by Laplace eigenfunctions. As described in the section below, the solutions are eigenfunctions of the two dimensional wave equation. The starting point for the solution begins with what may be called the Fourier eigenfunctions, but the boundary conditions (which are complicated in the case of a sheet which is free to move at the ends) change these solutions into something completely different.
What is the theory behind this fact?
This journal article has a systematic derivation of exactly the situation you are asking about. This note has a more easily parsed derivation of the similar situation of a stretched rectangular membrane. I can walk through the basics of the membrane case here, and point out where it is different from your case where necessary.
If you consider an elastic material and write down the forces on an infinitesimal piece of the material as a function of its height, you find that the equation describing the height of any piece is the two dimensional wave equation;
$$
\frac{\partial^2 u(x,y,t)}{\partial t^2}=\frac{T}{\rho}\left(\frac{\partial^2 u(x,y,t)}{\partial x^2} + \frac{\partial^2 u(x,y,t)}{\partial y^2}\right),
$$
where $T$ is the tension of the surface (units of force per length) and $\rho$ is the mass density (units of mass per area). A more simplified way of writing this equation which lends itself better to a solution by separation of variables is
$$
c^2\nabla^2 u=\frac{\partial^2}{\partial t^2}u
$$
where $c^2=\frac{T}{\rho}$. In the case of a rigid material such as your situation this coefficient is given by $c^2=\frac{m\omega^2}{D}$ where $D$ is the cylindrical stiffness of the material, $\omega$ is the resonant frequency, and $m$ is the mass.
Using separation of variables we can separate this equation into three independent one dimensional differential equations;
$$
\begin{align}
\frac{d^2}{dt^2}G(t)+(c\nu)^2G(t)&=0\\
\frac{d^2}{dx^2}H(x)+k^2H(x)&=0\\
\frac{d^2}{dy^2}Q(y)+p^2Q(y)&=0\\
p^2+k^2&=\nu^2,
\end{align}
$$
where $u(x,y,t)=G(t)H(x)Q(y)$. The solutions to these equation may be thought of as the Fourier eigenfunctions; i.e. $\sin$, $\cos$, $\sinh$, and $\cosh$.
As with all physical problems, we need to specify the boundary conditions in order to get physically meaningful solutions. As an example, a rubber sheet which is clamped at the edges has the boundary conditions,
$$
u(\pm a, y)=0\qquad\&\qquad u(x,\pm b)=0
$$
where $2a$ and $2b$ are the $x$ and $y$ dimensions of the sheet respectively. These conditions say that the sheet can't move up and down at the edges where it is clamped. This situation is rather straightforward to solve and is worked out in the note mentioned above.
You're situation is a bit more difficult because the edges are completely free to move. In this case you have to incorporate into the boundary conditions the fact that the plate is a rigid object which doesn't need to be under tension to be supported. Taking them from the article mentioned above, they are
$$
\begin{align}
\frac{\partial^2u}{\partial x^2}+\nu_p\frac{\partial^2 u}{\partial y^2}&=\frac{\partial^3u}{\partial x^3}+(2-\nu_p)\frac{\partial^3u}{\partial x\partial y^2}=0\qquad @ x=\pm a\\
\frac{\partial^2u}{\partial y^2}+\nu_p\frac{\partial^2 u}{\partial x^2}&=\frac{\partial^3u}{\partial y^3}+(2-\nu_p)\frac{\partial^3u}{\partial x^2\partial y}=0\qquad @ y=\pm b,
\end{align}
$$
where $\nu_p$ is a material parameter known as Poisson's ratio which describes how much a material expands in one direction when compressed in the other.
The nice thing about solution of the wave equation by separation of variables is that the solutions you find form a complete set of eigenfunctions for the problem. Furthermore, their eigenvalues tell you the resonant frequency of each particular mode. The image below shows the second, third, and fourth order solutions to your equation (taken from the referenced journal article). When you excite these modes with the sand on the plate, the sand is bounced away from the high points and stays at the un-moving points which are called nodal lines. These nodal lines are the places where the solution is still zero in the diagram below.