[This is an answer made of two comments and an example]

Since the laplacian is elliptic with real-analytic coefficients, a harmonic function $f$ is real-analytic in its domain of definition. Hence the set $C$ of critical points of $f$ is a real-analytic subset of $R^3$, and as such it admits a locally finite partition into real-analytic locally closed smooth submanifolds. Thus if $\dim C≤1$, it is locally a finite union of analytic open arcs and singular points (but the curves might not extend smoothly across those points).

A reference on real analytic functions (reedited in 2002) might be

S. Krantz, H. Parks, A primer of real analytic functions. Birkhäuser Verlag, 1992.

~~But maybe the "curve selection lemma" in Milnor's "Singular points on complex hypersurfaces" would be enough ~~. edit: it concerns real *algebraic* subsets.

As an example of a curve of critical points not extending throuh a singular point, take the harmonic polynomial
$f(x,y,z)=y^3-3x^2y+y^3z-yz^3$, which has critical locus $y=0,z^3=-3x^2$. But of course you have the singular parametrization $x=3t^3, z=-3t^2$. I don't know if they exist in general.

Addendum : in fact the critical locus of a harmonic can polynomial can have an arbitrary (real) plane algebraic curve as a union of irreducible components. Let $P(x,y)$ be a real two variable polynomial,
and define $$f(x,y,z)=\Sigma_k \frac{z^{2k+1}}{(2k+1)!}(-\Delta_{x,y})^k P(x,y) \; .$$
It is easy to check that $f$ is harmonic, and $df$ vanishes on $z=0$, $P(x,y)=0$.