I know that we can write Maxwell's equations in the covariant form, and this covariant form can be considered as a generalization of these equations in curved spacetime if we replace ordinary derivatives with covariant derivatives. But I have read somewhere that this generalization is not unique and it is just the simplest one. Can anyone introduce some resources about this subject and how electromagnetism and Maxwell's equations are generalized to curved spacetime?
Electromagnetism – Maxwell’s Equations Explained in Curved Spacetime
curvatureelectromagnetismgeneral-relativitymaxwell-equationsmetric-tensor
Related Solutions
In the language of differential forms, the Maxwell-Lorentz equations are simply $$\begin{eqnarray*}\mathrm{d}\!\star\!F = J/\lambda_0 &\text{,}\quad&\mathrm{d}F = 0\text{,}\end{eqnarray*}$$ where $1/\lambda_0$ is the characteristic impedance of free space, and can be fixed to $1$ in appropriate units.
From that point of view, all you need to "move" to curved spacetime is realize that the only place where the metric comes up is the in the Hodge dual operator $\star$.
This was what I expected initially, but I got scared by the existence of the exterior covariant derivative -- I can't figure out why it is needed if the exterior derivative already behaves geometrically.
As already said, the differential structure does not care about either the metric or the manifold-connection, and hence not the curvature. That's simply a matter of how they are defined, and therefore curved spacetime wouldn't care about the exterior derivative, but only have a different Hodge dual.
However, there is an informal way to motivate this by analogy, which is what I gather you're asking for (this is not clear, so I'm guessing). Recall that the notion of covariant derivative for scalar fields is particularly trivial and independent of curvature: $$\nabla_u \phi = u^\mu\phi_{,\mu}\text{.}$$ The same kind of thing happens with the covariant exterior derivative. Recall the $k$-forms are certain types of maps in the form $$\omega_p: (T_pM)^k\to\mathbb{R}\text{,}$$ and are therefore scalar-valued.
On the other hand, one can consider a kind of generalized form that takes values in some vector bundle over over the manifold, and try to have an "exterior derivative" with $\mathrm{d}$ acting on each component in an arbitrary basis. This generally does not work without a connection on that vector bundle, with a notable exception where the values taken are are actually scalars, since then all possible bases for this $1$-dimensional vector space are just proportional to one another.
So the short of it is that the exterior derivative makes sense regardless of curvature and the covariant exterior derivative doesn't do anything interesting for ordinary scalar-valued forms.
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.
Best Answer
Here is why I doubt there are other ways to generalize Maxwell's equations to curved spacetime.
Special relativity was obtained from the invariance of the speed of light. In special relativity, the electric field is not a vector field, and the magnetic field is not a pseudovector, but that they transform as the components of a two-form $F_{ab} = \partial_a A_b -\partial_b A_a$, where the four-vector $A_a$ contains the scalar and vector potentials.
Maxwell's equations become $$d F=0$$ $$d\ast F=J$$
When moving to curved spacetimes, they remain the same, since the Hodge dual $\ast$ is defined at each point $p$ of the manifold, on $\wedge^\bullet T^\ast_p$. When expressed in this form, the covariant derivative is not involved, although the metric is involved in the $\ast$ operator.
While I think the generalization of Maxwell's equations to curved spacetime is very rigid and I see no choice based on simplicity here, it is known that there are modified (nonlinear) versions, like the Born-Infeld theory. But they did not originate because of some freedom of generalizing Maxwell's equations to curved spacetimes.