Maxwell's equation in curved spacetime - how come? And experimental evidence? I'm trying to understand the generalization of Maxwell's equations to curved spacetime. 
In FLAT (Minkowski) SPACETIME:
If we define the "four-potential" as $$\ (\mathcal{A}^{0},\mathcal{A}^{1},\mathcal{A}^{2},\mathcal{A}^{3})=(\frac{1}{c} V,A_{1},A_{2},A_{3})\ $$ and the "four-current" as $$(\mathcal{J}^{0},\mathcal{J}^{1},\mathcal{J}^{2},\mathcal{J}^{3})=(c\rho,J_{1},J_{2},J_{3})\ .$$ We define the electromagnetic field tensor as:
$$\mathcal{F}^{ab}=\frac{\partial \mathcal{A}^{b}}{\partial x^{a}}-\frac{\partial \mathcal{A}^{a}}{\partial x^{b}}$$
Then we can write the original four Maxwell equations compactly as these two tensorial equations:
$$\partial_{a} \mathcal{F}^{ab}=\mu_{0} \mathcal{J^{b}}$$
$$\partial_{c} \mathcal{F}^{ab} + \partial_{a} \mathcal{F}^{bc} + \partial_{b} \mathcal{F}^{ac}=0$$
Where $\partial_{a}$ is the regular partial derivative.
IN CURVED SPACETIME:
Let's say we're in some curved spacetime, with some metric tensor $g_{ab}$ that is very different from the Minkowski one. 
Is it true that: We pretty much reproduce the above except using $covariant\ derivatives$ instead of the regular partial derivative?
Let's suppose we've got some vector with components $u^{a}$, then I write the covariant derivative of the components as $\nabla_{b}u^{a}=\partial_{b}u^{a}+\Gamma^{a}_{\ bc}u^{c}$, and you know how we can generalize this to tensors of higher rank. Obviously, $\Gamma^{a}_{\ bc}$ are the Christoffel symbols of the second kind, which encode the curvature of the spacetime.
So then our Maxwell equations are actually:
$$\nabla_{a} \mathcal{F}^{ab}=\mu_{0} \mathcal{J^{b}}$$
$$\nabla_{c} \mathcal{F}^{ab} + \nabla_{a} \mathcal{F}^{bc} + \nabla_{b} \mathcal{F}^{ac}=0$$
......in a curved spacetime?
I want to understand where this comes from - that we replace the derivatives with covariant derivatives - is this a $hypothesis$? I am assuming yes.
And if so, have there been some sort of large-scale (astronomical) tests done to verify these equations in curved spacetime? Maybe playing with electromagnetic fields somehow near some gravitating body - I have no idea how this would be done, I am just curious.
Since locally any curvature is negligible, we get back the Maxwell equations on ordinary length scales, so I have a feeling it would be enormously difficult to be able test out these equations on scales large enough to notice the effects of curvature.
 A: As ACuriousMind pointed out, writing the equations for electromagnetism using differential forms,
$$dF = 0, \quad d \, \star \, F = J $$
results in coordinate-independent expressions, which are valid on any manifold. Now, in flat space we can write,
$$\partial_\mu A^{\mu\alpha} - \partial^\mu A^\alpha_\mu = 4\pi J^\alpha.$$
Naively, to transition to curved space one would change $\partial_\mu \to \nabla_\mu$, and you are right to suspect this is not necessarily valid. Since covariant derivatives do not commute, that is,
$$(\nabla_\mu \nabla_\nu - \nabla_\nu \nabla_\mu) A_{\lambda} = A_\sigma R^\sigma_{\lambda \mu \nu}$$
there are essentially two versions we could produce by swapping to covariant derivatives which are,
$$\nabla_\mu A^{\mu\alpha} - \nabla^\mu A^\alpha_\mu = 4\pi J^\alpha, \quad R^\alpha_\mu A^\mu + \nabla_\mu A^{\mu\alpha} - \nabla^\mu A^\alpha_\mu = 4\pi J^\alpha.$$
There is no general recipe to guarantee which is right, but Misner, Thorne and Wheeler argue you should only expect curvature terms for expressions derived from or involving equations of double covariant derivatives. In addition, there should be a physical reason for coupling to curvature if applying the rule to physically measurable quantities, like the field-strength rather than $A_\mu$. 

It can be shown charge conservation requires $\nabla_\gamma F_{\alpha\beta} + \nabla_\alpha F_{\beta\gamma} + \nabla_\beta F_{\gamma \alpha} = 0$ without additional curvature terms; this equation is essentially $dF =0$. That being said,
$$R^\alpha_\mu A^\mu + \nabla_\mu A^{\mu\alpha} - \nabla^\mu A^\alpha_\mu = 4\pi J^\alpha$$
is the correct form, and it can be shown this must be true if $\nabla_\beta F^{\alpha\beta} = 4\pi J^\alpha$ using $F^{\alpha\beta} = 2A^{[\alpha\beta]}$.

As another example of this reasoning, consider the conservation of the angular momentum vector, $\nabla_u S = 0$ for the Earth in flat space, along a world line of its c.o.m. Transitioning to curved space, we know other bodies' curvatures give rise to tidal forces on the Earth, and since the Earth is a sort of ellipsoid, there is a torque, so one would expect $\nabla_u S$ possessing some Riemann terms. For completeness, it turns out
$$\nabla_u S^\alpha = \epsilon^{\alpha\beta\gamma\delta} \bar{I}_{\beta\mu}R^\mu_{\nu \gamma \xi} u_\delta u^\nu u^\xi$$
where $u$ is a four-velocity and $\bar{I}$ is the trace-free part of the second moment of mass distribution.
