$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\supp}{\operatorname{supp}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$ $\newcommand{\mes}{\operatorname{mes}}$
We are interested how eignevalues are distributed. We introduce eigenvalue distribution function $N(\lambda)$ which is the number of eigenvalues (of operator $-\Delta$) which are less than $\lambda$. In other words, \begin{equation} N(\lambda)=\max _{\lambda_n< b} n \label{eq-13.2.1} \end{equation} where $\lambda_n$ are eigenvalues, or, using Corollary 13.1.1 we reformulate it as
Definition 1. Eignevalue distribution function is \begin{equation} N(\lambda)=\max_{L\subset \mathsf{H}:\ Q(v)<\lambda \|v\|^2\ \forall 0\ne v\in L }\ \dim L \label{eq-13.2.2} \end{equation}
Remark 1. While these two definitions coincide in our particular case there is a huge difference: formula (\ref{eq-13.2.1}) counts only eigenvalues. However formula (\ref{eq-13.2.2}) takes into account also a continuous spectrum: this way defined $N(\lambda)$ is $+\infty$ if $(-\infty,\lambda)$ contains (some) points of continuous spectrum. This definition is one widely used.
In the same way as Theorem 13.1.5 and Theorem 13.1.6 one can prove
Theorem 1.
Consider rectangular box in $\mathbb{R}^d$: $\Omega =\{x:\, 0< x_1< a_1, 0< x_2< a_2, \ldots, 0< x_d< a_d\}$ with Dirichlet or Neumann boundary conditions. Then separation of variables brings us eigenfunctions and eigenvalues \begin{align} &X_m=\sin \Bigl(\frac{\pi m_1x_1}{a_1}\Bigr) \sin \Bigl(\frac{\pi m_2x_2}{a_2}\Bigr) \cdots \sin \Bigl(\frac{\pi m_dx_d}{a_d}\Bigr),\label{eq-13.2.3}\\ &\lambda_m= \pi^2\Bigl(\frac{m_1^2}{a_1^2}+\frac{m_2^2}{a_2^2}+\ldots +\frac{m_d^2}{a_d^2}\Bigr) \qquad m_1\ge 1, \ldots, m_d\ge 1 \label{eq-13.2.4} \end{align} for Dirichlet problem and \begin{align} &X_m=\cos \Bigl(\frac{\pi m_1x_1}{a_1}\Bigr) \cos \Bigl(\frac{\pi m_2x_2}{a_2}\Bigr)\cdots \cos \Bigl(\frac{\pi m_dx_d}{a_d}\Bigr),\label{eq-13.2.5}\\ &\lambda_m= \pi^2\Bigl(\frac{m_1^2}{a_1^2}+\frac{m_2^2}{a_2^2}+\ldots +\frac{m_d^2}{a_d^2}\Bigr) \qquad m_1\ge 0, \ldots, m_d\ge 0 \label{eq-13.2.6} \end{align} for Neumann problem where $m=(m_1,\ldots,m_d)\in \mathbb{Z}^d$ in both cases.
Exercise 1. Prove it by separation of variables.
Therefore
Theorem 2. For rectangular box $\Omega =\{x:\, 0< x_1< a_1, 0< x_2< a_2, \ldots, 0< x_d< a_d\}$
Calculation of the number of integer points in the "large" domains is an important problem of the Number Theory. For us the answer "the number of integer points inside of the large domain approximately equals to its volume" is completely sufficient but let us use the following
Theorem 3. As $d\ge 2$ the number of integer points inside ellipsoid \begin{equation*} \mathcal{E}(\lambda)=\Bigl\{m=(m_1,\ldots,m_d):\, \frac{m_1^2}{a_1^2}+\frac{m_2^2}{a_2^2}+\ldots +\frac{m_d^2}{a_d^2}<\pi^{-2}\lambda\Bigr\} \end{equation*} equals to the volume of $\mathcal{E}(\lambda)$ plus $o(\lambda^{(d-1)/2})$.
Remark 2. Actually much more precise results are known.
Observe that the volume of this ellipsoid $\mathcal{E}(\lambda)$ is $\pi^{-d}\omega_d a_1a_2\cdots a_d\lambda^{d/2}$ where $\omega_d$ is a volume of the unit ball in $\mathbb{R}^d$.
Exercise 2. Prove it observing that ellipsoid $\mathcal{E}(\lambda)$ is obtained by stretching of the unit ball.
Observe also that both $E_+(\lambda)$ and $E_-(\lambda)$ constitute $1/2^d$ part of $\mathcal{E}(\lambda)$ and therefore their volumes are $(2\pi)^{-d}\omega_da_1\cdots a_d \lambda^{d/2}$. However if we consider integer points we observe that
Actually, it is not completely correct as one needs to take into account intersections of $\mathcal{E}(\lambda)$ with two different planes $\Pi_j$ and $\Pi_k$ but the number of integer points there is $O(\lambda^{(d-2)/2})$.
Since intersections of $\mathcal{E}(\lambda)$ with $\Pi_j$ is an ellipsoid we conclude that the number of integer points there is $\pi^{1-d}\omega_{d-1}a_1a_2 \cdots a_d/a_j\lambda^{(d-1)/2}$.
Therefore the number of integer points in $E_\pm (\lambda)$ is equal to \begin{multline*} (2\pi)^{-d}\omega_d a_1a_2\cdots a_d \lambda^{d/2}\pm\\ \frac{1}{2}(2\pi)^{1-d}\omega_{d-1} \underbracket{a_1a_2\cdots a_d \bigl(a_1^{-1}+a_2^{-1}+\ldots a_d^{-1}\bigr) } \lambda^{(d-1)/2} + o\bigl(\lambda^{(d-1)/2}\bigr). \end{multline*} Observe that $a_1a_2\cdots a_d$ is a volume of $\Omega$ and $a_1a_2\cdots a_d \bigl(a_1^{-1}+a_2^{-1}+\ldots a_d^{-1}\bigr)$ is a half of area of its boundary $\partial\Omega$.
So we arrive to
Theorem 4. For rectangular box $\Omega =\{x:\, 0< x_1< a_1, 0< x_2< a_2, \ldots, 0< x_d< a_d\}$, $d\ge 2$ \begin{equation} N(\lambda)=(2\pi)^{-d} \mes_d(\Omega)\lambda^{d/2}+ O\bigl(\lambda^{(d-1)/2}\bigr) \label{eq-13.2.7} \end{equation} and more precisely \begin{equation} N (\lambda)=(2\pi)^{-d} \mes_d(\Omega)\lambda^{d/2} \pm \frac{1}{4} (2\pi)^{1-d} \mes_{d-1}(\partial\Omega)\lambda^{(d-1)/2} +o\bigl(\lambda^{(d-1)/2}\bigr) \label{eq-13.2.8} \end{equation} where "$+$" corresponds to Neumann and "$-$" to Dirichlet boundary conditions and $\mes_k$ means $k$-dimensional volume.
Remark 3.
Now we want to generalize these results for an arbitrary domain $\Omega$. To do so let us consider domain $\Omega^+_\varepsilon$ which includes $\Omega$ and all points on the distance $\le C\varepsilon$ from $\partial \Omega$ and partition it into cubic boxes such that
The boxes which are completely inside of $\Omega$ are called inner boxes and the boxes which intersect $\partial \Omega$ are called boundary boxes.
We consider now only Dirichlet conditions. Let us use Definition 1 for $N(\lambda)$; for Dirichlet boundary conditions we can consider $\mathsf{H}$ as the space of all functions which vanish outside of $\Omega$.
Let us tighten constrains to these functions: we assume that they are $0$ in the boundary boxes as well and that they vanish on the walls of all boxes. Then $N(\lambda)$ can only decrease: $N(\lambda)\ge N_*((\lambda)$ where $N_*(\lambda)$ is an eigenvalue counting function under these new constrains. But under these constrains calculation in each box becomes independent and we arrive to \begin{equation*} N_*(\lambda)= \Bigl[(2\pi)^{-d}\omega_d \lambda^{d/2}- c \lambda^{(d-1)/2} \varepsilon^{-1}\Bigr] \sum_\iota \mes_d (B_\iota) \end{equation*} where we sum over all inner boxes. As $\varepsilon \ge c_1 \lambda ^{-1/2}$ this expression is greater than $\Bigl[(2\pi)^{-d}\omega_d \lambda^{d/2}- c \lambda^{(d-1)/2} \varepsilon^{-1}\Bigr] \mes_d (\Omega^-_\varepsilon) $ where $\Omega^-_\varepsilon$ is the set of all points of $\Omega$ on the distance $> \varepsilon$ from $\partial\Omega$. Thus \begin{equation} N_D(\lambda)\ge \Bigl[(2\pi)^{-d}\omega_d\lambda^{d/2}- c\lambda^{(d-1)/2} \varepsilon^{-1}\Bigr] \bigl(\mes_d (\Omega) - \mes_d (\Omega \setminus \Omega^-_\varepsilon)\bigr). \label{eq-13.2.9} \end{equation}
Let us loosen constrains to these functions: we assume that they are arbitrary in the boundary boxes as well and that they can be discontinuous on the walls of all boxes (then $Q(u)$ is calculated as a sum of $Q_\iota (u)$ all over boxes). Then $N(\lambda)$ can only increase: $N(\lambda)\le N^*((\lambda)$ where $N^*(\lambda)$ is an eigenvalue counting function under these new constrains. But under these constrains calculation in each box becomes independent and we arrive to \begin{equation*} N_*(\lambda)= \Bigl[(2\pi)^{-d}\omega_d \lambda^{d/2}+ c \lambda^{(d-1)/2} \varepsilon^{-1}\Bigr] \sum_\iota \mes_d (B_\iota) \end{equation*} where we sum over all inner and boundary boxes. Thus \begin{equation} N_D(\lambda)\le \Bigl[(2\pi)^{-d}\omega_d \lambda^{d/2}+ c \lambda^{(d-1)/2} \varepsilon^{-1}\Bigr] \bigl(\mes_d(\Omega) + \mes_d(\Omega ^+_\varepsilon \setminus \Omega )\bigr). \label{eq-13.2.10} \end{equation} But what is the gap between the right-hand expressions of (\ref{eq-13.2.9}) and (\ref{eq-13.2.10})? It does not exceed \begin{equation*} C\lambda^{(d-1)/2}\mes_d \varepsilon^{-1} (\Omega ^+_\varepsilon ) + C\lambda^{d/2} \mes_d (\Sigma_\varepsilon ) \end{equation*} where $\Sigma_\varepsilon$ is the set of all points on the distance $\le C_0\varepsilon$ from $\Sigma=\partial \Omega$.
We definitely want $\mes_d (\Sigma_\varepsilon )\to 0$ as $\varepsilon\to +0$. According to Real Analysis we say that it means exactly that $\Sigma$ is a set of measure $0$.
Thus we arrive to
Theorem 5. If $\Omega$ is a bounded domain and $\partial \Omega$ is a set of measure $0$ then \begin{equation} N_D(\lambda)= (2\pi)^{-d}\omega_d \lambda^{d/2}+ o\bigl(\lambda^{d/2}\bigr) \label{eq-13.2.11} \end{equation} as $\lambda\to +\infty$.
In particular, it holds if $\Omega$ is a bounded domain and $\partial \Omega$ is smooth.
Remark 4.
Remark 5.
Remark 6. Let $\Omega$ be bounded domain and $\partial \Omega$ be smooth.
Remark 7. There are many generalizations of these asymptotics; let us mention only few. We consider Schrödinger operator \begin{equation} H=-\hbar^2\Delta + V(x). \label{eq-13.2.16} \end{equation} Then Weyl formula for this operator is \begin{equation} N(\lambda,\hbar)\approx (2\pi\hbar)^{-d}\omega_d \iiint (\lambda-V(x))_+^{d/2}\,dx \label{eq-13.2.17} \end{equation} where $z_+=\max(z,0)$. This formula is justified