\section{Introduction}
In \em The Grammar of Ornament, \em published in 1856, Owen Jones wrote \cite{jones}:
\begin{quote}
 Whenever any style of ornament commands universal admiration, it will always be found to be in accordance with the laws which regulate the distribution of forms in nature.
\end{quote}
In the case of crystals, the laws which regulate their shape are dictated by the crystallographic groups.

\subsection{Crystallographic groups}
A crystal or crystalline solid is a solid material whose constituents, such as atoms, molecules, or ions, are arranged in a highly ordered microscopic structure; for a two dimensional example see Figure \ref{fig:graph}.  Mathematically we can identify the locations of the constituents with the points of a lattice, also known as a \em crystal lattice.  \em  Every two dimensional crystal has a symmetry group which is a plane crystallographic group.  Such a group consists of isometries of the plane. There are three basic types of isometries:  translations, rotations and reflections.  These form a group under composition.  The patterns in Figure \ref{fig1} have symmetry groups which are \em plane crystallographic groups.  \em  These are subgroups of the group of isometries of the plane which are topologically discrete and contain two linearly independent translations.  Equivalently, a plane crystallographic group\footnote{These are also known as wallpaper groups.} is a co-compact subgroup of the group of isometries of the plane. A subgroup in this context is called co-compact if the quotient space $\mathbb{R}^2/\Gamma$ by the subgroup, $\Gamma$, is compact.  The classification of these groups, up to equivalence, was achieved at the end of the 19th century by E. S. Fedorov \cite{fed} and A. Schoenflies \cite{schoen}.  Two groups are equivalent if they are equal up to a translation. In two dimensions, up to this notion of equivalence, there are seventeen crystallographic groups.

