Solution of electromagnetic wave in frequency domain seems strange

Question

Bug introduced in 10, fixed in 11.1.

---

I want to solve the electromagnetic wave equation in frequency domain. A known solution, a plane wave Exp[ I k0 x] is used to set the Dirichlet conditions to all the boundaries. I expect to get the plane wave solution. The following code runs well ,but the result is strange. I don't know why. It seems high-frequency oscillation occurs. Is there some options to eliminate the unusual oscillation? thanks a lot!

<< NDSolve`FEM`  

λ = 0.53; k0 = 2 π/λ; R = λ;

mesh = ToElementMesh[FullRegion[2], {{0, R}, {0, R}}, 
   "MaxCellMeasure" -> 0.0005];

mesh["Wireframe"]

op = 
  Most[Curl[Curl[{u[x, y], v[x, y], 0}, {x, y, z}], {x, y, z}] - 
   k0^2 {u[x, y], v[x, y], 0}]

pde = op == {0, 0};

Subscript[Γ, D] = 
  DirichletCondition[{u[x, y] == 0., v[x, y] == Exp[I k0 x]}, True];

{us, vs} = 
  NDSolveValue[{pde, Subscript[Γ, D]}, {u, v}, {x, y} ∈ mesh]

DensityPlot[Re[vs[x, y]], {x, y} ∈ mesh, 
  ColorFunction -> "Rainbow", 
  PlotLegends -> Automatic, 
  PlotPoints -> 50, 
  PlotRange -> All]

Answer 1

Update

OP's code produces the correct result since v11.1, I think it's safe to say the "noisy" result is caused by a bug.

---

Here's a solution based on finite difference method (FDM). The definition of pdetoae (a general purpose function discretizing PDE to algebraic equations with finite difference formula) can be found here.

λ = 53/100; k0 = 2 π/λ; R = λ;

op = Most[Curl[Curl[{u[x, y], v[x, y], 0}, {x, y, z}], {x, y, z}] - 
    k0^2 {u[x, y], v[x, y], 0}];

pde = op == {0, 0};
bc = Function[{x, y}, {u[x, y] == 0, v[x, y] == Exp[I k0 x]}] @@@ {{0, y}, {R, y}, {x, 
      0}, {x, R}} // Flatten;

domain = {0, R};
points = 50;
grid = Array[# &, points, domain];
difforder = 4;
(*Definition of pdetoae isn't included in this code piece,
  please find it in the link above. *)
ptoa = pdetoae[{u, v}[x, y], {grid, grid}, difforder];
del = Most@Rest@# &;

ae = del /@ del@# & /@ ptoa@pde;
aebc = MapAt[del, ptoa@bc, List /@ Range@4];
{b, mat} = CoefficientArrays[{ae, aebc} // Flatten, 
   Outer[#[#2, #3] &, {u, v}, grid, grid] // Flatten];
sollst = LinearSolve[-mat, N@b];
solmat = ArrayReshape[sollst, {2, points, points}];
{solu, solv} = ListInterpolation[#, {grid, grid}] & /@ solmat;

Outer[Plot3D[#@First[#2][x, y], {x, 0, R}, {y, 0, R}, PlotLabel -> #@Last@#2] &, {Re, 
   Im}, {{solu, "u"}, {solv, "v"}}, 1] // Grid

Remark

1. Number of grid points i.e. points should be even, or the error of $\text{Im}(v)$ will be quite large, I'm not sure about the reason.

2. Difference order i.e. difforder should be even and greater than 2, or the error will be large and seems unlikely to be alleviated by increasing points. I guess this might explain why OP's code doesn't work well, because according to this document page, the default "MeshOrder" of ToElementMesh is exactly 2, but since I'm still in v9 and the Wolfram Cloud is too slow, I'd like to stop here.