Is there differential form notation for Maxwell's equation in curved spacetime? In special relativity, Maxwell's equations may be written as
\begin{align*}
dF = 0, \\ \star\, d\star F = J.
\end{align*}
In four-vector notation, this translates to $\partial_{\mu}F^{\mu\nu} = J^{\nu}$ and $\partial_{[\lambda}F_{\mu\nu]} = 0$ where the brackets mean that there is a sum with indices cyclically permuted. 
Generalizing this to curved spacetime requires us to use covariant derivatives instead of partial derivatives, so we instead write $\nabla_{\mu}F^{\mu\nu} = J^{\nu}$ and $\nabla_{[\lambda}F_{\mu\nu]} = 0$.
Is there a way to write this using differential forms? Is there any reference that talks about this?
 A: 
Is there a way to write this using differential forms (such that it reduces to equations involving covariant derivatives rather than partial derivatives)?

It turns out that the differential form equations written in my original post are already what I wanted. This was unclear to me, so I will show that the differential form Maxwell's equations reduce to the covariant Maxwell's equations.
Theorem. If $F$ is a differential $2$-form and $J$ is a differential $1$-form such that $dF = 0$ and $\star\, d\star F = J$, then in any coordinate system we have $\nabla_{[\lambda}F_{\mu\nu]} = 0$ and $\nabla_{\mu}F^{\mu\nu} = J^{\nu}$.

Proof 1. Given any point $p\in M$, choose Riemann normal coordinates $(x^{i})$ centered at $p$. These yield $g_{\mu\nu}(p) = \eta_{\mu\nu}$ and $\Gamma^{\mu}_{\nu\lambda}(p) = 0$. In these coordinates, the differential form version of Maxwell's equations reduce to
$$ \partial_{\mu}F^{\mu\nu}(p) = J^{\nu}(p) \quad\text{and}\quad \partial_{[\lambda}F_{\mu\nu]}(p) = 0. $$
We also have $\partial_{\mu} = \nabla_{\mu}$ at point $p$ in these coordinates, so this means the above is the same as
\begin{align}
\nabla_{\mu}F^{\mu\nu}(p) = J^{\nu}(p) \quad\text{and}\quad\nabla_{[\lambda}F_{\mu\nu]}(p) = 0 
\end{align}
in these coordinates.
But by general covariance, a transformation to any other coordinates must preserve the form of the above two equations, so they hold for all coordinates where $p$ is an arbitrary point.
$$\tag*{$\blacksquare$}$$

