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.