[WIP] Darwin: solve -lap(curl(curl(Z))) instead of bilaplacian(Z)

[RemiLehe]

The semi-implicit Darwin solver writes the vector-potential increment as dA = curl(Z) and solves the resulting fourth-order equation for Z. That equation was discretized with a bi-Laplacian, i.e.

bilaplacian(Z) + curl(chi curl(Z)) = 2 laplacian(B) + 2 mu_0 curl(J)

but bilaplacian(Z) is only equal to the fourth-order term -laplacian(curl(curl(Z))) where div(Z) = 0, which nothing in the scheme imposes on Z itself (only div(dA) = div(curl(Z)) = 0 is guaranteed, by construction). Evaluate the curl-curl form instead, which is what the derivation actually produces:

-laplacian(curl(curl(Z))) + curl(chi curl(Z)) = 2 laplacian(B) + 2 mu_0 curl(J)

The right-hand side is unchanged. The operator reuses the intermediate curl(Z) it was already computing for the plasma-response term, so the evaluation is one extra curl plus one vector Laplacian, and the scratch space it needs narrows from two ghost cells to one (every stencil in the chain is now nearest-neighbour, with guard cells refilled in between).

Since the bi-Laplacian is no longer used anywhere, remove FiniteDifferenceSolver::ComputeVectorBiLaplacian along with its Cartesian/Cylindrical/Spherical variants, and the fourth-derivative stencils Dxxxx/Dyyyy/Dzzzz/Dxxyy/Dyyzz/Dxxzz that only it used.

Both operators agree exactly on discretely divergence-free fields (the single-pass biharmonic stencil is the square of the 3-point Laplacian stencil, and curl curl = grad div - lap holds discretely on the Yee mesh), so results are unchanged where the iterates stay solenoidal: the 1D and 2D Darwin tests reproduce the previous GMRES history iterate for iterate, up to a 1-ULP difference in one particle-momentum sum.