w123jsc:
per your questions
1. Changing order of basis functions is critical in capturing accuracy, e.g. 0-th order basis function defines only 6 basis functions over a single tetrahedra. That means within that tet, there are six unknown vectors describing the EM field. That means you can only play with six coefficients to control the vector basis to model the EM field solution. With high-order basis, e.g. hierarchical curl-conforming, 1st order basis has around 36 (can't remember exactly) vector basis to describe the EM field, that means it can describe more pattern of the fields within that tet. This is why its more accurate. Of course, if you break your original tet with many smaller tets but stick to 0th order basis, you might get same phenomena. this is a well known topic called "hp-refinement" of FEM.
2. Typically, you want your element to be much smaller than the wavelength in your problem. e.g. your wavelength in your problem is 1m, a suggested value of max mesh size would be 0.1m, that means you have at least 10 elements describing the field within 1m. However, fine features and small gaps causes trouble, it create a lot of evanescent modes which need a lot of DoF to describe them accurately. In such case, especially when these evanescent modes plays a critical impact on results, you have to refine the mesh or increase basis # there at the cost of more computation time.
3. MoM face the same problem, I do not have data to compare against each other, but its similar for sure.