Consider the Neumann eigenvalue problem
\[\left\{ \begin{aligned} -\Delta \, u &= \mu u &&\text{in } \Omega, \\ \frac{\partial u}{\partial \nu} &= 0 &&\text{on } \partial \Omega, \end{aligned} \right.\]where $\Omega \subset \mathbb{R}^2$ is a bounded convex domain. Recently, Filonov in [1] obtained the following lower bound on the eigenvalue counting function:
\[N_{\mathcal{N}}(\Omega,\mu) \geq \frac{|\Omega| \lambda}{2 \sqrt{3} j_0^2},\]where $j_0^2$ is the first positive zero of the Bessel function $J_0$.
Sketchily, the approach of [1] is the following: first we densely pack equal disks in $\mathbb{R}^2$, and then choose those whose centers lie in $\Omega$. If a disk $B$ is a subset of $\Omega$, then we consider the first Dirichlet eigenfunction in $B$. If $B \setminus \Omega$ is nonempty, then we consider the restriction of the first Dirichlet eigenfunction in $B$ to $\Omega \cap B$. Using these functions, we construct the test subspace and estimate $\mu_k(\Omega)$ from above by a factor coming from the number of disks and the first Dirichlet eigenvalue $\lambda_1(B)$, which leads to a required lower bound for $N_{\mathcal{N}}(\Omega,\mu)$.
The tricky point here is the consideration of the case when $B \setminus \Omega$ is nonempty. In this case, roughly speaking, one needs to justify that
\[\tau_1(\Omega \cap B) \leq \lambda_1(B),\]where $\tau_1(\Omega \cap B)$ is the first eigenvalue in $\Omega \cap B$ under the zero Dirichlet boundary conditions on $\overline{\Omega} \cap \partial B$ and zero Neumann boundary conditions on the remaining part of $\partial (\Omega \cap B)$. This fact follows from Lemma 2.1 in [1] which states a certain integral property of Bessel functions. Here the convexity of $\Omega$ is employed. (Note that the fact remains true if $\Omega$ is merely star-shaped with respect to the center of $B$. However, it is hard to weaken the convexity in general, since the position of $B$ with respect to $\Omega$ is not given constructively.)
It is tempting to anticipate that one could substitute disks by hexagons in the approach above, and hence improve the upper bound for $\mu_k(\Omega)$, thereby improving the lower bound for $N_{\mathcal{N}}(\Omega,\mu)$. To do it rigorously, one has to prove the inequality
\[\tau_1(\Omega \cap H) \leq \lambda_1(H),\]where $H$ is a hexagon such that $H \setminus \Omega \neq \emptyset$.
Unfortunately, it seems that the inequality cannot be true, in general. Let $H$ be a hexagon with the side $1$ centered at $(0,0)$, and let $\Omega$ be a large triangle spanned on the points $(0,0)$, $(-20,-2)$, $(20,-2)$, see figure below.
Mathematica gives the following values for the corresponding eigenvalues:
\[\tau_1(\Omega \cap H) \approx 7.20569 \quad \text{and} \quad \lambda_1(\Omega) \approx 7.15548.\]I admit that calculations might be completely wrong for $\tau_1(\Omega \cap H)$, but $\lambda_1(\Omega)$ is calculated more-less ok. Playing with parameters (vertices of the triangle $\Omega$), I also observe continuous dependence of $\tau_1(\Omega \cap H)$ on them. The inequality holds for some parameters, and does not hold for others. This indirectly indicates that the calculation can be reliable.
Conclusion: it is not that easy to enhance the estimate from [1] using the same strategy.