In the framework of geometric optics, which sufficiently describes the effects in the near-earth environment, Faraday rotation is purely a reference frame effect. A simple encoding procedure could mitigate the Faraday phase error. However, geometric optics approximation is not sufficient to describe the propagation of waves of large but finite frequencies. So, we outline the technique to solve the equations for the propagation of an electromagnetic wave up to the subleading order geometric optics expansion in curved spacetimes. This could be achieved in two nontrivial steps. First, we need to construct a set of parallel propagated null tetrads in curved spacetimes. A general procedure exists for solving the parallel transport equations in Petrov-D spacetimes, which contain an extra constant of motion, also called the Carter constant. Two of the components of such tetrad give the propagation and polarization of an electromagnetic wave in geometric optics approximation. Then we should use the parallel propagated tetrad to solve the modified trajectory equation in curved spacetimes. The wavelength-dependent deviation of the electromagnetic waves is observed, which gives the mathematical description of the gravitational Faraday effect in curved spacetimes.