Different results in DEigensystem compared to NDEigensystem for Laplacian eigenvalue problem (-Δu=λu) on unit square

Question

I want to calculate the solution to the Laplacian eigenvalue problem on the unit square with trivial Dirichlet boundary conditions: $$- \Delta u(x,y) = \lambda u(x,y) \text{ on } {[0,1]}^2$$ with $u(0,y)=0$,$u(1,y)=0$,$u(x,0)=0$,$u(x,1)=0$.

However, Mathematica 12 reports different eigenfunctions when using NDEigensystem in contrast to DEigensystem using the following codes:

DEigensystem version:

{vals, funs} = 
DEigensystem[{-Laplacian[u[x, y], {x, y}], 
DirichletCondition[u[x, y] == 0, True]}, 
u[x, y], {x, y} ∈ Rectangle[], 2];
Table[ContourPlot[funs[[i]], {x, y} ∈ Rectangle[], 
PlotRange -> All, PlotLabel -> vals[[i]], PlotTheme -> "Minimal", 
Axes -> True], {i, Length[vals]}]

NDEigensystem version:

{vals, funs} = 
NDEigensystem[{-Laplacian[u[x, y], {x, y}], 
DirichletCondition[u[x, y] == 0, True]}, 
u[x, y], {x, y} ∈ Rectangle[], 2, 
Method -> {"PDEDiscretization" -> {"FiniteElement", 
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];
Table[ContourPlot[funs[[i]], {x, y} ∈ Rectangle[], 
PlotRange -> All, PlotLabel -> vals[[i]], PlotTheme -> "Minimal", 
Axes -> True], {i, Length[vals]}]

For the second eigenfunction, the DEigensystem reports the classical textbook eigenfunction, while the numerical solution with NDEigensystem is fundamentally different, although the mesh discretization is set to a very small value.

Why is that?

Answer 1

As already pointed out in the comments by @march and @xzczd, the second lowest state with eigenvalue $\lambda_{1,2} = \lambda_{2,1} = 5 \pi^2$ is doubly degenerate.

DEigensystem

NDEigensystem

This means that the corresponding eigenfunctions are not only determined up to a scaling (as for the lowest state). They are rather determined to be some orthogonal basis of the eigenspace

$$ E_{5 \pi^2} = \{a \phi_{1,2} + b \phi_{2,1} \mid a,b \in \mathbb{R}, \, -\Delta \phi_{1,2} = 5 \pi^2 \phi_{1,2}, \, -\Delta \phi_{2,1} = 5 \pi^2 \phi_{2,1}, \, \phi_{1,2} \perp \phi_{2,1}\} $$

We have $\text{dim}(E_{5 \pi^2}) = 2$. The results from NDEigensystem ($\phi_{1,2,\text{ND}}, \phi_{2,1,\text{ND}}$) are therefore also valid solutions because they span the same eigenspace:

$$ E_{5 \pi^2} = \text{span}\{\phi_{1,2,\text{NDEigen}}, \phi_{2,1,\text{NDEigen}}\} \\ = \text{span}\{\phi_{1,2,\text{DEigen}}, \phi_{2,1,\text{DEigen}}\} \\ = \text{span}\{\sin(1 \pi x)\sin(2 \pi y), \sin(2 \pi x)\sin(1 \pi y)\} $$