Proof 2. We will do explicit computations. Let $F$ be a $2$-form and $J$ be a $1$-form such that
$$ dF = 0 \qquad\text{ and }\qquad \star\, d\star F = J. $$
Take any coordinate system and write $F = \frac{1}{2}F_{\mu\nu}\,dx^{\mu}\wedge dx^{\nu}$ where $F_{\mu\nu}$ shall be antisymmetric in indices $\mu$ and $\nu$.
Then
\begin{align}
dF = \frac{1}{2}\partial_{\lambda}F_{\mu\nu}\, dx^{\lambda}\wedge dx^{\mu}\wedge dx^{\nu} = 0.
\end{align}
This implies that for a given set of distinct, fixed $\lambda, \mu, \nu$,
$$ \frac{1}{2}(\partial_{\lambda}F_{\mu\nu} - \partial_{\lambda}F_{\nu\mu} + \partial_{\mu}F_{\nu\lambda} - \partial_{\mu}F_{\nu\lambda} + \partial_{\nu}F_{\lambda\mu} - \partial_{\nu}F_{\mu\lambda}) = 0. $$
By antisymmetry of $F_{\mu\nu}$, the above reduces to
$$ \partial_{\lambda}F_{\mu\nu} + \partial_{\mu}F_{\nu\lambda} + \partial_{\nu}F_{\lambda\mu} = 0. $$
Given the last equation, we add & subtract "Christoffel terms" to the left-hand side and obtain the following:
\begin{align}
& \partial_{\lambda}F_{\mu\nu} - \Gamma^{\alpha}_{\lambda\mu}F_{\alpha\nu} - \Gamma^{\alpha}_{\lambda\nu}F_{\mu\alpha} + \Gamma^{\alpha}_{\lambda\mu}F_{\alpha\nu} + \Gamma^{\alpha}_{\lambda\nu}F_{\mu\alpha} \\
& \partial_{\mu}F_{\nu\lambda} - \Gamma^{\alpha}_{\mu\nu}F_{\alpha\lambda} - \Gamma^{\alpha}_{\mu\lambda}F_{\nu\alpha} + \Gamma^{\alpha}_{\mu\nu}F_{\alpha\lambda} + \Gamma^{\alpha}_{\mu\lambda}F_{\nu\alpha} \\
& \partial_{\nu}F_{\lambda\mu} - \Gamma^{\alpha}_{\nu\lambda}F_{\alpha\mu} - \Gamma^{\alpha}_{\nu\mu}F_{\lambda\alpha} + \Gamma^{\alpha}_{\nu\lambda}F_{\alpha\mu} + \Gamma^{\alpha}_{\nu\mu}F_{\lambda\alpha} = 0.
\end{align}
Using the symmetry of the Christoffel symbols $\Gamma^{\alpha}_{\mu\nu} = \Gamma^{\alpha}_{\nu\mu}$ and the antisymmetry of the Faraday tensor $F_{\mu\nu} = -F_{\nu\mu}$, the right two columns above cancel out. The remaining terms above reduce to
$$ \nabla_{\lambda}F_{\mu\nu} + \nabla_{\mu}F_{\nu\lambda} + \nabla_{\nu}F_{\lambda\mu} = 0. $$
This gives us one of the Maxwell's equations.
To get the other equation, we'll invoke two lemmas listed at the end of this post.
To start things off, we compute $\star\, d\star F$ using the explicit formula for the Hodge star operator:
\begin{align}
\star\, F &= \frac{\sqrt{|g|}}{2!} g^{\mu_{1}\nu_{1}}g^{\mu_{2}\nu_{2}} \, \epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}} \left(\frac{1}{2} F_{\mu_{1}\mu_{2}} \right) \, dx^{\nu_{3}}\wedge dx^{\nu_{4}}, \\
d\star F &= \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\nu_{1}\nu_{2}}}{4} \right) \epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}\, dx^{\lambda}\wedge dx^{\nu_{3}}\wedge dx^{\nu_{4}}, \\
\star\, d\star F &= \frac{\sqrt{|g|}}{1!} g^{\lambda\kappa_{1}}g^{\nu_{3}\kappa_{2}}g^{\nu_{4}\kappa_{3}}\, \epsilon_{\kappa_{1}\kappa_{2}\kappa_{3}\kappa_{4}} \cdot \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\nu_{1}\nu_{2}}}{4} \right) \epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}} \, dx^{\kappa_{4}}.
\end{align}
Given $\star\, d\star F = J$, we find
$$ \frac{\sqrt{|g|}}{1!} g^{\lambda\kappa_{1}}g^{\nu_{3}\kappa_{2}}g^{\nu_{4}\kappa_{3}}\, \epsilon_{\kappa_{1}\kappa_{2}\kappa_{3}\kappa_{4}} \cdot \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\nu_{1}\nu_{2}}}{4} \right) \epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}} \, dx^{\kappa_{4}} = J_{\kappa_{4}}dx^{\kappa_{4}}. $$
Now we equate the components and apply $g^{\sigma\kappa_{4}}$ to both sides:
\begin{align*}
\sqrt{|g|}\, g^{\lambda\kappa_{1}}g^{\nu_{3}\kappa_{2}}g^{\nu_{4}\kappa_{3}}g^{\sigma\kappa_{4}} \cdot \epsilon_{\kappa_{1}\kappa_{2}\kappa_{3}\kappa_{4}}\epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}} \cdot \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\nu_{1}\nu_{2}} }{4} \right) = J^{\sigma}.
\end{align*}
By Lemma 1, we have
\begin{align*}
\sqrt{|\det g_{\mu\nu}|}\, (\det g^{\mu\nu})\, {\epsilon}^{\lambda\sigma\nu_{3}\nu_{4}} \epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}} \cdot \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\nu_{1}\nu_{2}} }{4} \right) = J^{\sigma}.
\end{align*}
Note that $(g^{\mu\nu}) = (g_{\mu\nu})^{-1}$ so the prefactor is just $1/\sqrt{|g|} = 1/\sqrt{|\det g_{\mu\nu}|}$.
To handle the Levi-Civita symbols, one can reason that ${\epsilon}^{\lambda\sigma\nu_{3}\nu_{4}} \epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}} = 2\delta^{\lambda}_{\nu_{1}}\delta^{\sigma}_{\nu_{2}} - 2\delta^{\lambda}_{\nu_{2}}\delta^{\sigma}_{\nu_{1}}$, and so
\begin{align*}
\frac{1}{\sqrt{|g|}}\cdot (\delta^{\lambda}_{\nu_{1}}\delta^{\sigma}_{\nu_{2}} - \delta^{\lambda}_{\nu_{2}}\delta^{\sigma}_{\nu_{1}}) \cdot \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\nu_{1}\nu_{2}} }{2} \right) = J^{\sigma} \\
\frac{1}{\sqrt{|g|}}\, \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\lambda\sigma}}{2} \right) -  \frac{1}{\sqrt{|g|}}\, \partial_{\lambda}\left( \frac{\sqrt{|g|}\, F^{\sigma\lambda}}{2} \right) = J^{\sigma} \\
\frac{1}{\sqrt{|g|}}\, \partial_{\lambda}\left( \sqrt{|g|}\, F^{\lambda\sigma} \right) = J^{\sigma}.
\end{align*}
By Lemma 2, this gives the desired formula.
$$\tag*{$\blacksquare$}$$

