[proofplan]
We work in real coordinates, writing $z=s+it$ on $U$ and $f(z)=p(s,t)+iq(s,t)$. The second-order chain rule expresses $\Delta(u\circ f)$ in terms of the first and second derivatives of $p,q$ and the Hessian of $u$. Holomorphicity gives the [Cauchy-Riemann equations](/page/Cauchy-Riemann%20Equations), which make the mixed Hessian terms cancel and make the coefficients of $u_{xx}$ and $u_{yy}$ equal; harmonicity of $p$ and $q$ removes the remaining first-derivative terms. The result is the conformal Laplacian transformation formula $\Delta(u\circ f)=|f'|^2(\Delta u)\circ f$, and the right-hand side vanishes because $u$ is harmonic on $V$.
[/proofplan]
[step:Pass to real coordinates for the conformal map and the harmonic function]
Define the real-coordinate domains
\begin{align*}
\widetilde U := \{(s,t)\in \mathbb R^2 : s+it \in U\}
\end{align*}
and
\begin{align*}
\widetilde V := \{(x,y)\in \mathbb R^2 : x+iy \in V\}.
\end{align*}
Since $U$ and $V$ are open in $\mathbb C$ and the identification $(a,b)\mapsto a+ib$ is a homeomorphism $\mathbb R^2\to \mathbb C$, the sets $\widetilde U$ and $\widetilde V$ are open in $\mathbb R^2$.
Define the real-coordinate representative of $f$ by
\begin{align*}
F:\widetilde U &\to \widetilde V
\end{align*}
where
\begin{align*}
F(s,t) := (p(s,t),q(s,t))
\end{align*}
is determined by
\begin{align*}
f(s+it)=p(s,t)+iq(s,t).
\end{align*}
Because $f$ is conformal on $U$, it is holomorphic on $U$; hence $p,q\in C^2(\widetilde U;\mathbb R)$ and the Cauchy-Riemann equations hold:
\begin{align*}
p_s=q_t
\end{align*}
and
\begin{align*}
p_t=-q_s
\end{align*}
on $\widetilde U$. Also, the real and imaginary parts of a [holomorphic function](/page/Holomorphic%20Function) are harmonic, so
\begin{align*}
p_{ss}+p_{tt}=0
\end{align*}
and
\begin{align*}
q_{ss}+q_{tt}=0
\end{align*}
on $\widetilde U$.
Define the real-coordinate representative of $u$ by
\begin{align*}
\widetilde u:\widetilde V &\to \mathbb R
\end{align*}
where
\begin{align*}
\widetilde u(x,y):=u(x+iy).
\end{align*}
Then $\widetilde u\in C^2(\widetilde V;\mathbb R)$, and harmonicity of $u$ means
\begin{align*}
\widetilde u_{xx}+\widetilde u_{yy}=0
\end{align*}
on $\widetilde V$.
[/step]
[step:Compute the Laplacian of the pullback by the second-order chain rule]
Let
\begin{align*}
v:\widetilde U &\to \mathbb R
\end{align*}
be the real-coordinate representative of $u\circ f$, defined by
\begin{align*}
v(s,t):=\widetilde u(F(s,t))=\widetilde u(p(s,t),q(s,t)).
\end{align*}
Since $\widetilde u\in C^2(\widetilde V;\mathbb R)$ and $p,q\in C^2(\widetilde U;\mathbb R)$, the composition rule gives $v\in C^2(\widetilde U;\mathbb R)$.
Fix $(s,t)\in \widetilde U$, and evaluate all derivatives of $\widetilde u$ at $(p(s,t),q(s,t))$. The second-order chain rule gives
\begin{align*}
v_{ss}=\widetilde u_{xx}p_s^2+2\widetilde u_{xy}p_s q_s+\widetilde u_{yy}q_s^2+\widetilde u_x p_{ss}+\widetilde u_y q_{ss}.
\end{align*}
The same rule in the $t$-variable gives
\begin{align*}
v_{tt}=\widetilde u_{xx}p_t^2+2\widetilde u_{xy}p_t q_t+\widetilde u_{yy}q_t^2+\widetilde u_x p_{tt}+\widetilde u_y q_{tt}.
\end{align*}
Adding these two identities yields
\begin{align*}
\Delta v=(p_s^2+p_t^2)\widetilde u_{xx}+2(p_s q_s+p_t q_t)\widetilde u_{xy}+(q_s^2+q_t^2)\widetilde u_{yy}+\widetilde u_x(p_{ss}+p_{tt})+\widetilde u_y(q_{ss}+q_{tt}).
\end{align*}
[guided]
The goal is to compute $\Delta(u\circ f)$, but $u$ is naturally a function of the target coordinates while $u\circ f$ is a function of the source coordinates. That is why we introduced
\begin{align*}
v:\widetilde U &\to \mathbb R
\end{align*}
by
\begin{align*}
v(s,t):=\widetilde u(p(s,t),q(s,t)).
\end{align*}
Here $(s,t)$ are source variables, while $(x,y)$ are target variables. Since $\widetilde u\in C^2(\widetilde V;\mathbb R)$ and $p,q\in C^2(\widetilde U;\mathbb R)$, the ordinary [multivariable chain rule](/theorems/830) applies twice, so $v\in C^2(\widetilde U;\mathbb R)$.
Fix a point $(s,t)\in \widetilde U$. In the following formulas, the derivatives $\widetilde u_x,\widetilde u_y,\widetilde u_{xx},\widetilde u_{xy},\widetilde u_{yy}$ are evaluated at the target point $(p(s,t),q(s,t))$. First differentiating in the $s$-direction gives
\begin{align*}
v_s=\widetilde u_x p_s+\widetilde u_y q_s.
\end{align*}
Differentiating this identity once more in $s$, using the product rule and the chain rule on $\widetilde u_x(p,q)$ and $\widetilde u_y(p,q)$, gives
\begin{align*}
v_{ss}=\widetilde u_{xx}p_s^2+\widetilde u_{xy}q_s p_s+\widetilde u_x p_{ss}+\widetilde u_{yx}p_s q_s+\widetilde u_{yy}q_s^2+\widetilde u_y q_{ss}.
\end{align*}
Because $\widetilde u\in C^2$, Clairaut's theorem gives $\widetilde u_{xy}=\widetilde u_{yx}$, so this becomes
\begin{align*}
v_{ss}=\widetilde u_{xx}p_s^2+2\widetilde u_{xy}p_s q_s+\widetilde u_{yy}q_s^2+\widetilde u_x p_{ss}+\widetilde u_y q_{ss}.
\end{align*}
The same computation in the $t$-direction gives
\begin{align*}
v_{tt}=\widetilde u_{xx}p_t^2+2\widetilde u_{xy}p_t q_t+\widetilde u_{yy}q_t^2+\widetilde u_x p_{tt}+\widetilde u_y q_{tt}.
\end{align*}
The Laplacian in the source variables is $\Delta v=v_{ss}+v_{tt}$. Adding the two displayed identities and grouping the coefficients of each derivative of $\widetilde u$ gives
\begin{align*}
\Delta v=(p_s^2+p_t^2)\widetilde u_{xx}+2(p_s q_s+p_t q_t)\widetilde u_{xy}+(q_s^2+q_t^2)\widetilde u_{yy}+\widetilde u_x(p_{ss}+p_{tt})+\widetilde u_y(q_{ss}+q_{tt}).
\end{align*}
This formula is the exact point where conformality will enter: the Cauchy-Riemann equations will simplify the three Hessian coefficients, and harmonicity of $p$ and $q$ will eliminate the first-derivative terms.
[/guided]
[/step]
[step:Use the Cauchy-Riemann equations to reduce the Laplacian formula]
The Cauchy-Riemann equations $p_s=q_t$ and $p_t=-q_s$ imply
\begin{align*}
p_s q_s+p_t q_t=p_s q_s+(-q_s)q_t=q_s(p_s-q_t)=0.
\end{align*}
They also imply
\begin{align*}
q_s^2+q_t^2=(-p_t)^2+p_s^2=p_s^2+p_t^2.
\end{align*}
Since $p$ and $q$ are harmonic on $\widetilde U$, we have
\begin{align*}
p_{ss}+p_{tt}=0
\end{align*}
and
\begin{align*}
q_{ss}+q_{tt}=0.
\end{align*}
Substituting these identities into the formula for $\Delta v$ gives
\begin{align*}
\Delta v=(p_s^2+p_t^2)(\widetilde u_{xx}+\widetilde u_{yy}).
\end{align*}
Equivalently, since $f'(s+it)=p_s(s,t)+iq_s(s,t)$ and $q_s=-p_t$, the coefficient satisfies
\begin{align*}
p_s^2+p_t^2=|f'(s+it)|^2.
\end{align*}
Thus
\begin{align*}
\Delta v(s,t)=|f'(s+it)|^2(\Delta \widetilde u)(F(s,t)).
\end{align*}
[/step]
[step:Apply harmonicity of $u$ and return to complex notation]
Because $u$ is harmonic on $V$, its real-coordinate representative $\widetilde u$ satisfies
\begin{align*}
\Delta \widetilde u(x,y)=0
\end{align*}
for every $(x,y)\in \widetilde V$. Since $F(\widetilde U)\subset \widetilde V$, for every $(s,t)\in \widetilde U$ we obtain
\begin{align*}
\Delta v(s,t)=|f'(s+it)|^2(\Delta \widetilde u)(F(s,t))=0.
\end{align*}
Therefore $v$ is harmonic on $\widetilde U$. Translating back through the identification $(s,t)\mapsto s+it$, this says precisely that $u\circ f\in C^2(U;\mathbb R)$ and $\Delta(u\circ f)=0$ on $U$. Hence $u\circ f$ is harmonic on $U$.
[/step]