I would like to see either a proof or a reference for the following assertions:
Let $(M^3,g)$ be a closed, connected and oriented Riemannian $3$-manifold. For each homotopy class $\alpha \in [M, S^1]$, there exists a minimizer $u : M \to S^1$ of the energy $E(v) = \int_M \vert dv \vert^2$ restricted to $\alpha$. Moreover, if $h = u^*(d \theta)$, then $h$ is a harmonic $1$-form on $M$.
Let $(M^3, g)$ be a compact, connected and oriented Riemannian $3$-manifold with boundary. For each homotopy class $\alpha \in [M, S^1]$, there exists a minimizer $u : M \to S^1$ of the energy $E(v) = \int_M \vert dv \vert^2$ restricted to $\alpha$. Moreover, if $h = u^*(d \theta)$, then $h$ is a harmonic $1$-form on $M$ and $h$ satisfies Neumann boundary condition along $\partial M$.
I'm no expert in Hodge theory, but it seems the existence can be reduced to the Hodge isomorphism theorem, and the Harmonic condition can be derived from a direct computation.
Preliminaries
Given a metric and orientation, $S^1$ can be identified locally with $\mathbb{R}$ up to an additive constant, and so $TS^1$ has a canonical trivialization $TS^1\cong S^1\times\mathbb{R}$, and with this trivialization we may treat $du$ as a one-form (equicalent to $h$). $du$ is closed (since closure is a local property), but not necessarily exact, despite this abuse of notation. Also, we can still perform addition modulo $2\pi$, $+:S^1\times\mathbb{R}\to S^1$ independent of the local identification.
Existence
The existence of a minimum can at least be reduced to Hodge isomorphism theorem (for manifolds with boundary), which states that each De Rham cohomology class $\gamma\in H^k_{dR}(M;\mathbb{R})$ has a unique representative, which minimizes the $L^2$ norm on its class, and this representative is harmonic and satisfies the Neumann boundary conditions.
Choose arbitrary base points $p\in M,e\in S^1$ and let $\alpha'\in[(M,p),(S^1,e)]$ be a (smooth) homotopy class of pointed maps. Let $\beta=\{du:u\in\alpha'\}$ (interpreted as one forms). We wish to prove that $\beta$ is a (De Rham) cohomology class. To do this, let $u,v\in\alpha'$, and let $F:[0,1]\times M\to S^1$ be a smooth homotopy between them. Let $f:M\to\mathbb{R}$ be defined by $$ f(x)=\int_0^1 \frac{\partial F}{\partial t}(t,x)dt $$ identifying the derivative $\partial F/\partial t$ as a map $[0,1]\times M\to\mathbb{R}$ by the trivialization of $TS^1$. Since $v=u+f$, we have $dv=dv+df$ and thus $du$ and $dv$ are cohomologous.
Conversely If $\omega$ is cohomologous to $du$ for some $u\in\alpha'$, there is a unique $f\in C^\infty M$ such that $f(p)=0$ and $\omega=du+df$. This implies $\omega=d(u+f)$ where $u+f\in\alpha'$
The above gives a bijection $d:\alpha'\to\beta$ between the smooth homotopy class of pointed maps and a De Rham cohomology class of one forms. Under this identification, the energy is equivalent to the $L^2$ norm. Therefore, there is by Hodge theorem a unique one form, and thus a unique pointed map, which minimizes this energy.
Harmonicity
Assume a minimizing $u$ exists and is smooth. Rather than relying on Hodge theorem, we can consider the energy functional to be a Lagrangian $\mathcal{L}$ and write the corresponding Euler-Lagrange equation, which gives the desired results explicitly.
Let $u_\lambda:M\to S_1$ with $\lambda\in(-\epsilon,\epsilon)$ be a one parameter family of maps. We may write $\delta=\left(\frac{\partial}{\partial\lambda}\right)_{\lambda=0}$. Note that we may treat $\delta u$ as a function $M\to\mathbb{R}$, and that $\delta$ commutes with the differential/exterior derivative $d$ and the hodge star $\star$, and distributes over wedge products (evident again from working locally and identifying $S^1$ with $\mathbb{R}$). We may compute the variation of the Lagrangian directly. $$ \delta\mathcal{L}[u_\lambda]=\delta\int_Mdu_\lambda\wedge\star du_\lambda $$ $$ \delta\mathcal{L}[u_\lambda]=\int_M\left(d\delta u_\lambda\wedge\star du_\lambda+du_\lambda\wedge\star d\delta u\right) $$ $$ \delta\mathcal{L}[u_\lambda]=\int_M2d\delta u_\lambda\wedge\star du_\lambda $$ $$ \delta\mathcal{L}[u_\lambda]=\int_M2\left[d\left(\delta u_\lambda\wedge\star du_\lambda\right)-\delta u_\lambda\wedge d\star du_\lambda\right] $$ $$ \delta\mathcal{L}[u_\lambda]=2\int_{\partial M}\delta u_\lambda\wedge\star du_\lambda-2\int_M\delta u_\lambda\wedge d\star du_\lambda $$ Here, assuming $u_0$ minimizes $\mathcal{L}$, we can apply the fundamental lemma of variational calculus, obtaining a differential equation and boundary condition for $u=u_0$. $$ d\star du=0,\ \ \ \ \star du|_{\partial M}=0 $$ This implies that the codifferential $d^*du=0$, and since we already have $ddu=0$ we arrive at $\Delta du=0$. The boundary condition implies that $\star du\wedge n=0$, where $n$ is a unit normal covector field. By the definition of the Hodge star, this gives $\langle du,n\rangle=0$ on $\partial M$.