Puzzle concerning the Divergence Theorem
It's not due to the behavior at just the origin, it's due to the ill-defined behavior of $A$ along the line $\theta=\pi$. $d\phi$ doesn't make sense there and the coefficient multiplying it doesn't vanish. Integrate avoiding that half-line and your integral will work.
In flat space, I'll integrate $\theta$ from $0$ to $\theta_m$, $r_a$ to $r_b$, and over all $\phi$ and a time interval $T$. The normal to the boundary surface is in the $r$ direction everywhere except the surface at $\theta=\theta_m$ which is of course normal to the $\theta$ direction. This is exactly in the direction of $F^{\mu\nu}A_\mu$, and it will be the only contribution to the surface integral.
Your vector $v^\mu=F^{\mu\nu}A_\nu$ is in the $\theta$ direction and has magnitude $$|v|=\frac{Q^2}{r^3\sin\theta}(1-\cos\theta)$$ Now the surface integral is $$T \int_{r_a}^{r_b}|v|2\pi r\sin\theta_m dr= 2\pi T Q^2 (1-\cos\theta_m)\int_{r_a}^{r_b}\frac {dr}{r^2}\rightarrow 4\pi T Q^2 \int_{r_a}^{r_b}\frac {dr}{r^2}$$
Where in the last line I'm free to take $\theta_m=\pi$ now that the $A$ field is gone. This is exactly the volume integral of $Q^2/r^4$. We see the surface area $4\pi$ the time interval $T$ and a factor of $r^2$ coming from the measure of integration.
If you remove the machinery of curved spacetime, then you've just arrived at the usual subtlety involving magnetic monopoles. That is, if you assume $F = dA$ then you automatically have $dF = 0$, which includes Gauss's law for magnetism $\nabla \cdot \mathbf{B} = 0$ forbidding magnetic monopoles. Going to curved spacetime makes these equations look a little fancier but doesn't really change the logic.
To allow monopoles, we must either:
- switch to fiber bundle formalism, describing $A$ as a connection on a $U(1)$ bundle over your spacetime, or
- account for the "Dirac string", a singularity in $A$ that inevitably appears if you don't use fiber bundles, which in this case occurs at $\theta = \pi$.
Whatever option you choose, the technical fix is the same. If you use bundles, you'll have to cover your $r = \text{const}$ surface with two patches, and you'll pick up extra terms from the overlap. If you use the Dirac string, you have to cut out the part of the surface that the Dirac string passes through, and that yields a boundary which contributes on the right-hand side of the divergence theorem.