In special relativity one may show that a plane wave solution of Maxwell's equations (in a vacuum), of the form $A^a=C^a\mathrm{e}^{\mathrm{i}\psi}$ has the following properties: The normal $k:=\mathrm{d}\psi$ to the surfaces of constant $\psi$ is a null vector and the integral curves of $k$ are null geodesics. Here $A$ is the electromagnetic vector potential, $C$ is a constant vector and $\psi$ is some function.
This analysis is possible because of the relatively simple form of the Maxwell equations in flat space, $\partial^a\partial_a A^b=0$ (Lorenz gauge assumed). However, in curved spacetime, we have an extra term involving the Ricci tensor that is irrelevant for SR: $$\nabla^a\nabla_aA^b=R^b{}_aA^a,$$ where $\nabla$ is the Levi-Civita connection of our spacetime $(\mathcal{M},g)$ and $R_{ab}$ its Ricci curvature.
The textbook treatment is to now look at solutions of the form $A^a=C^a\mathrm{e}^{\mathrm{i}\psi}$ where the covariant derivatives of $C$ are "small." In order to obtain the condition for $\mathrm{d}\psi$ to be null and autoparallel ($\nabla_a\psi\nabla^a\psi=0$), one must ignore the term $\nabla_b\nabla^b C^a$ as well as the Ricci tensor term. The missing details may be found in [1], sections 4.2 and 4.3. This approximation is called the geometric optics approximation.
Ref. [2] gives the following characteristic lengths for ray optics (section 2.8):
The wavelength $\lambda$.
The typical length $L$ over which the amplitude, polarization and wavelength of the wave vary significantly.
A typical "radius of curvature," which can be taken to be $$R:=\lvert\text{typical component of the Riemann tensor in a typical local inertial system}\rvert^{-1/2}$$
The region of validity for geometric optics is then $\lambda\ll L$ and $\lambda \ll R$.
Question:
Since one must ignore terms in the preceding analysis, do light rays not actually follow null geodesics in GR? Is $\mathrm{d}\psi$ even null? Furthermore, how do wave-like solutions of Maxwell's equations behave over length-scales greater than those given in [2] (i.e., over length scales where the curvature of $\mathcal{M}$ can vary greatly and rapidly). In particular, to what extent do they travel along null geodesics?
Please note: The argument that massless particles travel along null geodesics in flat space so the same must be true (by the Equivalence Principle) in curved space is not an answer to this question. I'm asking about wave-like solutions to the Maxwell equations. Light, classically, is just a wave solution of the vacuum Maxwell equations. Any answer should include (or reference) a rigorous analysis of Maxwell's equations. This is not a question that can be answered well by stating a few equations everyone knows and nitpicking at the words in the OP. (Apparently one member of the site did not have this impression, so I'm making it more clear.)
References:
[1] Wald, R.M. General Relativity. Chicago University Press, 1984.
[2] Straumann, N. General Relativity. Springer, 2013.
For clarity I think is best to start with Minkowski spacetime.
The equation we are trying to solve to understand the radiation of a point particle is: $$\square A^{b}=j^{b}$$
with the gauge $\nabla_{a}A^{a}=0$ and $j^{b}$ is the current density.
The potential \begin{eqnarray} A^{b}(t,x)&=&\int G^{b}_{a}(t,x,t'x')j^{a}(t',x')dx'^{3}dt\\ &=&\int\delta_{a}^{b}\delta(t−t′−|x − x′|)∕|x − x′|j^{a}(t',x')dx'^{3}dt' \end{eqnarray}
where $G^{b}_{a}$ is the Green function with support in the past light cone. In fact, the potential $A^{b}(t,x)$ only depends in the single event $(t',x')$ in the past which is the intersection between the null cone from $(t,x)$ and the world line of the particle.
Now in curved spacetime the generalization \begin{eqnarray} A^{b}(t,x)&=&\int G^{b}_{a}(t,x;t'x')j^{a}(t,x)dV\\ &=&\int\delta_{a}^{b}\delta(\gamma(t,x,t'x'))∕\Gamma(x − x ′)|j^{a}(t',x')\sqrt(g)dx'^{3}dt' \end{eqnarray}
where $\gamma$ is the null geodesic between the two points $(t,x),(t',x')$ and $\Gamma$ is the distance with respect the induced metric of a suitable spacelike surface that contains $x,x'$ does not work.
In general curved spacetimes the retarded Green function would depend on the whole causal past cone and not only the past light cone. This dependence comes from the interaction with the curvature and is related with the extra terms that you point out that vanish for Minkowski.
Therefore the potential is not only defined by the information that travel along the null geodesics but depends on the whole past of the particle. Nevertheless, singularities of the field travel along null geodesics globally. This is the content of the propagation of singularity theorems for linear hyperbolic systems and is related with the geometric optics limit.
As you required rigorous analysis I will point you to some papers with appropriate calculations:
Section 1.4 of http://relativity.livingreviews.org/open?pubNo=lrr-2011-7&page=articlese1.html
http://arxiv.org/abs/1108.1825
http://arxiv.org/abs/gr-qc/0008047
Also notice that my answer is just about electromagnetism in curved spacetime. To talk about General Relativity we would need to solve also for the Einstein's Equations. The point particle will affect the metric as self force corrections to the background metric. These type of corrections are treated in depth in the first reference.
