Is there a way to 'derive' the equation of massless particle geodesics in curved space-time starting from the field equations of electromagnetism in curved space-time?
I'm imagining something along the following lines:
Choose/Define what it means to track a propagating disturbance of electromagnetic fields (some suitable pulse, "bump", wave-packet, etc.)
Start with the field equations for electromagnetism in a fixed potentially curved background space-time. (Evolving in a vacuum -- there can be non-zero electromagnetic fields, but no charges, nor polarizable material to travel through, etc.)
Make approximations to the field evolution as needed.
Apply definition of tracking a propagating "bump" of electromagnetic fields.
Arrive at equation of motion for a massless particle in curved space-time: $$ g_{\mu\nu} \frac{\partial x^\mu}{\partial \lambda} \frac{\partial x^\nu}{\partial \lambda} = 0$$ where $\lambda$ is a world-line parameter.
Or maybe it is easier to choose the definitions such that we can just show an approximation works at the level of the Lagrangian without need to solve for any equations of motion? But I'd think any definition of the path of some object would require at least discussing the equations of motion.
I'm really hoping that the right sequence of nice definitions and approximations make this not too difficult, but any attempt I make at the very first step seems inappropriate and suitable only for solving numerically.
Please Note:
Nothing here involves GR. I'd like to keep this 'derivation' of EM in curved space-time $\to$ massless particle geodesics as independent of any theory stating what the background is. Please consider the background as just static and taken as an input, for instance by specifying the metric $g_{\mu\nu}$. So for example, unless it is necessary for the equations to work out, please don't even assume that $R=0$ for the background as in GR vacuum. (Actually if it doesn't work out unless $R=0$, that would be very interesting to know in its own right.)
Some links:
Variational principle for a point particle (massive or massless) in curved space
Lagrangian for relativistic massless point particle
https://en.wikipedia.org/wiki/Maxwell%27s_equations_in_curved_spacetime
My attempt with user Someone's suggestion:
If I'm understanding correctly, the idea is to zoom in and look at a wave on a small patch of space-time. We can choose something like a plane wave with the frequency high enough to have the wavelength be small compared to the length scale of the curvature, and with this approximation we can see how the wavefront propagates and curves. This captures the idea of tracking a small wave packet or bump in the fields.
So consider $A_\mu = \epsilon_\mu e^{i\phi}$ on some patch of space-time. The parameters $\phi$ and $\epsilon_\mu$ are functions of position, and we define $k_\mu \equiv \nabla_\mu \phi$. Then plug this into the electrodynamic field equations to get constraints on $k_\mu$.
To start easy, I'll try this in a space-time with no Ricci curvature ($R_{ab} = 0$) and use the Lorenz gauge, so the equations of electrodynamics with no charges are $\nabla^a \nabla_a A_b = 0$.
The Lorenz gauge, $\nabla^a A_a = 0$, provides these constraints: $$\begin{equation} e^{i\phi} \nabla^a \epsilon_a + i \epsilon_a e^{i\phi} \nabla^a \phi = 0 \ \ \longrightarrow \ \ \nabla^a \epsilon_a = 0,\ \epsilon_a \nabla^a \phi = 0 \tag{1} \end{equation} $$
Putting in our form for $A_\mu$ into the electrodynamic equations results in: $$\begin{array} \ 0 = \nabla^a \nabla_a \left( \epsilon_b e^{i\phi} \right) &=& \nabla^a \left( e^{i\phi}\nabla_a \epsilon_b + i \epsilon_b e^{i\phi} \nabla_a \phi \right) \\ &=& \left( i e^{i\phi} \nabla^a \phi \right) \left( \nabla_a \epsilon_b \right) + e^{i\phi}\nabla^a \nabla_a \epsilon_b \\ & & + i e^{i\phi} \left( \nabla^a \epsilon_b \right) \left( \nabla_a \phi \right) - \epsilon_b e^{i\phi} \left( \nabla^a \phi \right) \left( \nabla_a \phi \right) +\left(i \epsilon_b e^{i\phi} \right) \left( \nabla^a \nabla_a \phi \right) \tag{2} \end{array}$$ Multiplying by $e^{-i\phi}$ and then collecting the real and imaginary terms: $$\begin{gather} \nabla^a \nabla_a \epsilon_b - \epsilon_b \left( \nabla^a \phi \right) \left( \nabla_a \phi \right)= 0 \tag{3a}\\ \left( \nabla^a \phi \right) \left( \nabla_a \epsilon_b \right) + \left( \nabla^a \epsilon_b \right) \left( \nabla_a \phi \right) + \epsilon_b \left( \nabla^a \nabla_a \phi \right) = 0 \tag{3b} \end{gather}$$
Now we use the "eïkonal approximation", noting that we can increase the frequency of our wave arbitrarily high so that the derivative terms on $\phi$ dominate. As long as the phase of the wave is oscillating faster than its amplitude is changing, the terms with a derivative of the amplitude can be considered negligible.
This approximation simplifies the results to: $$\begin{gather} \epsilon_b \left( \nabla^a \phi \right) \left( \nabla_a \phi \right)= 0 \ \ \longrightarrow \ \ k^a k_a = 0 \tag{4a}\\ \epsilon_b \left( \nabla^a \nabla_a \phi \right) = 0 \ \ \longrightarrow \ \ \nabla^a k_a = 0 \tag{4b} \end{gather}$$
Now noting that $$\nabla_a \nabla_b \phi = \nabla_b \nabla_a \phi \ \ \longrightarrow \ \ \nabla_a k_b = \nabla_b k_a \tag{5}$$ we can see that $k^a k_a = 0$ is equivalent to the geodesic equation (Form of the geodesic equation): $$0 = \frac{1}{2} \nabla_a \left( g^{bc} k_b k_c \right) = g^{bc} k_b \nabla_a k_c = g^{bc} k_b \nabla_c k_a = k^c \nabla_c k_a = 0 \tag{6}$$
If Ricci curvature was non-zero, the field equations would be $$\nabla^a \nabla_a A_b = {R^a}_b A_a \tag{7}$$ In the original constraints this would just change the zero term in eq 3a, but is still negligible with the eikonal approximation. So the result is the same geodesic equation.
