Let \(M\), \(N\) be smooth compact Riemannian manifolds without boundary, \(m= \dim M\) and the sectional curvature of \(N\) is nonpositive, and \(\phi:M\to N\) be a smooth map. This paper deals with the following conjecture: Any weakly harmonic map of finite energy from \(M\) to \(N\) is smooth provided that there are no harmonic spheres \(S^l\) in \(N\), for \(2\leq l\leq m-1\). Consider the penalized energy: \(I_\varepsilon(u)= \int_M\left({1\over 2} |Du|^2+{F(u)\over \varepsilon^2}\right) dx\), where \(F\) is a smooth function of \(u\) such that \(F(p)= \text{dist}^2(p, N)\) if \(\text{dist}(p, N)\leq \delta\); or \(=4\delta^2\) if \(\text{dist}(p, N)\geq 2\delta\). Here \(N\) is viewed as a submanifold of \(R^k\) and \(\delta\) is chosen so that \(\text{dist}^2(p, N)\) is smooth for \(p\in \{p:\text{dist}(P, N)\leq 2\delta\}\). In this paper, the authors use the gradient flow of \(I_\varepsilon(\cdot)\) to derive theorems similar to \textit{R. Schoen} and \textit{K. Uhlenbeck} [J. Differ. Geom. 17, 307-335 (1982; Zbl 0521.58021)] and \textit{F. H. Lin} [C. R. Acad. Sci., Paris, Sér. I 323, No. 9, 1005-1008 (1996; Zbl 0871.58029)], and thus recover the \textit{J. Eells} and \textit{J. H. Sampson's} theorem [Am. J. Math. 86, 109-160 (1964; Zbl 0122.40102)]. The main results of this paper are as follows: First, consider solutions of \((*)\) \(\Delta u_\varepsilon+{1\over \varepsilon^2}f(u_\varepsilon)= 0\) in \(B_1\). (A) Let \(\varepsilon_i\downarrow 0\), and \(u_{\varepsilon_i}\) be a sequence of solutions of \((*)\) with \(I_{\varepsilon_i}(u_{\varepsilon_i})\leq K<\infty\), and \(u_{\varepsilon_i}\to u\) weakly in \(H^1(B_1)\). Suppose that there is no harmonic \(S^2\) in \(N\). Then \(e(u_{\varepsilon_i})\equiv \left({1\over 2}|Du_{\varepsilon_i}|^2+ {1\over \varepsilon^2_i} F(u_{\varepsilon_i})\right) dx\rightharpoonup{1\over 2}|Du|^2 dx\), as Radon measures. In particular, \(u_{\varepsilon_i}\to u\) strongly in \(H^1_{\text{loc}}(B_1)\), and \(\int_{B_1} {1\over \varepsilon^2_i} F(u_{\varepsilon_i}) dx\to 0\). (B) Under the assumption that there is no harmonic \(S^2\) in \(N\), the map \(u\) obtained in (A) is a stationary harmonic map. In particular, the singular set of \(u\) has Hausdorff dimension at most \(m-4\). If, in addition, \(N\) has no harmonic \(S^l\) for \(3\leq l\leq m-1\), then \(u\) is smooth and \(u_{\varepsilon_i}\to u\) in \(C^k\) norm, for any \(k\geq 1\). Next, consider solutions of \((**)\) \(\partial_t u_\varepsilon-\Delta u_\varepsilon-{1\over \varepsilon^2} f(u_\varepsilon)= 0\), in \(B_1\times (0,1)\). (C) Let \(\varepsilon_i\downarrow 0\), \(u_{\varepsilon_i}\) be a sequence of solutions of \((**)\) with \(\int_{B_1\times [0,1]}(|\partial_t u_{\varepsilon_i}|^2+ e(u_{\varepsilon_i})) dx dt\leq K< \infty\). Suppose that there is no harmonic \(S^2\) in \(N\), and \(u_{\varepsilon_i}\to u\) weakly in \(H^1(B_1\times (0,1))\). Then \(e(u_{\varepsilon_i}) dx dt\rightharpoonup{1\over 2}|Du|^2 dx dt\), as Radon measures. In particular, \(u_{\varepsilon_i}\to u\) strong in \(H^1_{\text{loc}}(B_1\times (0,1))\). The limit map \(u\) is a weak solution of the equation \(\partial_t u= \tau(u)\) in \(M\times R_+\), where \(\tau(u)\) is the stress tension-field of \(u\), with \({\mathcal P}^m(\text{sing}(u))= 0\), and \(u\) satisfies both, energy inequality and monotonicity inequality. Here \({\mathcal P}^m\) denotes the \(m\)-dimensional Hausdorff measure with respect to the parabolic metric in \(R^{m+1}\). Moreover, some connections between critical points of \(I_\varepsilon(\cdot)\) and weakly harmonic maps from \(M\) into \(N\) are established. Finally, the authors propose the following questions: For a compact, smooth Riemannian manifold \(N\), are there any quasi-harmonic \(S^l\), \(l\geq 3\), of finite energy? Here a map \(\phi: R^l\to N\) is called to be quasi-harmonic if \(\phi\) is a nonconstant, smooth map from \(R^l\) to \(N\) such that it is a critical point of \(\int_{R^l}|Du|^2 e^{-{|y|^2\over 4}} dy\).