I have difficulties to understand how to solve the Maxwell equations on curved spacetime. I want to solve the equations in the weak regime $g_{\mu\nu}=\eta_{\mu\nu}+h{\mu\nu},~ h_{\mu\nu}\ll 1$ without sources $J_{\mu}=0$. I use $\eta_{\mu\nu}=\text{diag}(1,-1,-1,-1)$.

The Maxwell equations in curved spacetime are:

$F^{\mu\nu}_{~~~~~||\mu}=0 ~(*)$

$\epsilon^{\mu\nu\lambda\kappa}F_{\lambda\kappa||\nu}=0 ~(**)$

where $||$ means covariant derivative. In the following I focus on $(*)$. I can rewrite $(*)$ as:

$0=F^{\mu\nu}_{~~~~~||\mu}$=$F^{\mu\nu}_{~~~~~|\mu}+\Gamma^{\mu}_{\mu\rho}F^{\rho\nu}=\frac{1}{\sqrt{g}}\partial_{\mu}(\sqrt{g}F^{\mu\nu})$,

where I have used the definition of the Christoffel symbols $\Gamma^{\mu}_{\mu\rho}=\frac{1}{\sqrt{g}}\partial_{\rho}\sqrt{g}$ and $g=-\text{det}(g^{\mu\nu})$. Summarized I get

$\partial_{\mu}(\sqrt{g}F^{\mu\nu})=0,~~~ (***)$

I furthermore use the transverse traceless gauge

$h_{\mu\nu}=\begin{pmatrix}0&0&0&0\\0&h_+&h_\times&0\\0&h_\times&-h_+&0\\0&0&0&0\end{pmatrix}$

and assume that $h_{\mu\nu}$ is a background field. The coupling to the Einstein Equations is neglected. I can expand

$g=1-h_+^2-h_\times^2+\mathcal{O}(h^4)$

I now can divide $F^{\mu\nu}$ in powers of $h$: $F^{\mu\nu}=F_{(0)}^{\mu\nu}+F_{(1)}^{\mu\nu}+\cdots$ and expand all fields in equation $(***)$

The zeroth order yields normal electrodynamics:

$\partial_{\mu}F^{\mu\nu}_{(0)}=0$

However the first order field strength tensor also yields normal electrodynamics:

$\partial_{\mu}F^{\mu\nu}_{(1)}=0$

$\textbf{Now comes my question}$

If I rewrite $(***)$ to

$\partial_{\mu}(\sqrt{g}g^{\mu\alpha}g^{\nu\beta}F_{\alpha\beta})=0,~~~ (****)$

and do the expansions I get non trivial first order equations for $F^{\mu\nu}_{(1)}$:

$\partial_{\beta}F^{\mu\beta}_{(1)}=-(\partial_{\nu}h^{\mu}_{~\alpha})F^{\alpha\nu}_{(0)}-\partial_\nu(h^{\nu}_{~\beta}F^{\mu\beta}_{(0)})-h^{\mu}_{~\alpha}\partial_{\beta}F^{\alpha\beta}_{(0)}$

How can this be? Does it have to do something with the choice of the system where we do the calculation? I am interested in the case that a gravitational wave hits the Earth. What would be the physical electromagnetic fields? Would they be components of $F^{\mu\nu}_{(1)}$ or components of $F_{(1)\mu\nu}$? I think this should not matter since we have to raise and lower the indices with $\eta$ in order not to change the order of the tensor!? My question also arises due to this paper. They solve my equation $(****)$ what is their equation (9) in the paper. However to get to the physical electromagnetic fields they say that one has to do another transformation, c.f., their equation (A.9)-(A.11). I simply do not understand this.

I have difficulties to understand how to solve the Maxwell equations on curved spacetime. I want to solve the equations in the weak regime $g_{\mu\nu}=\eta_{\mu\nu}+h{\mu\nu},~ h_{\mu\nu}\ll 1$ without sources $J_{\mu}=0$. I use $\eta_{\mu\nu}=\text{diag}(1,-1,-1,-1)$.

The Maxwell equations in curved spacetime are $$F^{\mu\nu}_{~~~~~||\mu}=0 \tag{$*$}$$ and $$\epsilon^{\mu\nu\lambda\kappa}F_{\lambda\kappa||\nu}=0 \tag{$**$}$$ where $||$ means covariant derivative. In the following I focus on $(*)$. I can rewrite $(*)$ as $$0=F^{\mu\nu}{}_{||\mu}=F^{\mu\nu}{}_{|\mu}+\Gamma^{\mu}{}_{\mu\rho}F^{\rho\nu}=\frac{1}{\sqrt{g}}\partial_{\mu}(\sqrt{g}F^{\mu\nu}),$$ where I have used the definition of the Christoffel symbols $\Gamma^{\mu}{}_{\mu\rho}=\frac{1}{\sqrt{g}}\partial_{\rho}\sqrt{g}$ and $g=-\text{det}(g^{\mu\nu})$. In summary, I get $$\partial_{\mu}(\sqrt{g}F^{\mu\nu})=0, \tag{$***$}$$ and I furthermore use the transverse traceless gauge $$h_{\mu\nu}=\begin{pmatrix}0&0&0&0\\0&h_+&h_\times&0\\0&h_\times&-h_+&0\\0&0&0&0\end{pmatrix}$$ and assume that $h_{\mu\nu}$ is a background field. The coupling to the Einstein Equations is neglected. I can expand $$g=1-h_+^2-h_\times^2+\mathcal{O}(h^4)$$

I now can split $F^{\mu\nu}$ in powers of $h$: $F^{\mu\nu}=F_{(0)}^{\mu\nu}+F_{(1)}^{\mu\nu}+\cdots$ and expand all fields in equation $(***)$

The zeroth order yields normal electrodynamics, $$\partial_{\mu}F^{\mu\nu}_{(0)}=0.$$ However the first order field strength tensor also yields normal electrodynamics, $$\partial_{\mu}F^{\mu\nu}_{(1)}=0$$

Now comes my question:

If I rewrite $(***)$ to $$\partial_{\mu}(\sqrt{g}g^{\mu\alpha}g^{\nu\beta}F_{\alpha\beta})=0, \tag{$****$}$$ and do the expansions I get non-trivial first-order equations for $F^{\mu\nu}_{(1)}$, $$\partial_{\beta}F^{\mu\beta}_{(1)}=-(\partial_{\nu}h^{\mu}{}_{\alpha})F^{\alpha\nu}_{(0)}-\partial_\nu(h^{\nu}_{~\beta}F^{\mu\beta}_{(0)})-h^{\mu}{}_{\alpha}\partial_{\beta}F^{\alpha\beta}_{(0)}$$

How can this be? Does it have to do something with the choice of the system where we do the calculation? I am interested in the case that a gravitational wave hits the Earth. What would be the physical electromagnetic fields? Would they be components of $F^{\mu\nu}_{(1)}$ or components of $F_{(1)\mu\nu}$? I think this should not matter since we have to raise and lower the indices with $\eta$ in order not to change the order of the tensor!? My question also arises due to this paper. They solve my equation $(****)$ what is their equation (9) in the paper. However to get to the physical electromagnetic fields they say that one has to do another transformation, c.f., their equation (A.9)-(A.11). I simply do not understand this.