I am assuming $\Omega$ is an open, connected subset of $\Bbb C$. More remarks pertinent to this assumption are made near the end of his answer.
With holomorphic
$f(x + iy) = u(x, y) + iv(x, y), \tag 1$
we have
$f^2(x + iy) = (u(x, y) + iv(x, y))^2 = u^2(x, y) - v^2(x, y) + 2i u(x, y) v(x, y) \tag 2$
is also holomorphic. Thus we may apply the Cauchy-Riemann equations to $ u^2(x, y) - v^2(x, y)$ and $2u(x, y) v(x, y)$ and see that
$(u^2(x, y) - v^2(x, y))_x = (2u(x, y) v(x, y))_y, \tag 3$
$(u^2(x, y) - v^2(x, y))_y = -(2u(x, y) v(x, y))_x. \tag 4$
If, then,
$u(x, y) v(x, y) = \beta, \tag 5$
a constant, we see from (3)-(4) that
$(u^2(x, y) - v^2(x, y))_x = 0, \tag 6$
$(u^2(x, y) - v^2(x, y))_y = 0, \tag 7$
i.e., that
$\nabla(u^2(x, y) - v^2(x, y)) = 0. \tag 8$
If now $(x_0, y_0)$ and $(x, y)$ are any two points in a given connected, open set $\Omega \subset \Bbb C$, we may join them by a differentiable path
$\gamma:[0, 1] \to \Omega, \tag 9$
such that
$\gamma(t) = (\gamma_x(t),\gamma_y(t)), \; \gamma(0) = (x_0, y_0), \; \gamma(1) = (x, y); \tag{10}$
then
$(u^2(x, y) - v^2(x, y)) - (u^2(x_0, y_0) - v^2(x_0, y_0)) = \displaystyle \int_0^1 \dfrac{d(u^2(\gamma_x(t), \gamma_y(t)) - v^2(\gamma_x(t), \gamma_y(t))}{dt} dt$
$= \displaystyle \int_0^1 \nabla(u^2(\gamma_x(t), \gamma_y(t) - v^2(\gamma_x(t), \gamma_y(t)) \cdot \gamma'(t) dt = \int_0^1 0 \cdot \gamma'(t) dt = 0;\, \tag{11}$
or
$u^2(x, y) - v^2(x, y) = u^2(x_0, y_0) - v^2(x_0, y_0) \tag{12}$
for any $(x, y)$ in the same component of $\Omega$ as $(x_0, y_0)$. It follows that $u^2(x, y) - v^2(x, y)$ is constant on components of $\Omega$, hence $f^2(x + iy) = u^2(x, y) - v^2(x, y) + 2iu(x, y) v(x, y)$ is also, hence $f(x + iy)$ as well.
We see upon scrutiny of the preceding argument that it requires $\Omega$ to be path connected to be effecitve; thus we have in actuality shown that $f(x + iy)$ is constant on the path components of $\Omega$; but these are the same as the components of $\Omega$ for open $\Omega \subset \Bbb R^2$. It is possiible for $f(x + iy)$ to take on different values on different components of $\Omega$.
As far as our OP Heuristics' proof, the only feature which bothers me is the division by $v'_x$, which may vanish on some subset of $\Omega$. For example, taking $f(x + iy) = (x + iy)^2 = (x^2 - y^2) + 2ixy$, we see that $v = 2xy$ hence $v'_x = 2y = 0$ on the entire real axis. I think there may be a way to work around this, but I'm not sure what it is at the moment.