Helmholtz equations with strict Neumann boundary conditions [RESOLVED] I've been looking for answers to this for some time with no avail. I have come across this boundary value problem in my research which I solved in a spherically symmetric case, but I'm curious about its more general properties. I'm looking for solutions to the Helmholtz equation in a closed domain $\Omega$ with Neumann boundary conditions:
$\nabla^2 u = -k^2 u$, $\nabla u |_{\partial \Omega}= 0$
Notice that this isn't the usual Neumann problem $\partial u/\partial n |_{\partial \Omega} = 0$, in this case all derivatives are $0$ at the boundary. In the spherical case, this simply forces $u = \sum_n a_n j_0(r\alpha_{0n}/R)$, where $j_0$ is the 0th order spherical Bessel function of the first kind, and $\alpha_{0n}$ is the $n$th zero of its derivative. This gives me the eigenvalue $k_n = \alpha_{0n}/R$. I was surprised to find that there is no angular dependence at all here, the solutions are effectively one-dimensional along the radial direction.
After trying this in another easy geometry (between two infinite parallel plates) and finding that similar things happen (the gradient of $u$ is purely normal to the plates), I was wondering if there was a nice way to state that if I impose $\nabla u = 0$ on the boundary of a simple closed domain (let's say it has to be homeomorphic to a sphere), it essentially kills all dependence except for the "radial" direction. Is there such a theorem that I am missing?
Thanks!
Edit: After considering the method shown below by @Jap88, I realized that this problem can be re-cast in a way which may be more familiar to the usual mathematics literature. The usual Neumann boundary condition is $\partial u/\partial n |_{\partial \Omega} = 0$. In our BC, we have the added condition that the tangential derivative must also be 0 everywhere, which is equivalent to making $u$ a constant on the boundary. So we can consider the problem:
$\nabla^2 u = -k^2 u,\; \hat n\cdot \nabla u|_{\partial \Omega} = 0,\; u|_{\partial \Omega} = C$
Which I think in a generic geometry is an overdetermined problem. I guess in that case it turns out that there are some special geometries (like the sphere) where it just so happens that some solutions meet both boundary conditions.
Edit 2
After more time spent on this problem, I have convinced myself the boundary conditions I am imposing here are not only mathematically over-constraining, they are simply unphysical for the system of interest. In reality the normal boundary condition is a bit more complex and would require a deeper explanation of the underlying model, while the tangential BCs remain simple Neumann boundary conditions.
 A: Perhaps casting the equation on weak form can be illuminating.
Multiply the equation with a test function $v$, see https://en.wikipedia.org/wiki/Weak_formulation.
This gives:
$$\int_\Omega (\nabla \cdot \nabla u) v d V +\int_\Omega k^2 u v d V=0$$
Now, using partial integration (divergence theorem), we get:
$$\int_{\partial \Omega} (\mathbf{n} \cdot \nabla u) v d S -\int_\Omega (\nabla  u \cdot \nabla v) d V+ \int_\Omega k^2 u v dV=0 $$
This is one way to see how the surface normal appears naturally in this context.
So, the only part of the gradient at the boundary that directly influences what is going on in the volume is the normal component of the gradient.
Intuitively, this is since the Laplace operator operating on a field in a control volume encodes what is flowing in and out of that control volume only in the normal direction. The Laplace operator doesn't "sense" the tangential component.
More Details:
I think one can reason like this, but the arguments are incomplete and more of an "engineering" argument I think:
Ignoring the Helmholtz term, since this discussion is only related to the Laplace operator, one can see the entire Laplace problem as minimization of potential energy (Dirichlet's principle):
$$L(u)=\frac{1}{2}\int_\Omega (\nabla u \cdot \nabla u) d V -\int_{\partial \Omega} uq d S$$
Now, we can add a general constraint to this by means of Lagrange multipliers:
$$L(u)=\frac{1}{2}\int_\Omega (\nabla u \cdot \nabla u) d V -\int_{\partial \Omega} uq d S + \int_{\partial \Omega}\lambda \cdot \nabla u dS$$
where $\lambda$ is a vector of same length as $\nabla u $. By modifying $\lambda$, which is a vector field, we can get any constraint we want.
To find the minimum we seek a stationary point to this system. To find this we set the first variation (first "functional derivative") to zero:
$$\delta L=0$$
("Differentiating" with respect to $u$.)
This will lead us something similar to the above weak form with $v=\delta u$ and $\nabla v = \nabla \delta  u=\delta \nabla u$.
The first two terms in the energy expression will, after minimizing, lead back to Laplace's equation:
$$\nabla \cdot \nabla u =0$$
with the Neumann boundary condition
$$\mathbf{n} \cdot \nabla u =q$$
The Lagrange term may participate in the divergence theorem and affect the solution.
To see this we can split the Lagrange term in a normal and tangential component:
$$L(u)=\frac{1}{2}\int_\Omega (\nabla u \cdot \nabla u) d V -\int_{\partial \Omega} uq d S + \int_{\partial \Omega}\lambda_n \cdot \nabla u dS+\int_{\partial \Omega}\lambda_t \cdot \nabla u dS$$
If we now focus on the boundary we have:
$$-\int_{\partial \Omega} uq d S + \int_{\partial \Omega}\lambda_n \cdot \nabla u dS+\int_{\partial \Omega}\lambda_t \cdot \nabla u dS$$
or
$$-\int_{\partial \Omega} u \mathbf{n} \cdot \nabla u d S + \int_{\partial \Omega}|\lambda_n| \mathbf{n} \cdot \nabla u dS+\int_{\partial \Omega}\lambda_t \cdot \nabla u dS$$
To build intuition, if we now restrict ourself to the part of the boundary where $u \neq 0$ we can scale the Lagrange multiplier with $u$ (I think this can be motivated by shifting $u$ by a constant away from 0 but not sure about this one.)
$$-\int_{\partial \Omega} u \mathbf{n} \cdot \nabla u d S + \int_{\partial \Omega}u\frac{|\lambda_n|}{u} \mathbf{n} \cdot \nabla u dS+\int_{\partial \Omega}\lambda_t \cdot \nabla u dS$$
Which we can write as:
$$-\int_{\partial \Omega} u q_1 d S + \int_{\partial \Omega}u q_2 dS+\int_{\partial \Omega}\lambda_t \cdot \nabla u dS$$
So that the Lagrange contribution in the normal direction acts as an additional source term contribution on the boundary which naturally gets absorbed by the divergence theorem.
In conclusion, we are left with
$$\nabla \cdot \nabla u =0$$ with
$$\mathbf{n} \cdot \nabla u=q_1-q_2$$ on the boundary.
and
$$\delta(\lambda_t \cdot \nabla u)=0 $$
as an "orphan" constraint on the boundary with seemingly no influence on what is going on inside the domain.
In other words, the normal component of the gradient constraint on the boundary will manifest itself as an equivalent source term in a Neumann boundary condition.
Note that the shape of the boundary will determine how to choose the fields $\lambda_n$ and $\lambda_t$ in order to fulfil $\nabla u=0$.
For Helmholtz equation, it doesn't seem that setting $\nabla u =0$ would give only the trivial solution, unless it is covering the entire boundary. It rather seems as if it can be seen as a convoluted way of adding a Neumann condition source term that varies on the boundary. However, if $\nabla u =0$ holds on the entire boundary then you should get the trivial solution, since then this corresponds to $q_1=q_2=0$, or, rather, the two boundary integrals in the normal direction become identially zero + the third tangential term which doesn't affect the solution in the interior.