\begin{figure} \centering \includegraphics[width=8cm]{graphen} \caption{Graphene is an allotrope of carbon in the form of a two-dimensional, atomic-scale hexagonal lattice such that each point in the lattice corresponds to an atom.  This image is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license at \url{https://commons.wikimedia.org/wiki/File:Graphen.jpg}.} \label{fig:graph}  \end{figure}

One can also consider crystals in three dimensions, and mathematically we may generalize all of these notions to $\mathbb{R}^n$.  An $n$-dimensional crystallographic group is a discrete group of isometries of $\mathbb{R}^n$ which is co-compact.  Fedorov \cite{fed} and Schoenflies \cite{schoen} proved that there are, up to equivalence, 219 crystallographic groups in $\mathbb{R}^3$.  In 1910, Bieberbach proved that for any $n$, there are only finitely many $n$-dimensional crystallographic groups up to equivalence \cite{bieber1, bieber2}, thereby solving Hilbert's 18th problem.

\subsection{Strictly tessellating polytopes and our main result}
Crystal lattices create a perfectly regular pattern, and the fundamental domain of the associated full-rank lattice tessellates space by translation.  Another way to create a perfectly regular pattern is by `strict tessellation.'  This is a notion specific to polytopes.

\begin{definition} \label{def:polytope}
The set of all one dimensional polytopes is the set of all bounded open intervals,
$$\wp_1 := \{ (a,b) : -\infty < a < b < \infty \}.$$
Inductively we define the set of polytopes in $\mathbb{R}^n$ for $n\geq 2$ to be the set of connected bounded domains in $\mathbb{R}^n$ such that $\Omega \in \wp_n$ if and only if
$$\partial \Omega = \bigcup_{j=1} ^m \overline{P_j}, \quad P_j \cong Q_j \in \wp_{n-1}.$$
Above, the boundary of $\Omega$ consists of the closures of $n-1$ dimensional polytopes, $P_j$.  Each $P_j$ is contained in an $n-1$ dimensional hyperplane, which is a set of the form
$$\{ \bx \in \mathbb{R}^n : \textbf{M} \cdot \bx = b \},$$
for some fixed $\textbf{M} \in \mathbb{R}^n$ and $b \in \mathbb{R}$.  Such a hyperplane is canonically identified with $\mathbb{R}^{n-1}$ so that with this identification $P_j$ is canonically identified with $Q_j \in \wp_{n-1}$, which is the meaning of $P_j \cong Q_j$.
\end{definition}

\begin{figure} \centering \includegraphics[width=5cm]{egypt} \includegraphics[width=5cm]{wallpaper2} \caption{Two Egyptian patterns whose symmetry groups are planar crystallographic groups.  These patterns were documented by Owen Jones in 1856 \cite[Egyptian No 7 (plate 10), images 8 and 13]{jones}.}  \label{fig1}
\end{figure}

Next we introduce the notion of strict tessellation. We are not aware of the term `strict tessellation' in the literature, but it may be known under a different name.  An example of a strict tessellation of the plane is given in Figure \ref{fig:tesselating}.  For an example of a tessellation of the plane which is not strict, see Figure \ref{fig:nontesselating}.

\begin{figure}[h]
\centering
        \includegraphics[width=0.48\textwidth]{tiling.pdf}
        \caption{Equilateral triangles are shown here to strictly tessellate the plane.}
        \label{fig:tesselating}
\end{figure}

\begin{definition} \label{defst} A polytope $\Omega \in \wp_n$ strictly tessellates $\mathbb{R}^n$ if
$$\mathbb{R}^n = \bigcup_{j\in \mathbb{Z}} \overline{\Omega_j}.$$
Above, each $\Omega_j$ is isometric to $\Omega$, and they are obtained by reflecting $\Omega$ across its boundary faces.  Furthermore, the hyperplanes which contain the boundary faces of each $\Omega_j$ have empty intersection with (the interior of) $\Omega_k$, for all $j$ and $k$.
\end{definition}

One dimensional polytopes always strictly tessellate because for any real numbers $a<b$,
$$\mathbb{R} = \bigcup_{j \in \mathbb{Z}} \left[ j(b-a)+a, j(b-a) + b \right].$$

\begin{remark}
Any polytope covers $\mathbb{R}^n$ by reflection.  If there were to be a gap, that is a region of $\mathbb{R}^n$ not covered by reflecting the polytope across its boundary faces, such a gap will border a face of a reflected copy of the polytope. Then we simply reflect the polytope again and do this for all bordering faces until we have filled the gap.  Thus one can always cover $\mathbb{R}^n$ by repeatedly reflecting any polytope across its boundary faces, but there may be significant overlap.  If the (interior of the) reflected copies of the polytope have empty intersections, then we obtain a tessellation.
\end{remark}

\begin{figure}[h]
\centering
        \includegraphics[width=0.48\textwidth]{hex.pdf}
        \caption{Although it is well known that regular hexagons tessellate the plane by reflection, the tessellation is not strict, because the lines which contain the edges of the hexagon cut through the interior of the reflected copies.}
        \label{fig:nontesselating}
\end{figure}

In 2008, McCartin proved a remarkable classification theorem \cite{mccartin}, connecting geometry and analysis.  Recall
the Laplacian on $\mathbb{R}^n$ is the partial differential operator
$$\Delta := - \sum_{k=1} ^n \frac{\partial^2}{\partial x_k ^2}.$$
The Laplace eigenvalue problem for a domain $\Omega \subset \mathbb{R}^n$ with the Dirichlet boundary condition is to find all functions $u$ which are not identically zero and satisfy
$$\Delta u = \lambda u, \quad \left . u\right|_{\partial \Omega} = 0, \textrm{ for some constant } \lambda.$$
This is a difficult problem, because in general it is impossible to compute these numbers $\lambda$.  However, using the tools of functional analysis \cite{courhil} one can prove that these eigenvalues are discrete and positive and therefore can be ordered
$$0 < \lambda_1 \leq \lambda_2 \leq \ldots \uparrow \infty.$$
In this way we may speak of the first eigenfunction which has eigenvalue $\lambda_1$.  In one dimension,
by definition \ref{def:polytope}, a polytope is a bounded open interval, $(a,b)$ for some real numbers $a<b$.
The Laplace eigenvalue equation with the Dirichlet boundary condition on such a polytope is to find all functions $u$ such that there exists $\lambda \in \mathbb{C}$ with
$$-u''(x) = \lambda u(x), \quad a < x < b, \quad u(a) = u(b) = 0.$$
Solutions to this equation are
$$u_k (x) = \sin \left( \frac{x-a}{b-a} k \pi \right), \quad \lambda_k = \frac{k^2 \pi^2}{(b-a)^2}, \quad k \in \mathbb{N}.$$
Moreover, Fourier analysis \cite[Chapter 4]{folland} can be used to show that these are \em all \em the solutions to the equation.  These eigenfunctions are all trigonometric functions.  We can also define trigonometric functions on $\mathbb{R}^n$.

\begin{definition} An eigenfunction $u$ for the Laplacian is \em trigonometric \em if it can be expressed as a finite sum of trigonometric functions,
$$u(\bx) = \sum_{j=1} ^m a_j \sin(\textbf{L}_j \cdot \bx) + b_j \cos(\textbf{M}_j \cdot \bx).$$
Above, $a_j, b_j, \in \mathbb{C}$ and $\textbf{L}_j, \textbf{M}_j \in \mathbb{R}^n$ satisfy $||\textbf{L}_j||^2 = ||\textbf{M}_j||^2 = \lambda$ for all $j=1, \ldots, m$.
\end{definition}

\begin{remark}  \label{r:trigef} Since
$$\cos(t) = \sin(t+\pi/2), \quad \forall t \in \mathbb{R},$$
it is equivalent to define a trigonometric eigenfunction to be a function of the form
$$u(\bx) = \sum_{j=1} ^m a_j \sin(\textbf{L}_j \cdot \bx + \phi_j).$$
Above, $a_j \in \mathbb{C}$, $\textbf{L}_j \in \mathbb{R}^n$, $\phi_j \in \mathbb{R}$, and $||\textbf{L}_j||$ are the same for all $j=1, \ldots, m$.
\end{remark}

In general, it is impossible to compute the eigenfunctions of an arbitrary polygonal domain.  Nonetheless, McCartin proved the following classification theorem which shows the equivalence of the analytic property, having trigonometric eigenfunctions, with the geometric property, strictly tessellating.

\begin{theorem}[McCartin \cite{mccartin}] The only polygonal domains in the plane which have a complete set of trigonometric eigenfunctions for the Laplacian with the Dirichlet boundary condition are those which strictly tessellate the plane.  There are precisely four types:  rectangles, isosceles right triangles, equilateral triangles, and hemiequilateral triangles as shown in Figure \ref{fig:strikt_poly}.
\end{theorem}

\begin{remark} \label{r:rectangle} The Laplace eigenfunctions for a rectangular domain with vertices at the points $(0, 0)$, $(a, 0)$, $(0, b)$ and $(a, b)$ with the Dirichlet boundary condition can be computed using separation of variables which reduces the problem to two one-dimensional problems.  The resulting eigenfunctions are indexed by $m, n \in \mathbb{N}$.  For Cartesian coordinates $\bx = (x,y) \in \mathbb{R}^2$, the eigenfunctions are
$$u_{m,n} (x,y) = \sin \left( \frac{ m \pi x}{a} \right) \sin \left( \frac{n \pi y}{b} \right).$$
Using trigonometric identities,
$$u_{m,n} (x,y) = \frac 1 2 \left[ \cos\left( \begin{bmatrix} \frac{m \pi}{a} \\ - \frac{n \pi}{b} \end{bmatrix} \cdot \bx \right)  - \cos\left( \begin{bmatrix} \frac{m \pi}{a} \\ \frac{n \pi}{b} \end{bmatrix} \cdot \bx \right) \right].$$
Consequently, these are trigonometric eigenfunctions.
\end{remark}

\begin{figure}[h]
\includegraphics[width=0.8\textwidth]{strikt_poly.pdf}
\caption{From left to right above we have an equilateral triangle, rectangle, hemi-equilateral triangle, and an isosceles right triangle.  These are the four types of polygons which strictly tessellate the plane.}
\label{fig:strikt_poly}
\end{figure}

Our main result is a generalization to all dimensions.

\begin{theorem} \label{thm1} Assume that $\Omega$ is a polytope in $\mathbb{R}^n$. Then the following are equivalent:
\begin{enumerate}
\item The first eigenfunction for the Laplace eigenvalue equation with the Dirichlet boundary condition extends to a real analytic function on $\mathbb{R}^n$.
\item $\Omega$ strictly tessellates $\mathbb{R}^n$.
\item $\Omega$ is congruent to a fundamental domain of a crystallographic Coxeter group as defined in Bourbaki \cite[VI.25 Proposition 9 p. 180]{bourbaki}, and is also known as an alcove \cite[p. 179]{berard}; see also \S 3.
\end{enumerate}
\end{theorem}

The three equivalent statements in Theorem \ref{thm1} are respectively analytic, geometric, and algebraic.  Our work therefore reveals an intimate connection between analysis, geometry, and algebra.  Moreover, combining our theorem with B\'erard's Proposition \cite[Proposition 9, p. 181]{berard}, we obtain the following rather remarkable result.

\begin{cor}\label{cor0}
Assume that $\Omega$ is a polytope in $\mathbb{R}^n$.  If the first eigenfunction for the Laplace eigenvalue equation with the Dirichlet boundary condition extends to a real analytic function on $\mathbb{R}^n$, then it is a trigonometric eigenfunction.  Moreover, \em all \em the eigenfunctions of $\Omega$ are trigonometric.
\end{cor}

\begin{remark}  Every trigonometric eigenfunction satisfies the first condition of Theorem \ref{thm1}.  However, there are many functions which satisfy this condition but are \em not \em trigonometric.  Examples include the eigenfunctions for a disk in $\mathbb{R}^2$ which are products of Bessel functions and trigonometric functions.  There is no contradiction with the above corollary because a disk is not a polygonal domain.
\end{remark}

\begin{figure}[h]
    \centering
    \scalebox{0.8}{\input{Holy_Trinity3.tex}}
    \caption{This diagram shows the three statements of Theorem \ref{thm1} and how they were proven.}
     \label{fig:Holy_Trinity} \end{figure}

\subsection{Organization}
In \S \ref{s:st}, we prove that if the first eigenfunction of a polytope satisfies the hypotheses of Theorem \ref{thm1}, then the polytope strictly tessellates $\mathbb{R}^n$.  We prove this by generalizing classical results by Lam\'e \cite{lame}.  In \S \ref{s:proof}, we introduce the notions of root systems and alcoves and prove that all polytopes which strictly tessellate $\mathbb{R}^n$ are alcoves.  We then recall the result of B\'erard \cite{berard}:  all alcoves have a complete set of trigonometric eigenfunctions for the Laplace eigenvalue equation with the Dirichlet boundary condition.  These results together complete the proofs of Theorem \ref{thm1} and Corollary \ref{cor0}.
In \S \ref{s:fug} we discuss connections to the Fuglede and Goldbach conjectures.  We make our own conjecture and conclude with a purely geometric conjecture which is equivalent to the strong Goldbach conjecture.

\section{The first eigenfunction and strict tessellation} \label{s:st}
There is no known method to explicitly compute the eigenfunctions and eigenvalues for an arbitrary polytope.  However, using the tools of functional analysis, one can prove general facts about them.  We summarize briefly here.  For the Dirichlet boundary condition for the Laplace eigenvalue equation on a bounded domain, $\Omega \subset \mathbb{R}^n$, the eigenvalues form a discrete positive set which accumulates only at infinity.  We can therefore order the eigenvalues as they increase
$$0<\lambda_1 \leq \lambda_2 \ldots \uparrow \infty.$$
We may correspondingly order the eigenfunctions.  In this way, we may speak of the ``first'' eigenfunction, which is the eigenfunction whose eigenvalue is equal to $\lambda_1$.  The eigenfunctions form an orthogonal basis of the Hilbert space $\mathcal{L}^2(\Omega)$.  We shall require the following well-known fact about the first eigenfunction.  The proof of this theorem can be found in the classical PDE textbook of Evans \cite[\S 6.5]{evans}.

\begin{theorem} \label{t:lam1} Let $\Omega$ be a bounded domain in $\mathbb{R}^n$.  Then the first eigenfunction for the Laplace eigenvalue equation with the Dirichlet boundary condition does not vanish in $\Omega$.
\end{theorem}

The following result is originally due to Lam\'e \cite{lame} in two dimensions and specified to trigonometric eigenfunctions.  Here we generalize the result to $\mathbb{R}^n$ for all $n$ as well as to all real analytic functions.  We include a short proof for completeness.

\begin{lemma}[Vanishing planes]  Let $u$ be a real analytic function on $\mathbb{R}^n$.  Assume that $u$ vanishes on an open, non-empty, subset of a hyperplane,
$${\mathcal{P}} := \{ \bx \in \mathbb{R}^n : \textbf{M} \cdot \bx = b \}.$$
Then $u$ vanishes on all of ${\mathcal{P}}$.
\end{lemma}
\begin{proof}
Let $\textbf{y}$ be any point of ${\mathcal{P}}$.  Let $\bx$ be a point in the interior of the open neighborhood on which $u$ vanishes.  We parametrize a line segment to join the points $\bx$ and $\textbf{y}$, defining
$$l(t) := t\bx + (1-t) \textbf{y}, \quad 0 \leq t \leq 1.$$
We note that
$$\textbf{M} \cdot l(t) = t \textbf{M} \cdot \bx + (1-t) \textbf{M} \cdot \textbf{y} \equiv b, \quad \forall t \in [0,1].$$
Consequently this line segment is contained entirely in ${\mathcal{P}}$.  Now we consider the function
$$u(t) := u(l(t)).$$
Since $u$ is a real analytic function, $u(t)$ is also a real analytic function of $t$.  Since $\bx$ is in the interior of the open neighborhood on which $u$ vanishes, and this neighborhood is contained in ${\mathcal{P}}$, there is $\varepsilon > 0$ such that $l(t)$ is contained in this neighborhood for all $t \in [1-\varepsilon, 1]$.  The function $u(t)$ vanishes identically on $[1-\varepsilon, 1]$.  Consequently it vanishes identically for all $t$, and in particular
$$u(l(0)) = u(\textbf{y}) = 0.$$
Since $\textbf{y} \in {\mathcal{P}}$ was arbitrary, we obtain that $u$ vanishes at every point of ${\mathcal{P}}$.
\end{proof}

We will also generalize Lam\'e's Fundamental Theorem, which was originally proven in two dimensions and for trigonometric functions, to $n$ dimensions and real analytic eigenfunctions.

\begin{theorem}[Lam\'e's Fundamental Theorem] Assume that $u$ is a real analytic function on $\mathbb{R}^n$ which satisfies the Laplace eigenvalue equation with the Dirichlet boundary condition on a polytope $\Omega \in \wp_n$.  Then $u$ is anti-symmetric with respect to all $(n-1)$ dimensional hyperplanes on which $u$ vanishes.
\end{theorem}

\begin{proof}
Let $\lambda$ be the eigenvalue corresponding to $u$, so that on $\Omega$ we have
$$\Delta u(\bx) = \lambda u(\bx) \quad \forall \bx \in \Omega.$$
Then, since $u$ is real analytic, $\Delta u$ is also real analytic on $\Omega$.  The function
$$\Delta u - \lambda u$$
is real analytic and vanishes on $\Omega$ which is an open subset of $\mathbb{R}^n$.  Consequently, it vanishes on all of $\mathbb{R}^n$, and therefore $u$ satisfies the same Laplace eigenvalue equation on all of $\mathbb{R}^n$.

Now, let $H$ be an $(n-1)$ dimensional hyperplane on which $u$ vanishes.  We consider coordinates
$$\textbf{y} = (r, {\bf z}) \in \mathbb{R}^n,$$
where ${\bf z} \in \mathbb{R}^{n-1}$ and
$$r=0 \iff \textbf{y} \in H.$$
Let us now define the function
$$\widetilde u(r, {\bf z}) := \begin{cases} u(r, {\bf z}) & (r, {\bf z}) \in \overline \Omega \\ - u(-r, {\bf z}) & (r, {\bf z}) \not \in \overline \Omega. \end{cases}$$
With this definition, $\widetilde u$ is anti-symmetric with respect to $H$.  Since $u$ satisfies the Laplace eigenvalue equation on $\mathbb{R}^n$ and is real analytic, the same is true for $\widetilde u$.  Moreover, on $\partial \Omega$, $u$ and $\widetilde u$ have the same Cauchy data; both functions vanish and have the same normal derivative.  Consequently, by standard uniqueness theory of partial differential equations \cite{courhil}, \cite{evans}, $u = \widetilde u$.  Therefore $u$ is also anti-symmetric with respect to $H$.
\end{proof}

We are now poised to prove the first implication in Theorem \ref{thm1}.

\begin{prop} \label{prop1thm1} Assume that $\Omega$ is a polytope in $\mathbb{R}^n$, and the first eigenfunction satisfies the first condition of Theorem \ref{thm1}.  Then $\Omega$ strictly tessellates $\mathbb{R}^n$.
\end{prop}

\begin{proof}
Let $\Omega$ be a polytope in $\mathbb{R}^n$ as in the statement of the proposition.  If $n=1$, then $\Omega$ is a segment and may be written as $(a,b)$ for some real numbers $a<b$.  We have computed the eigenfunctions explicitly in this case.  They are
$$u_k (x) = \sin \left( \frac{x-a}{b-a} k \pi \right).$$
The first eigenfunction in particular satisfies the hypotheses of Theorem \ref{thm1}, and we have also shown that all one dimensional polytopes strictly tessellate $\mathbb{R}^1$.  Hence the proposition is proven in one dimension.  So let us assume that $n \geq 2$.  By the Vanishing Planes Lemma, for an affine subset ${\mathcal{P}}$ which contains a boundary face of $\Omega$, all eigenfunctions of $\Omega$ vanish on ${\mathcal{P}}$.  Since the first eigenfunction never vanishes in the interior of $\Omega$ by Theorem \ref{t:lam1}, it follows that all of the hyperplanes which contain the boundary faces of $\Omega$ have empty intersection with the interior of $\Omega$.

Near a boundary face of $\Omega$ which is contained in the hyperplane ${\mathcal{P}}$, we use the coordinates $(r, {\bf z})$ as in the previous Theorem.   Without loss of generality, we may assume that for $r>0$ sufficiently small, there are points $(r,{\bf z})$ inside $\Omega$.  Since $\Omega$ is a domain, it is an open set.  By Lam\'e's Fundamental Theorem, $u_1(-r, {\bf z}) = -u_1(r,{\bf z})$.  Consequently, reflecting $\Omega$ across ${\mathcal{P}}$, we see that $u_1$ is also an eigenfunction in this reflected copy of $\Omega$ which satisfies the Dirichlet boundary condition.  Since the first eigenfunction does not vanish on the interior of $\Omega$, it also does not vanish on the interior of the reflection of $\Omega$ across ${\mathcal{P}}$.  This shows that every hyperplane which contains a boundary face of $\Omega$ has empty intersection with this reflected copy of $\Omega$.  We repeat this argument, reflecting across all of the boundary faces of $\Omega$ and then reflecting across the boundary faces of the reflected copies of $\Omega$. This shows that the hyperplanes which contain the boundary faces of the reflected copies of $\Omega$ must always have empty intersection with the interior of all reflected copies of $\Omega$.  Consequently, repeating indefinitely, we obtain a strict tessellation of $\mathbb{R}^n$ by $\Omega$.
\end{proof}

\section{Root systems, alcoves, and strictly tessellating polytopes} \label{s:proof}
In 1980, Pierre B\'erard showed that a certain type of bounded domain in $\mathbb{R}^n$, known as an \em alcove, \em always has a complete set of trigonometric eigenfunctions for the Laplace eigenvalue equation with the Dirichlet boundary condition.  To define alcoves, we must first define root systems.  The concept of a root system was originally introduced by Wilhelm Killing in 1888 \cite{killing}.  His motivation was to classify all simple Lie algebras over the field of complex numbers.  In this section, we will see how our analytic problem, the study of the Laplace eigenvalue equation, is connected to these abstract algebraic concepts from Lie theory and representation theory.

\begin{definition}\label{def:rootsystem}
A \em root system \em in $\mathbb{R}^n$ is a finite set $R$ of vectors which satisfy:
\begin{enumerate}
\item $0$ is not in $R$.
\item The vectors in $R$ span $\mathbb{R}^n$.
\item The only scalar multiples of $\textbf{v} \in R$ are $\pm \textbf{v}$.
\item $R$ is closed with respect to reflection across any hyperplane whose normal is an element of $R$, that is
$$\textbf{v} - 2 \frac{ \textbf{u} \cdot \textbf{v}}{||\textbf{u}||^2} \textbf{u} \in R, \quad \forall \textbf{u}, \textbf{v} \in R;$$
\item If $\textbf{u}, \textbf{v} \in R$ then the projection of $\textbf{u}$ onto the line through $\textbf{v}$ is an integer or half-integer multiple of $\textbf{v}$.  The mathematical formulation of this is that
$$2 \frac{ \textbf{u} \cdot \textbf{v}}{||\textbf{v}||^2} \in \mathbb{Z}, \quad \forall \textbf{u}, \textbf{v} \in R.$$
\end{enumerate}
The elements of a root system are often referred to as \em roots.  \em
\end{definition}

\begin{remark}  There are different variations of the definition given above for a root system depending on the context.  Sometimes only conditions 1--4 are used to define a root system.  When the additional assumption (5) is included then the root system is said to be \em crystallographic.  \em  In other contexts, condition 3 is omitted, and they would call a root system which satisfies condition 3 \em reduced.  \em
\end{remark}

We will need the dual root system to define the eigenvalues of the polytope which will be naturally associated to the root system.

\begin{definition} \label{def:rootinverse}
Let $R$ be a root system.  Then for $\textbf{v} \in R$ the \em coroot \em $\textbf{v}^\vee$ is defined to be
$$\textbf{v}^\vee = \frac{2}{||\textbf{v}||^2} \textbf{v}.$$
The set of coroots $R^\vee := \{ \textbf{v}^\vee \}_{\textbf{v} \in R}$.  This is called the \em dual root system, \em and may also be called the inverse root system.  The dual root system is itself a root system.
\end{definition}

We associate a Weyl group to a root system.  These Weyl groups are subgroups of the orthogonal group $O(n)$.

\begin{definition} \label{def:weylgroup}
For any root system $R\subset \mathbb{R}^n$ we associate a subgroup of the orthogonal group $O(n)$ known as its \em Weyl group. \em  This is the subgroup $W < O(n)$ generated by the set of reflections by hyperplanes whose normal vectors are elements of $R$.   For $\textbf{v} \in R$ reflection across the hyperplane with normal vector equal to $\textbf{v}$ is explicitly
$$\sigma_\textbf{v} : \mathbb{R}^n \to \mathbb{R}^n, \quad \sigma_\textbf{v} (\bx) = \bx - 2 \frac{(\textbf{v}\cdot \bx)}{||\textbf{v}||^2} \textbf{v}.$$
\end{definition}

\begin{figure}[h]
    \centering
    \scalebox{0.8}{\input{Rotsystem_2D.tex}}
    \caption{Here are four root systems in $\mathbb{R}^2$.  Below each root system is the name of its Weyl group. The name of the Weyl group may also be used as the name of the root system.}
    \label{fig:rotsystem}
\end{figure}

By the definition of a root system, the associated Weyl group is finite.  To explain what was proven in \cite{berard} by B\'erard, we require the notion of Weyl chamber.

\begin{definition} \label{def:chamber}
For a root system $R \subset \mathbb{R}^n$ for each $\textbf{v} \in R$, let $H_\textbf{v}$ denote the hyperplane which contains the origin and whose normal vector is $\textbf{v}$.  In particular
$$H_\textbf{v} := \{ \bx \in \mathbb{R}^n : \bx \cdot \textbf{v} = 0 \}.$$
Let $\mathcal H = \{ H_\textbf{v} \}_{\textbf{v} \in R}$.  Then $\mathbb{R}^n \setminus \mathcal H$ is disconnected, and each connected open component is known as a \em Weyl chamber. \em
\end{definition}

\begin{figure}[h]
    \centering
    \scalebox{0.8}{\input{Weyl_chamb.tex}}
    \caption{Extending the shaded area to infinity it is a Weyl chamber of the Weyl group $A_2$.}
    \label{fig:Weyl_chamb}
\end{figure}

To obtain an alcove, we not only wish to consider reflecting across the hyperplanes $H_v$, but also across a discrete set of parallel translations of these hyperplanes.

\begin{figure}[h]
   \centering
    \scalebox{0.8}{\input{alcoveB2.tex}}
    \caption{This shows an alcove, $A$, corresponding to the root system with Weyl group $B_2$. For $\alpha \in B_2$, the hyperplanes $H_{\alpha,k}$ for $k \in \mathbb{Z}$ are the parallel hyperplanes which have normal vector equal to $\alpha$. Note that $A$ is an isosceles right triangle.}
    \label{fig:alcove}
    \end{figure}

\begin{definition} \label{def:alcove} Let $R$ be a root system.  Denote by $H_\textbf{v}$ the hyperplane in $\mathbb{R}^n$ which contains the origin and whose normal vector is equal to $\textbf{v}$ for $\textbf{v} \in R$.  Let
$$H_{\textbf{v},k} = \{ \bx \in \mathbb{R}^n : \textbf{v} \cdot \bx = k\},$$
for $k \in \mathbb{Z}$.  Then $H_{\textbf{v},0} = H_\textbf{v}$.  For $k \neq 0$, the hyperplane $H_{\textbf{v},k}$ is parallel to $H_\textbf{v}$.  We define an \em alcove \em to be a connected component of
$$\mathbb{R}^n \setminus \left\{ \bigcup_{\textbf{v} \in R, k \in \mathbb{Z}} H_{\textbf{v},k} \right\}.$$
\end{definition}

We note that the definition of an alcove immediately implies that it is a polytope in $\mathbb{R}^n$. An example of an alcove is shown in Figure \ref{fig:alcove}.

\begin{prop}[Proposition 9, p. 181 \cite{berard}] \label{prop1} Let $\Omega \subset \mathbb{R}^n$ be an alcove.  Then $\Omega$ has a complete set of trigonometric eigenfunctions for the Laplace eigenvalue equation with the Dirichlet boundary condition.
\end{prop}

For readers who understand French and read \cite{berard}, you may notice that the statement of Proposition \ref{prop1} is not the English translation of  \cite[Proposition 9, p. 181]{berard}.  B\'erard proved a stronger result; he specified the eigenvalues and corresponding eigenfunctions.  To understand what B\'erard proved, let $R$ be a root system.  Let $C(R)$ denote a Weyl chamber, and let $D(R)$ denote an alcove which is contained in the Weyl chamber $C(R)$.  Consider the dual root system $R^\vee$.  The vertices of the closures of the alcoves associated to $R^\vee$ create a lattice.  Let us denote this lattice by $\Gamma$.  The dual lattice
$$\Gamma^* := \{ \bx \in \mathbb{R}^n : \bx \cdot \boldsymbol\gamma \in \mathbb{Z}, \quad \forall \boldsymbol\gamma \in \Gamma\}.$$
B\'erard referred to the points contained in this dual lattice as `the group of weights of $R$' (le groupe des poids de R) \cite{berard}.  He proved that the eigenvalues for the alcove $D(R)$ are given by
$$\{ 4 \pi^2 ||\textbf{q}||^2 : \textbf{q} \in \Gamma^* \cap C(R)\}.$$
The multiplicity of the eigenvalue $\lambda = 4 \pi^2 ||\textbf{q}||^2$ is equal to the number of vectors $\textbf{q} \in \Gamma^* \cap C(R)$ which satisfy $\lambda = 4 \pi^2 ||\textbf{q}||^2$.  The eigenfunctions are certain linear combinations of $e^{2\pi i \bx \cdot w(\textbf{q})}$, where $w(\textbf{q})$ is in the affine Weyl group of $R$.  The affine Weyl group of $R$ is the semi-direct product of the Weyl group and the lattice $\Gamma$.  Combining our Proposition \ref{prop1thm1} with B\'erard's Proposition \ref{prop1} we obtain the following corollary which states that every alcove is a strictly tessellating polytope.

\begin{cor} \label{cor1} Let $\Omega \subset \mathbb{R}^n$ be an alcove.  Then $\Omega$ is a polytope which strictly tessellates $\mathbb{R}^n$.
\end{cor}

In the following proposition we prove the converse: every strictly tessellating polytope is an alcove of a root system.

\begin{prop} \label{prop2thm1} Let $\Omega \subset \mathbb{R}^n$ be a polytope which strictly tessellates $\mathbb{R}^n$.  Then $\Omega$ is an alcove.
\end{prop}

\begin{proof}
We will build a root system, $R$, using the fact that $\Omega$ strictly tessellates space.   The main idea is to collect all the normal vectors to all faces of $\Omega$ in the tessellation.  Since we may translate the entire tessellation, we will collect these normal directions by translating each vertex to the origin.  We will specify the lengths of the vectors later, using the fact that the tessellation is strict.

So, begin with a vertex at the origin.  For each of the boundary faces which contains this vertex at the origin, reflection across the boundary face is an element of the orthogonal group $O(n)$.  Since the tessellation is strict, and reflection fixes the origin, after finitely many reflections across the boundary face and its images under reflection, we must return to the original copy of $\Omega$.  We begin to build the root system $R$ by defining it to contain $\pm \textbf{v}$, the normal vectors to this boundary face, and to also contain all the images of $\pm \textbf{v}$ under the repeated reflection until we return to the original copy of $\Omega$.  The lengths of the vectors will be specified later.  Since this repeated reflection returns to the original copy of $\Omega$ after finitely many reflections, the set of vectors defined in this way is finite.  Repeat this for all the boundary faces which contain the origin, and include any new normal directions $\pm \textbf{u}$ in $R$.  Next, consider a different vertex in $\Omega$.  Translate the entire tessellation so that this vertex is at the origin and repeat, including any new normal directions $\pm \textbf{w}$ in the set $R$.  Do this for all vertices.  Each time we include at most finitely many vectors.  Since the number of vertices is also finite, the total resulting set of vectors, denoted by $R$ is a finite set of vectors.  We define this set so that if $\textbf{v} \in R$ then $-\textbf{v} \in R$ as well.  This specifies the directions of all vectors, but not their lengths.  The polytope is a bounded, connected, open set with boundary consisting of flat faces.  The collection of normal vectors to the faces of $\Omega$ therefore span $\mathbb{R}^n$, for if not, $\Omega$ would be contained in a $k$-dimensional hyperplane in $\mathbb{R}^n$ and thus not an open set in $\mathbb{R}^n$.  Since the set $R$ contains all the normal vectors to all boundary faces of $\Omega$, the set $R$ spans $\mathbb{R}^n$.  As we have defined it, $0 \not \in R$.

Now fix a tessellation by $\Omega$.  Without loss of generality at least one vertex in $\overline \Omega$ is at the origin.  The tessellation defines hyperplanes in $\mathbb{R}^n$ which contain the boundary faces of the copies of $\Omega$ in the tessellation.  The way we have defined the set $R$, it contains all the normal directions of all these hyperplanes.  Since $R$ is a finite set, and there are countably many hyperplanes defined by the tessellation, this means that for each $\textbf{v} \in R$, we may enumerate the hyperplanes whose normal direction is $\pm \textbf{v}$ by $H_{\textbf{v}, k}$ for $k \in \mathbb{Z}$.  These hyperplanes are parallel.  Consider those $\textbf{v} \in R$ such that there is a hyperplane, denoted by $H_{\textbf{v}, 0}$, whose normal vector is $\pm \textbf{v}$ and which contains the origin.  We will define the length of $\textbf{v}$ so that the closest parallel hyperplanes to $H_{\textbf{v}, 0}$ are
$$H_{\textbf{v}, \pm 1} = \{ \bx \in \mathbb{R}^n : \bx \cdot \textbf{v} = \pm 1\}.$$

Since the hyperplanes $H_{\textbf{v}, 1}$ and $H_{\textbf{v},0}$ are parallel, there is a point $\textbf{y}_\textbf{v} \in \mathbb{R}^n$ with $\textbf{y}_\textbf{v}=c_\textbf{v} \textbf{v}$ for some $c_\textbf{v} \in \mathbb{R}$, where $\textbf{v}$ is the normal vector of (currently) unknown length.  This point $\textbf{y}_\textbf{v}$ is the unique point in $H_{\textbf{v}, 1}$ which is closest to the origin.  Since $\textbf{y}_\textbf{v} \in H_{\textbf{v},1}$ we have by definition
$$\textbf{y}_\textbf{v} \cdot \textbf{v} = c_\textbf{v} ||\textbf{v}||^2 = 1.$$
We know the location of $H_{\textbf{v}, 1}$ since it is given by the tessellation, so we also know the length of $\textbf{y}_\textbf{v}$.
Using the above equation,
$$||\textbf{y}_\textbf{v}|| = |c_\textbf{v}| ||\textbf{v}||, \quad \textbf{y}_\textbf{v} \cdot \textbf{v} = c_\textbf{v} ||\textbf{v}||^2 = 1 \implies c_\textbf{v} = \frac{1}{||\textbf{v}||^2} \implies ||\textbf{y}_\textbf{v}|| = \frac{1}{||\textbf{v}||},$$
so we conclude that
$$||\textbf{v}|| = \frac{1}{||\textbf{y}_\textbf{v}||}.$$
Geometrically,
$$||\textbf{y}_\textbf{v}|| = \textrm{ the distance between the hyperplanes $H_{\textbf{v},0}$ and $H_{\textbf{v}, \pm1}$}.$$
We also have
$$||\textbf{v}|| = \frac{1}{||\textbf{y}_\textbf{v}||}, \quad \textbf{y}_\textbf{v} = \frac{\textbf{v}}{||\textbf{v}||^2}.$$
A schematic image of this is given in Figure \ref{fig:prop_pic}.
\begin{figure}[h]
   \centering
    \scalebox{0.8}{\input{prop_pic.tex}}
    \caption{Given a polytope $\Omega$, we construct the hyperplanes $H_{\textbf{v},0}$, here in black dotted lines, and the normal vectors $\textbf{v}$. The set $\{H_{\textbf{v},k}\}$ includes the gray dotted lines.}
    \label{fig:prop_pic}
    \end{figure}

To define the lengths of all the elements of $R$, take each vertex in a single copy of $\Omega$, and translate the entire tessellation so that the vertex is at the origin.  Repeat the above argument to define the lengths of all the vectors in $R$.

For $\textbf{v} \in R$, reflection across a hyperplane in the tessellation whose normal direction is that of $\textbf{v}$ may only permute the ensemble of hyperplanes by virtue of the strict tessellation.  Let $\textbf{w} \in R$.  By possibly translating the entire picture, assume that there is a hyperplane in the tessellation with normal direction $\pm \textbf{w}$ and which contains the origin, such that the origin is a vertex of a copy of $\Omega$ in the tessellation.  Thus $H_{\textbf{w}, 0}$ is a hyperplane in the tesssellation.  Consider the reflection with normal direction $\textbf{v}$, denoted by $\sigma_\textbf{v}$, that is
$$\sigma_\textbf{v} (\bx) = \bx - 2 \frac{ \bx \cdot \textbf{v}}{||\textbf{v}||^2} \textbf{v}.$$
Then $\sigma_\textbf{v} (0) = 0$.  Consequently, $\sigma_\textbf{v} (H_{\textbf{w}, 0})$ is another hyperplane in the strict tessellation which also contains the origin, thus it is $H_{\textbf{u}, 0}$ for some $\textbf{u} \in R$.  Similarly, we also have $\sigma_\textbf{v} (H_{\textbf{w}, 1}) = H_{\textbf{u}, j}$ for some $j \in \mathbb{Z}$.  Since $\sigma_\textbf{v}$ preserves the scalar product, for $\bx \in H_{\textbf{w}, 1}$, by definition
$$\bx \cdot \textbf{w} = 1 \implies \sigma_\textbf{v} (\bx) \cdot \sigma_\textbf{v} (\textbf{w}) = 1.$$
Since $\sigma_\textbf{v}$ sends $\bx$ to a point in $H_{\textbf{u}, j}$ we also have
$$\sigma_\textbf{v} (\bx) \cdot \textbf{u} = j.$$
Since $\sigma_\textbf{v} (H_{\textbf{w},0}) = H_{\textbf{u},0}$, we must have that $\sigma_\textbf{v} (\textbf{w}) = \alpha \textbf{u}$ for some $\alpha \in \mathbb{R}$.  Therefore combining with the above we obtain
$$1 = \bx \cdot \textbf{w} = \sigma_\textbf{v} (\bx) \cdot \sigma_\textbf{v} (\textbf{w}) = \alpha \sigma_\textbf{v} (\bx) \cdot \textbf{u} = \alpha j \implies \alpha = \frac 1 j.$$
So we have proven that
$$\sigma_\textbf{v} (\textbf{w}) = \frac 1 j \textbf{u}.$$

Now let us consider what happens to $\textbf{y}_\textbf{w}$ when we reflect by $\sigma_\textbf{v}$.  The vector $\textbf{y}_\textbf{w}$ is orthogonal to the hyperplanes $H_{\textbf{w},0}$ and $H_{\textbf{w}, 1}$.  The vector goes from the origin in $H_{\textbf{w},0}$ and since its length is the distance between $H_{\textbf{w},0}$ and $H_{\textbf{w}, 1}$, the end of this vector sits in $H_{\textbf{w}, 1}$.  When this vector is reflected by $\sigma_\textbf{v}$, it will again start from the origin and have its endpoint lying on one of the parallel hyperplanes, by virtue of the strict tessellation.  We compute explicitly that
$$\sigma_\textbf{v} (\textbf{y}_\textbf{w}) = \textbf{y}_\textbf{w} - 2 \textbf{y}_\textbf{w} \cdot \textbf{v} \frac{\textbf{v}}{||\textbf{v}||^2} = \textbf{y}_\textbf{w} - 2 (\textbf{y}_\textbf{w} \cdot \textbf{v}) \textbf{y}_\textbf{v}.$$
On the other hand, since $\sigma_\textbf{v} (\textbf{w}) = \frac 1 j \textbf{u}$, we compute using the definitions of $\textbf{y}_\textbf{w}$, $\textbf{y}_\textbf{v}$, and $\textbf{y}_\textbf{u}$ that
$$\sigma_\textbf{v} (\textbf{y}_\textbf{w}) = \sigma_\textbf{v} \left( \frac{\textbf{w}}{||\textbf{w}||^2} \right) = \frac{1}{||\textbf{w}||^2} \sigma_\textbf{v} (\textbf{w}) = \frac{1}{||\textbf{w}||^2} \frac 1 j \textbf{u}.$$
Now, since
$$||\textbf{u}||^2 = j^2 ||\textbf{w}||^2 \implies \sigma_\textbf{v} (\textbf{y}_\textbf{w}) = j \left( \frac{\textbf{u}}{||\textbf{w}||^2 j^2} \right) = j \textbf{y}_\textbf{u}.$$
Combining these calculations we obtain
$$\sigma_\textbf{v} (\textbf{y}_\textbf{w}) = \textbf{y}_\textbf{w} - 2 (\textbf{y}_\textbf{w} \cdot \textbf{v}) \textbf{y}_\textbf{v} = j \textbf{y}_\textbf{u} \implies 2 (\textbf{y}_\textbf{w} \cdot \textbf{v}) \textbf{y}_\textbf{v} = \textbf{y}_\textbf{w} - j \textbf{y}_\textbf{u}.$$
The vector $\textbf{y}_\textbf{w}$ goes from the origin to $H_{\textbf{w}, 1}$, while the vector $-j \textbf{y}_\textbf{u}$ goes from the origin to $H_{\textbf{u}, -j}$.  By vector addition and the strict tessellation, the sum $\textbf{y}_\textbf{w} - j \textbf{y}_\textbf{u}$ must go from the origin and end precisely at one of the parallel hyperplanes.  Consequently, the vector
$$2 (\textbf{y}_\textbf{w} \cdot \textbf{v}) \textbf{y}_\textbf{v}$$
must be an integer multiple of $\textbf{y}_\textbf{v}$ because it goes from the origin in the direction of $\textbf{y}_\textbf{v}$ and lands at one of the parallel hyperplanes $H_{\textbf{v}, k}$ for some $k \in \mathbb{Z}$.  Therefore,
$$2 (\textbf{y}_\textbf{w} \cdot \textbf{v}) = k \in \mathbb{Z}.$$
By the definitions of $\textbf{y}_\textbf{w}$ and $\textbf{v}$,
$$2 (\textbf{y}_\textbf{w} \cdot \textbf{v}) = 2 \frac{\textbf{w} \cdot \textbf{v}}{||\textbf{w}||^2} = k \in \mathbb{Z}.$$
In a similar way, reversing the roles of $\textbf{w}$ and $\textbf{v}$, we also obtain
$$2 \frac{\textbf{v} \cdot \textbf{w}}{||\textbf{v}||^2} \in \mathbb{Z}.$$
Since $\textbf{w}, \textbf{v} \in R_\Omega$ were arbitrary, this shows the final condition needed for $R_\Omega$ to be a root system in Definition \ref{def:rootsystem} is satisfied.  We conclude that $R_\Omega$ is a root system and that $\Omega$ is one of its alcoves.
\end{proof}

The proofs of Theorem \ref{thm1} and Corollary \ref{cor0} will now follow from Propositions \ref{prop1thm1} and \ref{prop2thm1} and B\'erard's Proposition \ref{prop1}.

\subsubsection{Proof of Theorem \ref{thm1}}
In Proposition \ref{prop1thm1} we proved that if $\Omega$ is a polytope, then the statements in Theorem \ref{thm1} satisfy:  $1 \implies 2$.  In Proposition \ref{prop2thm1}, we proved that the statements in Theorem \ref{thm1} satisfy:  $2 \implies 3$.  Finally, by B\'erard's Proposition \ref{prop1}, statement $3$, that $\Omega$ is an alcove, implies that all the eigenfunctions of $\Omega$ are trigonometric.  All trigonometric eigenfunctions are real analytic on $\mathbb{R}^n$ and satisfy the Laplace eigenvalue equation on $\mathbb{R}^n$.  Consequently, $3 \implies 1$.  The three statements are therefore equivalent.
\qed

\subsubsection{Proof of Corollary \ref{cor0}}
If the first eigenfunction of a polytope in $\mathbb{R}^n$ satisfies the hypotheses of Theorem \ref{thm1}, then the polytope is an alcove.  By B\'erard's Proposition \ref{prop1}, all of the eigenfunctions of the polytope are trigonometric.
\qed

\section{Concluding remarks and conjectures}  \label{s:fug}
We have now answered the analysis question:  when does a polytope in $\mathbb{R}^n$ have a complete set of trigonometric eigenfunctions for the Laplace eigenvalue equation?  In geometric terms, the necessary and sufficient condition for a polytope to have a complete set of trigonometric eigenfunctions is that the polytope strictly tessellates $\mathbb{R}^n$.  In algebraic terms, in the language of Bourbaki,  the equivalent necessary and sufficient condition is that the polytope is congruent to a fundamental domain of a crystallographic Coxeter group \cite[VI.25 Proposition 9 p. 180]{bourbaki}, \cite[p. 179]{berard}.  Returning to the analysis problem, it is interesting to note that it is enough to know that the first eigenfunction is real analytic and satisfies the Laplace eigenvalue equation on $\mathbb{R}^n$ to conclude that it is a trigonometric function and moreover, \em all \em the eigenfunctions are trigonometric.  This is a remarkable fact.  Moreover, the equivalence of analytic, geometric, and algebraic statements shows that these different areas of mathematics are intimately connected.  The Fuglede conjecture similarly brings together different areas of mathematics in the study of a single question.

\subsection{The Fuglede Conjecture}
To state the Fuglede Conjecture, we introduce a few concepts.
\begin{definition}
A domain $\Omega \subset \mathbb{R}^d$ is said to be a \em spectral set \em if there exists $\Lambda \subset \mathbb{R}^n$ such that the functions
$$\{ e^{2\pi i \lambda \cdot \bx} \}_{\lambda \in \Lambda}$$
are an orthogonal basis for $\mathcal{L}^2 (\Omega)$.  The set $\Lambda$ is then said to be \em a spectrum of $\Omega$, \em  and $(\Omega, \Lambda)$ is called \em a spectral pair.  \em
\end{definition}

To relate these notions to our work here, we observe that if a domain $\Omega$ were to have all its eigenfunctions for the Laplace eigenvalue equation of the form $e^{2\pi i \lambda \cdot \bx}$, then these functions would comprise an orthogonal basis for $\mathcal{L}^2 (\Omega)$. Consequently, knowing that the eigenfunctions are precisely of this form implies that the domain is a spectral set.  However, the converse is not true, in the sense that if $\Omega$ is a spectral set, then its eigenfunctions are not necessarily individual complex exponential functions.  If $\Omega$ is a spectral set, then the eigenfunctions must be linear combinations of the $e^{2\pi i \lambda \cdot \bx}$, since these are a basis for $\mathcal{L}^2 (\Omega)$.  However, the linear combinations could have countably infinitely many terms, so it is not clear what precise form the eigenfunctions will take.

\begin{conj}[Fuglede \cite{fug}] Every domain of $\mathbb{R}^n$ which has positive Lebesgue measure is a spectral set if and only if it tiles $\mathbb{R}^n$ by translation.
\end{conj}

Fuglede proved in 1974 that the conjecture holds if one assumes that the domain is the fundamental domain of a lattice \cite{fug}.  Only several years later, in 2003, was further progress made by Iosevich, Katz, and Tao \cite{ios} who proved that the Fuglede conjecture is true if one restricts to convex planar domains.  In the following year, Tao proved that the Fuglede Conjecture is false in dimension 5 and higher \cite{tao}.  In 2006, the works of Farkas, Kolounzakis, Matolcsi and Mora \cite{fkmm}, \cite{km}, \cite{km2}, \cite{mat} proved that the conjecture is also false for dimensions 3 and 4.  In 2017, Greenfeld and Lev proved that Fuglede's Conjecture is true if one restricts attention to domains which are convex polytopes, but only in $\mathbb{R}^3$ \cite{green}.  In 2019, Lev and Matolcsi proved that Fuglede's Conjecture is true if one restricts attention to convex domains, in any dimension \cite{lm}.  Interestingly, the Fuglede Conjecture is still an open problem for arbitrary domains in dimensions one and two.  Here we make the following conjecture which is related to yet independent from Fuglede's.

\begin{conj} Let $\Omega$ be a domain in $\mathbb{R}^n$.  Then $\Omega$ has a complete set of trigonometric eigenfunctions for the Laplace eigenvalue equation with the Dirichlet boundary condition if and only if $\Omega$ is a polytope which strictly tessellates $\mathbb{R}^n$.  Equivalently, $\Omega$ has a complete set of trigonometric eigenfunctions for the Laplace eigenvalue equation with the Dirichlet boundary condition if and only if $\Omega$ is an alcove.
\end{conj}

\begin{figure}[h] \centering
\includegraphics[width=8cm]{zero.pdf}
\caption{This figure shows the null set of the function
$u(x,y) = \sin(x)+\sin(y)+\sin((x+y)/\sqrt{2})$ in a square shaped region of $\mathbb{R}^2$.  The null set includes the line $y=-x$ as well as the other curves in the region.  Consequently,  by uniqueness, this function is the first eigenfunction of the connected, open domains which are bounded by these curves, since it vanishes on the boundary but not on the interior and satisfies the Laplace eigenvalue equation.  Hence the first eigenfunction satisfies the first condition of Theorem \ref{thm1}, but we do not obtain any further conclusions because the domain is not a polytope.}
\label{fig:nullset}
\end{figure}

The difficulty in treating arbitrary domains is that we do not have a replacement for Lam\'e's results which are central to our proof.  Moreover, it is possible to construct linear combinations of trigonometric functions which vanish on curved regions; an example is given in Figure \ref{fig:nullset}. Consequently, we cannot immediately conclude that domains which have trigonometric eigenfunctions have flat boundary faces, and hence they are polytopes.  A domain with a curved boundary could have a few trigonometric eigenfunctions.  What is reasonable to expect, however, is that it does not have a \em complete \em set of trigonometric eigenfunctions.
\subsection{The crystallographic restriction theorem and a geometric approach to the Goldbach Conjecture}
The vertices of the strict tessellation given by a polytope which is an alcove are in fact the set of points in a full rank lattice.  We note that two different polytopes may give rise to the same lattice; for example an isosceles right triangle and the square obtained by two copies of that triangle will produce the same lattice.   For any discrete group of isometries of $\mathbb{R}^n$, an element $g$ of such a group has finite order if there is an integer $k>0$ such that $g$ composed with itself $k$ times is the identity.  The minimal such $k$ is the order of $g$.  To state the crystallographic restriction theorem, we define a function which is like an extension of the Euler totient function.  For an odd prime $p$ and $r\geq 1$,
$$\psi(p^r) := \phi(p^r), \quad \phi(p^r) = p^r - p^{r-1}.$$
Above, $\phi$ is the Euler totient function.  The Euler totient function of a positive integer $n$ counts the positive integers which are relatively prime with and less than or equal to $n$.  So, for example, for an odd prime $p$, the positive integers which are not relatively prime with $p^r$ are $p$, $2p$, $3p$, \ldots $p^{r-1} p = p^r$.  There are $p^{r-1}$ of these.  All other positive integers are relatively prime with $p^r$, hence $\phi(p^r) = p^r - p^{r-1}$.  The function $\psi$ is further defined as follows:
$$\psi(1) = \psi(2) = 0, \,  \psi(2^r) :=\phi(2^r) \textrm{ for } r>1,$$
and
$$\textrm{ for } m = \prod_i p_i ^{r_i}, \quad  \psi(m) := \sum_i \psi(p_i ^{r_i}).$$

\begin{theorem}[Crystallographic Restriction I]  For any discrete group $G$ of isometries of $\mathbb{R}^n$, for $n \geq 2$, the set of orders of the elements $G$ which have finite order is equal to
$$\operatorname{Ord}_n = \{ m \in \mathbb{N} : \psi(m) \leq n \}.$$
\end{theorem}

The crystallographic restriction theorem is connected to the mathematics of crystals when we reformulate the theorem in the context of lattices.  A full-rank lattice is a set of points in $\mathbb{R}^n$ of the form
\begin{equation} \label{lattice1} \Gamma = \{ \bp \in \mathbb{R}^n : \bp = L \bx, \quad L \in \operatorname{GL}(n, \mathbb{R}), \quad \bx \in \mathbb{Z}^n \}. \end{equation}
Above, $\operatorname{GL}(n, \mathbb{R})$ is the set of $n \times n$ invertible matrices with real entries, and $\mathbb{Z}^n$ are the elements of $\mathbb{R}^n$ whose entries are integers.  We say that the matrix $L$ generates the lattice $\Gamma$.  The generating matrix $L$ is not unique, because for any $M \in \operatorname{GL} (n , \mathbb{Z})$ the set of points in \eqref{lattice1} is equal to
$$ \{ \bp \in \mathbb{R}^n : \bp = LM \bx, \quad \bx \in \mathbb{Z}^n \}.$$
Here $\operatorname{GL}(n, \mathbb{Z})$ is the group of invertible $n \times n$ matrices whose entries are integers. Note that to be a group, this requires the determinant of all elements of $\operatorname{GL}(n, \mathbb{Z})$ to be equal to $\pm 1$.  Two matrices $L_1, L_2 \in \operatorname{GL}(n, \mathbb{R})$ generate the same lattice if and only if there is an $M \in \operatorname{GL}(n, \mathbb{Z})$ such that $L_1 = L_2 M$.  For a matrix $M \in \operatorname{GL}(n, \mathbb{Z})$, we identify it with the isometry of $\mathbb{R}^n$ which maps $\bx \in \mathbb{R}^n$ to $M\bx$. The matrices in $\operatorname{GL}(n, \mathbb{Z})$ can therefore be identified with the group of symmetries of the crystal whose atoms lie on the points of the lattice.   Hence, the order of $M$ is equal to the smallest positive integer $k$ such that $M^k$ is the identity matrix.  It turns out that the set of orders of the elements of any discrete group $G$ of isometries of $\mathbb{R}^n$ which have finite order is equal to the set of orders of the elements of $\operatorname{GL}(n, \mathbb{Z})$.  Consequently, the crystallographic restriction may be reformulated as follows.

\begin{theorem}[Crystallographic Restriction II] \label{th:crystal} For any $n \geq 2$, the set of orders of the elements of $\operatorname{GL}(n, \mathbb{Z})$ is equal to
$$\operatorname{Ord}_n = \{ m \in \mathbb{N} : \psi(m) \leq n \}.$$
\end{theorem}

In \cite{goldbach}, Bamberg, Cairns and Kilminster proved that one may reformulate the Strong Goldbach Conjecture in terms of the orders of elements of $\operatorname{GL} (n, \mathbb{Z})$.

\begin{conj}[Strong Goldbach] Every even natural number greater than six can be written as the sum of two distinct odd primes.
\end{conj}

\begin{theorem}[Theorem 3 of \cite{goldbach}] \label{th:gold} The following statements are equivalent:
\begin{enumerate}
\item The strong Goldbach conjecture is true;
\item For each even $n \geq 6$ there is a matrix $M \in \operatorname{GL}(n, \mathbb{Z})$ which has order $pq$ for distinct primes $p$ and $q$, and there is no matrix in $\operatorname{GL}(k, \mathbb{Z})$ of order $pq$ for any $k < n$.
\end{enumerate}
\end{theorem}

The Goldbach Conjecture is an extremely difficult problem.  Difficult, long-standing open problems have sometimes been solved by translating the problem into a different field of mathematics.  The proof of Fermat's Last Theorem, also a statement in number theory, was achieved using newly developed techniques in algebraic geometry \cite{wiles1}, \cite{wiles2}.  To approach the Goldbach Conjecture geometrically, we ask
\begin{question} \em Is there a geometric reason for the existence of a symmetry for full rank lattices in $\mathbb{R}^n$, with $n\geq 6$ an even number, such that this symmetry is of order $pq$ for two odd primes $p \neq q$ such that $p+q=n+2$?
\em
\end{question}
The condition that there is no matrix in $\operatorname{GL}(k, \mathbb{Z})$ of order $pq$ for any $k<n$ is equivalent to requiring $p+q = n+2$.  This follows from the Crystallographic Restriction Theorem \ref{th:crystal} which states that the orders of the elements of $\operatorname{GL}(k, \mathbb{Z})$ is equal to the set of non-negative integers $m$ with $\psi(m) \leq k$.  In order to guarantee that
$$\psi(pq) = p+q-2 > k \textrm{ for all } k < n, \textrm{ but } \psi(pq) \leq n \implies \psi(pq) = p+q-2 = n.$$

Consequently, the symmetry of order $pq$ would correspond to a matrix $M \in \operatorname{GL}(n, \mathbb{Z})$ which does not admit a diagonal decomposition into two matrices of smaller dimensions.  Geometrically, this matrix would not arise as a product of symmetries of $\mathbb{R}^k$ and $\mathbb{R}^{n-k}$ for any $k=1, \ldots, n-1$.  It would be a new symmetry occurring first in $\mathbb{R}^n$.  Since \cite{goldbach} already realized the connection between the Goldbach Conjecture and the Crystallographic Restriction Theorem, this geometric approach is pretty unlikely to lead to any new developments.  Nonetheless, it is interesting that a famous number theoretic conjecture can be equivalently phrased as a simple question about the orders of symmetries of full-rank lattices in $\mathbb{R}^n$.


\section{Introduction}
Neutrinos produced in high-energy extraterrestrial acceleration sites travel large distances not scattered by either interstellar or intergalactic medium and undeflected by cosmic magnetic fields. Information about neutrino direction, energy and time of arrival at Earth combined with results of gamma-ray, optical, cosmic ray and gravitational wave observations can constrain models of high-energy cosmic particle production mechanisms. Back in 1960 M.A. Markov proposed to instrument natural water reservoirs with sparse arrays of photodetectors to detect upgoing charged particles by means of Cerenkov radiation \cite{markov}. Neutrino interactions near or inside the instrumented water volume can be detected by reconstructing the cascade of charged particles accompanying the interaction. A larger detection range is accessible for $\nu_{\mu}$ and partly $\nu_{\tau}$ interactions mediated by W$^{\pm}$ which are accompanied by the muon track due to long range of high energy muons. This technique allows to increase the sensitive volume of the detector to the order of 1~km$^3$ making it a good choice for the search of small flux of high-energy astrophysical neutrinos. There are three large-scale neutrino telescopes operating presently: ANTARES telescope in the Mediterranean sea, IceCube in the Antarctic Ice Sheet, and Baikal-GVD in Lake Baikal. In addition work on constructing ANTARES successor KM3NET telescope is ongoing in the Mediterranean. IceCube was the first detector to reach instrumented volume of 1km$^3$. Sensitivity of these devices extends from 10's - 100's of GeV to multi-PeV neutrino energies and beyond with lower bound depending on density of photodetectors and ice or water light absorption length. In this energy range most of neutrinos reaching the detectors are produced in cosmic ray interactions with Earth's atmosphere. The IceCube experiment demonstrated an excess of neutrino events over the atmospheric neutrino spectrum starting from $\sim$100 TeV in 2011 - 2013 \cite{ic_diff_discov}, in 6 year data analysis released in 2016 the significance of an excess has reached 8$\sigma$ \cite{ic_review}. Presently these events are attributed to neutrinos produced in high-energy processes in remote cosmic objects, however any indications of astrophysical neutrino sources have not yet reached the 5$\sigma$ discovery level. An important step towards understanding of cosmic neutrino origin was the detection of neutrino alert event associated with blazar TXS0506+056 \cite{ic_txs}, analysis of IceCube events prior to alert had shown 3.5$\sigma$ evidence of neutrino emission from that blazar direction \cite{ic_txs_prior}. The recent point source search with 10 year data of IceCube track-like events indicate 3.3$\sigma$ inconsistency with background for the full Nothern Hemisphere source catalog \cite{ic_10y}. Moreover recent studies indicate 4.1$\sigma$ association of IceCube track-like events with radio-bright blazars selected on the basis of very long baseline interferometry observations \cite{plavin2}.

Work on construction of the neutrino telescope in Lake Baikal began in 1980. In the mid-1990s first underwater neutrinos were detected consecutively at the NT-36 and NT-96 neutrino telescopes with 36 and 96 photodetectors respectively \cite{nt96_perf,nt96_neutrino}. Work on construction of the Baikal-GVD telescope of the cubic-kilometer scale began in 2011, and the first demo cluster named Dubna was commissioned in 2015. Since April 2020, the telescope includes seven clusters carrying 2016 optical modules (OMs). The present work reports on the design and status of the Baikal-GVD detector and on first results of data analysis.

\section{Baikal-GVD detector}
The experiment site is located in the southern part of Lake Baikal 3.6 km away from shore where the lake depth is nearly constant at 1366-1367 m. Telescope consists of independent structural units - clusters (Fig.~\ref{fig:det}). Cluster consists of 8 strings each carrying 36 optical modules (OM) and calibration systems.
OMs are located at depths from 750 to 1275 (m) with a vertical step of 15 m. Optical module is a glass sphere which hosts photomultiplier tube (PMT) Hamamatsu R7081-100 with a 10-inch hemispherical photocathode oriented towards the lake floor along with various sensors, readout and HV control electronics, and calibration LED. OM readout is organised in sections each including 12 OMs. Readout controllers are housed in a section's Central Module (CM) in which digitisation of PMT signal occurs with a step of 5 ns. Time of signal arrival and charge deposited in PMT is derived from the pulse shape analysis. CM generates local trigger signal if pulses with deposited charge above adjustable thresholds $Q_{high}$ and $Q_{low}$ coincident within 100 ns time window are found in adjacent modules within one section.
\begin{wrapfigure}[20]{r}{0.4\textwidth}
\centering
    \includegraphics[width=0.4\textwidth]{baikalGVD_detector.pdf}
    \caption{General view of Baikal-GVD detector in 2020. Cluster radius is 60 m and the distance between neighbouring cluster centers is $\sim$ 300 m. A yearly progression of the detector deployment is shown in the legend.}
    \label{fig:det}
\end{wrapfigure}
Typical values for $Q_{high}$ are 3-5 photoelectrons (p.e.) and 1-2 p.e. for $Q_{low}$ depending on the cluster and season. If a local trigger is generated within one of CM an event frame of 5~$\mu$s is read out from all CMs of the cluster and sent to the shore center via optoelectric cable attached to a central cluster control module. Data are further transmitted to JINR (Dubna) for event reconstruction and long-term storage.
Baikal water currents cause devations of topmost OMs up to 50 m from their median position with average speed of 0.5 cm/s. For precise OM positioning each string carries system of acoustic modems, a cross-poll of modems if performed each 1 - 6 minutes. OM coordinates are reconstructed with the precision better than 20~cm \cite{acoustic}. Single cluster time calibration is performed with the precision of 2.5 ns using OM integrated LEDs and LED matrices which produce light flashes propagating to up to 100 m radius~\cite{tcalib}. A system of three technological strings located between clusters carries 5 dedicated lasers producing isotropical flashes which are used for intercluster calibration and water properties monitoring.
The light absorption length of Baikal water amounts to $\sim$22~m.
The PMT noise rate is dominated by chemiluminescence, with substantial seasonal variations. Noise pulses at the level of 1 p.e. are produced at the rate of 20-50 kHz in quiet period in April - June and up to > 100 kHz at topmost OMs for some periods in the rest of the year \cite{baikal_noise_icrc}. The noise rate is monitored and is taken into account in the detector Monte-Carlo (MC) simulations.

\section{Physics results}

Presently the data are analysed for each of 7 clusters idependently what is referred to as single cluster analysis while the multi cluster analysis is in development.
The muon track reconstruction is performed in two stages. At first stage PMT pulses are clustered using the time and distance constraints with respect to the seed pulse with high charge deposition. A preliminary muon track parameters estimation is performed and outlier pulses are removed. As a results of this procedure the collection of pulses with average noise contribution at the level of 1\% is selected.
At the second stage the muon track parameters are reconstructed by means of minimisation of the quality function $Q=\chi^{2}(t) + w*f(q,r)$ where $\chi^{2}(t)$ is the chi-square sum of time residuals with respect to direct Cerenkov light from the muon, $f(q,r)$ is the sum of products of charges deposited in OMs and their distances from the track and $w$ is the relative weight of the second term.
\begin{wrapfigure}[18]{l}{0.4\textwidth}
\centering
  \includegraphics[width=0.4\textwidth]{muon_flux.pdf}
  \caption{Reconstructed muon zenith angle distribution in cluster 1 of 2016 season. Data are compared to MC simulation based on CORSIKA event generator~\cite{corsika}.}
\label{fig:muons}
\end{wrapfigure}
Events with at least 6 pulses at 2 strings are used for the reconstruction.
Distribution of zenith angle for tracks reconstructed in quiet period of season 2016 with minimal quality requirements is shown at Fig.~\ref{fig:muons}. A few~\% of events are reconstructed as upgoing muons forming a background to upgoing neutrino events exceeding the atmospheric neutrino flux by a factor $10^5-10^6$. To suppress the background from misreconstructed tracks a set of cuts on 13 reconstructed track parameters is used with the most powerful variable being $Q/ndf$.
For the single cluster analysis the zenith angle is constrained to $\Theta_{zenith} > 120^{\circ}$ to ensure sufficient length of the track for high quality reconstruction. Median angular resolution for events passing the analysis cuts is $\sim$1.0$^{\circ}$.
The developed cuts were applied to a collection of runs from April 1st to June 30th 2019 (quiet period).
The total live time of selected sample is 323 days of single cluster data-taking. A total of 57 neutrino candidate events were selected while the expectation from atmospheric neutrino MC simulation is 54.3 events (Fig.~\ref{fig:neutrino}).
\begin{figure}[ht]
\centering
 \subfigure[]
 {
    \includegraphics[width=0.4\textwidth]{cands_theta_sep20_clusSum.pdf}
    }
 \subfigure[]
 {
    \includegraphics[width=0.17\textwidth]{cand_cl4.pdf}
	}
 \subfigure[]
 {
    \includegraphics[width=0.17\textwidth]{cand_cl5.pdf}
	}
    \caption{(a): Zenith angle distribution of 57 neurino candidate events selected in April - June 2019. An expectation for backgound atmospheric muons was estimated to < 1 events and thus is not shown. (b): Neutrino candidate reconstructed on June 11th 2019 in cluster 1, $\Theta_{zenith}=161.8^{\circ}$. (c): Neutrino candidate reconstructed on June 18th 2019 in cluster 5, $\Theta_{zenith}=148.0^{\circ}$. On event displays early hits are shown as red and late hits in blue. }
    \label{fig:neutrino}
\end{figure}
\newpage

Cascade event energy and direction are reconstructed using maximum likelihood technique where time and charge of pulses are taken into account \cite{baikal_casc_icrc}. For the cascade of energy 100 TeV the uncertainty of cascade direction reconstruction amounts to $\sim$4$^{\circ}$ while the energy reconstruction uncertainty is $\delta E/E \sim 30\%$.
The telescope effective volume for detecting high-energy cascades was 0.25 km$^{3}$ in 2019 and became 0.35 km$^3$ in 2020. Based on the IceCube fit of the astrophysical neutrino spectrum \cite{ic_spectrum_modern} 0.6 neutrino events are expected per year per cluster or 4.2 events per year for detector configuration in 2020. The high-energy cascade search was performed using the data of years 2015, 2016, 2018 and 2019. The sum live time for the used data is 2294 days of single cluster data-taking.
Events with reconstructed energies above 100 TeV with at least 19 fired OMs were selected as neutrino candidate events. In 2015-2019 data twelve such events were found, all of these events were reconstructed as downgoing. In addition one upgoing event with 19 hits and reconstructed energy $\sim$91 TeV was found in 2019 data.

On August 17, 2017 the LIGO and VIRGO detectors of gravitational waves recorded a signal named GW170817 which was folllowed by gamma and optical bursts \cite{gw170817}.
This event was associated with the merger of neutron stars in the NGC 4993 galaxy.
The ANTARES and IceCube telescopes, the Pierre-Auger observatory and the Super-Kamiokande detector performed searches for neutrino associated with this event \cite{gw170817_neutrino,gw_superk}.
In Baikal-GVD the search for neutrinos from GW170817 was performed in the cascade channel \cite{gw_gvd}. In 2017 Baikal-GVD consisted of two clusters. At the time of the event the NGC 4993 galaxy was approximately at the horizon ($\Theta_{zenith}=93.3^{\circ}$). The search was conducted in the time window of $\pm$500 s and $\pm$14 days around the event to take into account all neutrino production models. The requirement on the number of fired optical modules was relaxed but no event in the direction of of NGC 4993 was detected. A limit  on the neutrino flux at a 90\% C.L. was set under the assumption of $E^{-2}$ spectrum and identical fluxes of all neutrino flavors (Fig.~\ref{fig:gw}).
\begin{wrapfigure}[18]{r}{0.45\textwidth}
\centering
  \includegraphics[width=0.45\textwidth]{gw.pdf}
  \caption{Limits on the neutrino flux from GW170817 event. Figure is taken from~\cite{gw_gvd}.}
\label{fig:gw}
\end{wrapfigure}

\section{Conclusions}
The Baikal-GVD deployment is ongoing with a pace of 2 clusters per year. First results from uncompleted detector are available. A first set of track-like neutrino candidate events obtained using the data from April - June of 2019 was presented. An acceptable agreement of track-like neutrino event rate with MC expectation was demonstrated. A set of twelve downgoing high energy cascade events with estimated energy above 100 TeV and an upgoing cascade event with estimated energy of $\sim$91 TeV were selected in 2015-2019 data. A limit on neutrino flux from neutron star merger event GW170817 was set.

This work is supported by the Ministry of Science and Higher Education of Russian Federation under the contract 075-15-2020-778 in the framework of the Large scientific projects program within the national project "Science".
We acknowledge the support from RFBR grant 20-02-00400 and grant 12-29-11-029, as well as support by the JINR young scientist and specialist grant 20-202-09. We also acknowledge the technical support of JINR staff for the computing facilities (JINR cloud).