Hi, harmonicity in 2d is preserved under mappings that satisfy the Cauchy-Riemann equations.

What about 3D? What conditions should a mapping satisfy to preserve harmonicity?

is there a general characterization a la CR for 3D?

Here is an example of non-trivial such mapping

Let $u(x,y,z)=U(X,Y,Z)$ where $$X=xy+z,~~~~ Y= \frac{\sqrt3}4 (x^2-y^2) - \frac{xy}{2}+z ,~~~~ Z= -\frac{\sqrt3}4 (x^2-y^2) - \frac{xy}{2}+z$$

(I found this example by first assuming $X=xy+z$ then guessing for Y,Z from the overdetermined system that they satisfy... hope it's right...)