Something is puzzling me concerning the divergence theorem. Usually, one writes the divergence theorem as \begin{equation} \int_\mathcal{M} d^4x \sqrt{-g} \nabla_\mu v^\mu=\int_{\partial \mathcal{M}} d\Sigma_\mu v^\mu \tag{3.22} \end{equation} for some vector field $v^\mu$, where $d\Sigma^\mu$ is the directed surface element. See for example equation $(3.22)$ in A Relativist's Toolkit by Eric Poisson. I was considering turning the Maxwell Lagrangian into a boundary term, assuming that I am on-shell. In order to do that, we have \begin{equation} \begin{split} \int_\mathcal{M}d^4x \sqrt{-g} F^{\mu \nu}F_{\mu \nu}&=2\int_\mathcal{M}d^4x \sqrt{-g}F^{\mu \nu}\nabla_\mu A_\nu\\ &=2\int_\mathcal{M}d^4x\sqrt{-g}\nabla_\mu \left(F^{\mu \nu}A_\nu \right)\\ &-2\int_\mathcal{M}d^4x\sqrt{-g} \nabla_\mu F^{\mu \nu}A_{\nu} \, \end{split} \end{equation} where I just used symmetry of the Christoffel symbols and the anti-symmetry of the Faraday tensor. Now, using the equation of motion \begin{equation} \nabla_\mu F^{\mu \nu}=0 \end{equation} and the divergence theorem, I should be able to write \begin{equation}\label{puzzle} \int_\mathcal{M}d^4x \sqrt{-g} F^{\mu \nu}F_{\mu \nu}=2\int_{\partial \mathcal{M}} d\Sigma_\mu F^{\mu \nu} A_\nu\,. \end{equation} Now we can consider a magnetic solution given by \begin{equation} A=Q(1-\cos{\theta}) d\phi \end{equation} \begin{equation} F=dA=Q\sin{\theta} d\theta \wedge d\phi\,. \end{equation}
Given this solution, we have \begin{equation} F_{\mu \nu}F^{\mu \nu}=2\frac{Q^2}{r^4}\, \end{equation} assuming a metric of the form \begin{equation} ds^2=-f(r)dt^2+\frac{1}{f(r)}dr^2+r^2(d\theta^2 + \sin^2 \theta d\phi^2)\,. \end{equation} I'm actually dealing with Reissner-Nordstrom black holes with magnetic charge but I think my confusion holds more generically. My confusion is that if I take some boundary at constant $r$, or just some boundary such that its normal looks like \begin{equation} n_{\mu}=\alpha dt + \beta dr\,, \end{equation} then the right-hand side of \begin{equation} \int_\mathcal{M}d^4x \sqrt{-g} F^{\mu \nu}F_{\mu \nu}=2\int_{\partial \mathcal{M}} d\Sigma_\mu F^{\mu \nu} A_\nu\,. \end{equation} seems to yield zero trivially because $F$ has only angular components. However, I think it is very easy to come up with a region of spacetime $\mathcal{M}$ such that the left-hand side does not vanish. What did I do wrong in this reasoning? I thought it might have to do with the smoothness of the vector field $v^{\mu}=F^{\mu \nu}A_{\nu}$ in which we are applying the Stokes' theorem, but since I saw the divergence theorem stated above without any assumptions on the vector field, I am not sure if this is the problem. The vector field yields \begin{equation} v=\frac{Q^2}{r^4 \sin{\theta}}(1-\cos{\theta})\frac{\partial}{\partial \theta}\,, \end{equation} It is possible to see that $v$ has a singularity at $\theta=\pi$. This could be a coordinate artifact but it is easy to see that $g_{\mu \nu}v^\mu v^\nu$ also has this singular point. For this reason, I think it might be problematic to apply the divergence theorem to $v$. Am I getting the problem right?
EDIT1: I am considering a region $\mathcal{M}$ of spacetime that does not contain the origin. Since I am thinking of this in the context of a black hole, I don't want to include the singularity $r=0$ in the patch over which I am integrating. As an example of this procedure, see Eq. $(5.15)$ of https://arxiv.org/abs/1606.08307 where an electric solution is assumed. Now if we assume an electric solution, we have \begin{equation} A=\left(\frac{Q}{r_{+}}-\frac{Q}{r} \right)dt \end{equation} \begin{equation} F=dA=\frac{Q^2}{r^2}dr \wedge dt \end{equation} where $r_{+}$ is the event horizon. Hence \begin{equation} F_{\mu \nu}F^{\mu \nu}=-2\frac{Q^2}{r^4}\,, \end{equation} which means the bulk integration on $\mathcal{M}$ just changes by a sign with respect to the magnetic solution. This is related to electromagnetic duality. However, in this case using the Stokes theorem leads to a non-zero boundary contribution, as it should, I think. Now my puzzle is why in one case we can use Stokes and the other it seems to fail. It turns out that in this case the vector field in which we apply Stokes is \begin{equation} v=\frac{Q}{r^2}\left(\frac{1}{r}-\frac{1}{r_{+}} \right) \frac{\partial}{\partial r} \end{equation} which seems to be smooth everywhere if we look at $g(v,v)$, using a good coordinate system.
EDIT2: According to Introduction to Smooth Manifolds by John M. Lee, the vector field must be smooth. But if that is the case, when we use the variational principle to derive Maxwell's equations, we assume that $F^{\mu \nu}\delta A_\nu$ is smooth to make it into a boundary term that vanishes. Of course here we are off-shell, but it still seems strange to me that $F^{\mu \nu}A_\mu$ turns out to not be smooth for a purely magnetic solution. Do you have any insights?