Lemma 1. Given any $\lambda_{1}, \lambda_{2}, \lambda_{3}, \lambda_{4}\in \{0, 1, 2, 3\}$ and any pseudo-Riemannian metric $g$, we have
$$ g^{\lambda_{1}\kappa_{1}}g^{\lambda_{2}\kappa_{2}}g^{\lambda_{3}\kappa_{3}}g^{\lambda_{4}\kappa_{4}}\epsilon_{\kappa_{1}\kappa_{2}\kappa_{3}\kappa_{4}} =  (\det g^{\mu\nu}){\epsilon}^{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}} $$
where $\epsilon_{\kappa_{1}\kappa_{2}\kappa_{3}\kappa_{4}}$ and ${\epsilon}^{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}$ are the signs of permutations $(0, 1, 2, 3)\rightarrow (\kappa_{1}, \kappa_{2}, \kappa_{3}, \kappa_{4})$ and $(0, 1, 2, 3)\rightarrow (\lambda_{1}, \lambda_{2}, \lambda_{3}, \lambda_{4})$, respectively (and either one is zero if it has repeated indices).
Lemma 2. For any rank-$2$ tensor $F^{\mu\nu}$, we have
$$ \frac{1}{\sqrt{|g|}}\partial_{\lambda} \left(\sqrt{|g|}\, F^{\lambda \sigma} \right) = \nabla_{\lambda} F^{\lambda\sigma}. $$
A: Differential forms are natural objects that depend solely on the smooth structure, therefore they are valid in any spacetime.
Note that the Hodge star $\star$ does, however, depend on a metric tensor. However, any metric tensor works with it. Therefore, Maxwell's equations in differential form notation are valid in any spacetime without modifications, assuming the Hodge star in the equations refer always to the metric tensor that defines the geometry of your spacetime.
Note that even in standard tensor calculus notation, Maxwell's equations can be cast in the form $$ \partial_{[\kappa} F_{\mu\nu]}=0 \\ \partial_\nu\mathfrak F^{\mu\nu}=\mu_0\mathfrak{j}^\mu, $$ where $\mathfrak F^{\mu\nu}=F^{\mu\nu}\sqrt{-g}$ and $\mathfrak j^\mu=j^\mu\sqrt{-g}$, which do not involve covariant derivatives.
