Let $S$ be a punctured Riemann surface with Euler characteristic $χ(S)<0$. For any unitary representation $ρ: π_1(S) \to U(N)$, we introduce its renormalized energy and its harmonic representatives, which are equivariant harmonic maps from the universal cover of $S$ to the unit sphere in $\mathbb{C}^N$. Our main result is that if a sequence of unitary representations $ρ_j$ strongly converges, then their renormalized energies converge to $\fracπ{4}|χ(S)|$ and the shape of their harmonic representatives converges to a unique rescaled hyperbolic metric. Combining this statement with examples of strongly converging representations provided by random matrix theory, we derive the following applications. (1) If $π_1(S)$ is a free group, then for a random $ρ: π_1(S) \to U(N)$, the shape of its harmonic representatives concentrates around a rescaled hyperbolic metric with high probability as $N\to \infty$. (2) For any closed hyperbolic surface, a finite covering admits a harmonic immersion into some Euclidean unit sphere, which is almost isometric after rescaling. (3) There are closed, branched, minimal surfaces $\mathfrak{S}_j$ in some Euclidean unit spheres such that $\mathfrak{S}_j$ Benjamini-Schramm converges to a rescaled hyperbolic plane as $j\to \infty$, and the Gaussian curvature $K_j$ of $\mathfrak{S}_j$ satisfies $\lim_{j\to \infty} \frac{1}{\mathrm{Area}(\mathfrak{S}_j)}\int_{\mathfrak{S}_j} |K_j+8|=0.$