To answer my own question; this is my conclusion after reading and investigating.
According to Wikipedia a harmonic conjugate can be defined as follows:
A real-valued function $u(x,y)$ defined on a connected open set
$\Omega \subset \mathbb R^2$ is said to have a conjugate (function)
$v(x,y)$ if and only if they are respectively the real and imaginary
parts of a holomorphic function $f(z)$ of the complex variable
$z:=x+iy \in \Omega$
That is, $v$ is conjugate to $u$ if $f(z):=u(x,y)+iv(x,y)$ is
holomorphic on $\Omega$.
A holomorphic function is defined by Wikipedia as:
A complex-valued function of one or more complex variables that is, at
every point of its domain, complex differentiable in a neighborhood of
the point.
It also says that "all holomorphic functions are complex analytic functions, and vice versa, is a major theorem in complex analysis".
So for $v(x,y)$ to be a conjugate to $u(x,y)$ the function $f(z)=u(x,y)+iv(x,y)$ must holomorphic or analytic.
According to Theorem 5.2 in this PDF from MIT the following is true:
$f(z)=u(x,y)+iv(x,y)$ is analytic on a region $A$ then both $u$ and
$v$ are harmonic functions on $A$.
So in conclusion; since a harmonic conjugate is defined as $v(x,y)$ being the imaginary part of a holomorphic function $f(z)$ and a holomorphic function must have harmonic parts $u$ and $v$, then $u(x,y)$ must be harmonic for being able to find a harmonic conjugate.
All this is perhaps obvious since it's all in the definitions, but I wanted to post an answer regardless so it becomes obvious.
Best Answer
To the function to be harmonic its Laplacian should be zero $$ \Delta f = f_{xx}(x,y)+f_{yy}(x,y) = 0 $$ Just check it $$ f_{xx} = \left(\frac {x}{x^2+y^2}\right)_{xx} = \left [\frac {-x^2+y^2}{\left( x^2+y^2\right)^2}\right]_x = \frac {2x(x^2-3y^2)}{\left( x^2+y^2\right)^3} \\ f_{yy} = \left(\frac {x}{x^2+y^2}\right)_{yy} = \left [-\frac {2xy}{\left ( x^2+y^2\right)^3} \right ]_y = -\frac {2x(x^2-3y^2)}{\left( x^2+y^2\right)^3} $$ which means $\Delta f = 0$ hence harmonic.