\section{Introduction} \label{sec:intro}

\paragraph{Beyond geometrical optics.}

Cosmological observations and hence our conclusions about the universe at large are mainly based on light (although neutrinos, charged particles and, most recently, gravitational waves are also observed). The analysis of light propagation is usually based on the geometrical optics approximation \cite{PhoneBook} (page 570), \cite{Lensing} (page 93). It assumes that the wavelength is much smaller than the other scales in the problem, namely the length over which the amplitude and polarisation change, and the spacetime curvature radius. In this case light propagation can be described in terms of a local plane wave orthogonal to a null geodesic ray. This is a version of the WKB approximation and is also called the eikonal approximation.

The geometrical optics approximation is known to fail in gravitational lensing when the lens size is comparable to the wavelength, as the wavefront is not constant over a few wavelengths. Instead, the wave-optics approximation is used (see e.g.\  \cite{WaveOptics_4, WaveOptics_5, WaveOptics_2,WaveOptics_3, Turyshev:2021eoh} and references therein\footnote{We focus on light; for discussion of geometrical optics and wave optics for gravitational waves, see \cite{GW_1,GW_2,GW_3, WaveOptics_GW, GW_4,GW_5}.}). A different issue is related to spacetime curvature. In general, vacuum solutions of the Maxwell equation in curved spacetime do not follow null geodesics, as is well known \cite{DeWitt:1960fc, MaxwellEFE_1,MaxwellEFE_4, Copi:2020qur, critique_1}. Expanding in the short wavelength approximation beyond the leading terms \cite{HigherOrder_1, Dolan:2018ydp}, the covariant derivatives in the equation of motion lead to spin-dependent light propagation, known as the gravitational Hall effect \cite{Oancea:2019pgm, Oancea:2020khc, Frolov:2020uhn}. However, the Maxwell equation in curved spacetime also involves the Ricci tensor, and the related curvature radius $R$ can be much smaller than the wavelength.

If dark matter consists of well-separated particles with mass $m$, approximated as spheres of constant density with radius equal to the Compton wavelength $\lambda_\text{c}=2\pi/m$, we have $R^{-2}\sim\rho/M_{{}_{\mathrm{Pl}}}^2\sim 10^{-3} m^4/M_{{}_{\mathrm{Pl}}}^2$, where $M_{{}_{\mathrm{Pl}}}$ is the Planck mass and $\rho$ is the local energy density. The lowest energy of observed cosmic photons is determined by the fact that the ionosphere is opaque to light frequencies below $\sim$ 10 MHz, corresponding to light wavelength $\sim10$ m and energy $E\sim10^{-8}$ eV. (None of the radio telescopes in space have been sensitive to such small frequencies.) Setting $\rho/M_{{}_{\mathrm{Pl}}}^2=E^2$ gives $m\sim10$ GeV. If dark matter is heavier than this, the curvature is larger than the smallest observed photon energy, so curvature effects may be relevant for some electromagnetic observations. By this estimate, the effect of electrons is negligible. The density of nuclei is $\rho\sim10^{-3}$ GeV$^{4}$, which gives $\rho/(E^2 M_{{}_{\mathrm{Pl}}}^2)\sim10^{-6}$, also negligible. Also, light waves cannot pass through charged particles in the same way as through dark matter. The highest energy photons observed have $E\sim10^5$ GeV, and the corresponding dark matter mass is $m\gtrsim10^{12}$ GeV. Cosmic microwave background (CMB) photons at redshift $z$ have energy $\sim10^{-3}(1+z)$ eV, so for $m\gtrsim10^4\sqrt{1+z}$ GeV, the curvature radius at dark matter particle position is smaller than the photon wavelength. The cosmological average number density of dark matter particles is $\sim (1+z)^3 (m/\text{GeV})^{-1}$m$^{-3}$, so for all these masses the distance between dark matter particles is many orders of magnitude larger than their size, even when local clustering is taken into account. If the curvature radius is smaller than the photon wavelength, the photon travels through an environment characterised by well-separated, narrow and high curvature spikes, rather than a smooth background as assumed in geometrical optics. We will see that the naive estimate above based on the curvature radius is too conservative, as the angular diameter distance involves not just the density but also its derivative.

While the effect of clumping on light propagation beyond perturbation theory has been considered before, including with discrete distributions of matter \cite{Zeldovich:1964, Dashevskii:1965, Dashevskii:1966, Bertotti:1966, Gunn:1967, Dyer:1974, Dyer:1981, Kasai:1990hd, Kibble:2004tm, Clifton:2009jw, Clifton:2009bp, Clifton:2010fr, Clarkson:2011br, Clifton:2011mt, Bruneton:2012cg, Bruneton:2012ru, Larena:2012vn, Liu:2015bya, Sanghai:2015wia, Fleury:2015rwa, Bentivegna:2016fls, Sanghai:2017yyn, Bentivegna:2018koh, Fleury:2017owg, Fleury:2018cro, Fleury:2018odh}, these studies have been in the context of the geometrical optics approximation, and have also considered the possible effect of inhomogeneities on the average expansion rate \cite{Buchert:2011sx}. We are concerned only with light propagation, and start from the Maxwell equation. In \sec{sec:post-geometric} we derive a post-geometrical approximation that includes curvature and calculate the effect of dark matter clumpiness on redshift and angular diameter distance. We summarise our findings and highlight open questions in \sec{sec:conclusion}. In appendix \ref{sec:area} we calculate the relation between the angular diameter distance and light bundle area expansion rate, and in appendix \ref{sec:duality} we show that distance duality holds in our post-geometrical approximation.

\section{Post-geometrical optics} \label{sec:post-geometric}

\subsection{Post-geometrical approximation}

\paragraph{The equation of motion.}

We start from the action $S=-\int\mathrm{d}^4 x\sqrt{-g}\frac{1}{4}F_{\a\b} F^{\a\b}$ in curved spacetime, where the electromagnetic field strength $F_{\a\b}$ is written in terms of the vector potential $A^\a$ as $F_{\a\b}=\partial_\a A_\b-\partial_\b A_\a$. Varying the action with respect to $A^\a$ gives the equation of motion (we consider an uncharged medium):
\begin{align} \label{eom}
  \Box A^{\a} - \nabla^\a \nabla_\b A^\b - R^\a{}_\b A^\b = 0\ ,
\end{align}
where $\Box\equiv g^{\a\b}\nabla_\a \nabla_\b$ and $R_{\a\b}$ is the Ricci tensor. We take the Lorenz gauge condition
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{Lorenz}
  \nabla_\a A^\a=0 \ ,
\eea
so \re{eom} becomes a curved spacetime generalisation of the wave equation, which depends explicitly on the Ricci tensor.

As in the geometrical optics approximation, we consider the local plane wave form \cite{PhoneBook} (page 570), \cite{Lensing} (page 93)
\begin{align}\label{eq:A_expansion}
	A^{\a} = \sum_{n=0}^\infty \text{Re} ( A^\a_n \epsilon^n e^{i S/\epsilon} ) \ ,
\end{align}
where the amplitude is expanded in a series of the constant parameter $\epsilon\ll1$, which represents the ratio of the wavelength to other relevant scales that are considered large. Unlike in the geometrical optics approximation, we do not assume that the spacetime curvature radius is large compared to the wavelength.

Inserting the expansion \re{eq:A_expansion} into the equation of motion \re{eom} gives, in the Lorenz gauge,
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{eomexp}
    \!\!\!\!\!\!\!\!\! 0 &=& \sum_{n=0}^\infty \text{Re} \left[ \epsilon^n \left( - \epsilon^{-2} k^2 A^\a_n + i \epsilon^{-1} \nabla_\b k^\b A^\a_n + 2 i \epsilon^{-1} k^\b \nabla_\b A^\a_n + \Box A^\a_n - R^\a{}_\b A^\b_n \right) e^{i S/\epsilon} \right] \ ,
\eea
where $k_\a\equiv\partial_\a S$ and $k^2\equiv k_\a k^\a$. In the geometrical optics approximation, the derivatives of $A^\a_n$ are taken to be of the same order of smallness as $A^\a_n$, meaning that amplitude and polarisation change slowly compared to the wavelength. We make the same assumption. In the geometrical optics approximation, $R^\a{}_\b$ is taken to be independent of $\epsilon$. However, if the curvature is large, it will instead be the dominant term in the equation of motion, and will set the scale of the derivatives. We correspondingly assume that $R^\a{}_\b$ is of order $\epsilon^{-2}$. We can still formally smoothly take the limit $R^\a{}_\b\to0$ to recover geometrical optics. However, in between, $R^\a{}_\b$ would contribute at order $\epsilon^{-1}$, so the approximation may not cover regimes of intermediate curvature properly. Expanding the Lorenz gauge condition \re{Lorenz} in the same way, we have
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{Lorenzexp}
    0 &=& \sum_{n=0}^\infty \text{Re} \left[ \epsilon^n \left( \epsilon^{-1} i k_\a A^\a_n + \nabla_\a A^\a_n \right) e^{i S/\epsilon} \right] \ .
\eea

Considering the equation of motion \re{eomexp} at orders $\epsilon^{-2}$ and $\epsilon^{-1}$ gives
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray}
  \label{go1} R_{\a\b} a^\b + k^2 a_\a &=& 0 \\
  \label{go2} R_{\a\b} b^\b + k^2 b_\a &=& a_\a \nabla_\b k^\b + 2 k^\b \nabla_\b a_\a \ ,
\eea
where we have denoted $A^\a_0\equiv a^\a$, $A^\a_1\equiv i b^\a$. Expanding the Lorenz gauge condition \re{Lorenz} to orders $\epsilon^{-1}$ and $\epsilon^0$ gives (for any vectors $v^\a$ and $w^\a$, we denote $v\cdot w\equiv g_{\a\b} v^\a w^\b$)
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{Lorenz2}
  \label{go3} a\cdot k &=& 0 \\
  \label{go4} b \cdot k &=& \nabla_\a a^\a \ .
\eea
This is the same result as in geometrical optics and shows that the light ray is transverse. Contracting \re{go2} with $a^\a$ and using \re{go1}, we obtain
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{pn}
  0 &=& \nabla_\b ( a \cdot a \, k^\b ) \ ,
\eea
so photon number is conserved, as in geometrical optics. In geometrical optics $k^\a$ is null, so it follows from the transversality condition \re{go3} that $a^\a$ is spacelike or null. Here that is not the case, but we assume, as in geometrical optics, that $a^\a$ is spacelike and decompose it as $a^\a=a f^\a$, where $a>0$ and $f\cdot f=1$. Using this decomposition in \re{go2} and applying \re{pn}, we find
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{pol}
  k^\b \nabla_\b f_\a &=& - \frac{1}{2 a} ( R_{\a\b} b^\b + k^2 b_\a ) \ .
\eea
This shows that the polarisation vector is not necessarily parallel transported along the light ray. Using the decomposition of $a^\a$ in \re{go1} and contracting with $f^\b$, we obtain
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{null}
  k^2 &=& - R_{\a\b} f^\a f^\b \ ,
\eea
so in general the light rays are not null. Taking a covariant derivative of \re{null} and using the fact that $\nabla_\b k_\a=\nabla_\a k_\b$ (since $k_\a=\partial_\a S$) we get
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{knk}
  k^\b \nabla_\b k_\a &=& \frac{1}{2} \nabla_\a k^2 = - \frac{1}{2} \nabla_\a ( R_{\a\b} f^\a f^\b ) \ ,
\eea
showing that in general the light rays are not geodesic.

\paragraph{The energy-momentum tensor.}

Consider a unit timelike vector $u^\a$ ($u\cdot u=-1$). Without loss of generality, we can decompose the energy-momentum tensor as
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{Tab}
  T_{\a\b} = ( \rho + p ) u_\a u_\b + p g_{\a\b} + 2 q_{(\a} u_{\b)} + \Pi_{\a\b} \ ,
\eea
where $\rho$ is the energy density, $p$ is the pressure, $q_\a$ is the energy flux, and $\Pi_{\a\b}$ is the anisotropic stress. Both $q_\a$ and $\Pi_{\a\b}$ are orthogonal to $u^\a$, and $\Pi_{\a\b}$ is traceless and symmetric. From the Einstein equation $R_{\a\b}-\frac{1}{2} R g_{\a\b}=M_{{}_{\mathrm{Pl}}}^{-2} T_{\a\b}$, where $R\equiv g^{\a\b} R_{\a\b}$, we get
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{Rab}
  M_{{}_{\mathrm{Pl}}}^2 R_{\a\b} &=& ( \rho + p  ) u_\a u_\b + \frac{1}{2} ( \rho - p ) g_{\a\b} + 2 q_{(\a} u_{\b)} + \Pi_{\a\b} \ .
\eea
Inserting \re{Rab} into \re{go1} and using \re{null}, we obtain
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{fu}
  [ ( \rho + p ) u \cdot f + q \cdot f ] u_\a &=& [ ( \rho + p )  (u\cdot f)^2 + 2 u \cdot f q \cdot f  + \Pi_{\c\d} f^\c f^\d ] f_\a \nonumber \\
  && - u \cdot f q_\a - \Pi_{\a\b} f^\b \ .
\eea
As $u_\a$ is timelike and the right-hand side is a spacelike vector, their coefficients have to vanish separately. From the left-hand side we get $( \rho + p ) u \cdot f + q \cdot f=0$. As there is no reason for the directions of the energy flux and the polarisation to be locally related, the most reasonable way to satisfy this condition is $u\cdot f=0$, $q_\a=0$. The right-hand side condition then reduces to $\Pi^\a{}_\b f^\b=\Pi_{\c\d} f^\c f^\d f^\a$. By the same argument, the most reasonable solution is $\Pi_{\a\b}=0$. These constraints are a limitation of our post-geometrical approximation. When the curvature (and hence the energy-momentum tensor) is directly coupled to light propagation, the local plane wave approximation is not valid if the energy-momentum tensor introduces new privileged spatial directions, because in the leading order plane wave approximation only the propagation direction given by $k^\a$ and the polarisation direction $f^\a$ are relevant. (Although $u^\a$ is timelike, the projection of $f^\a$ orthogonal to it introduces a new spatial direction unless $u\cdot f=0$.) With $u\cdot f=0$, $q_\a=0$, $\Pi_{\a\b}=0$, \re{null} reduces to
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{null2}
  k^2 &=& - \frac{\rho - p}{2 M_{{}_{\mathrm{Pl}}}^2} \ .
\eea
For matter satisfying $\rho\geq p$, the tangent vector $k^\a$ is timelike or null, so the gravitational effect of matter slows down light, analogously to refraction when passing through charged matter.

Applying \re{Tab} to the polarisation propagation equation \re{pol} gives
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{pol2}
  k^\b \nabla_\b f_\a &=& - \frac{1}{2 M_{{}_{\mathrm{Pl}}}^2 a} ( \rho + p ) u \cdot b \, u_\a \ .
\eea
The right-hand side is in general non-zero, so the polarisation vector rotates, i.e.\  there is gravitational circular birefringence, unlike in geometrical optics. This is required by the condition $u\cdot f=0$, because $u^\a$ is in general not parallel transported along the light ray. The magnitude of the birefringence can be calculated by finding $u\cdot b$ from the requirement that the condition $u\cdot f=0$ is preserved along the light ray, although we should then check overall consistency to next order in $\epsilon$.

The propagation equation \re{knk} for $k^\a$ reduces to the simple expression
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{knk2}
  k^\b \nabla_\b k_\a &=& - \frac{1}{4 M_{{}_{\mathrm{Pl}}}^2} \partial_\a ( \rho - p ) \ .
\eea
The results \re{null2} and \re{knk2} allow us to calculate the post-geometrical correction to redshift and angular diameter distance.

\subsection{Redshift} \label{sec:redshift}

\subsubsection{Expression for the redshift}

\paragraph{Dispersion relation.}

We decompose the light ray tangent vector as
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{kdec}
  k^\a &=& E ( u^\a + v e^\a ) \ ,
\eea
where $E=-u \cdot k$, $u \cdot e=0$, $e \cdot e = 1$, and $v\equiv\sqrt{1+k^2/E^2}=\sqrt{1-\frac{\rho-p}{2 M_{{}_{\mathrm{Pl}}}^2 E^2}}$. Correspondingly, the dispersion relation is
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray}
  E^2 = \vec k^2 + \frac{\rho - p}{2 M_{{}_{\mathrm{Pl}}}^2} \ ,
\eea
where $|\vec k|=v E$ is the spatial wavenumber, and as usual the wavelength is $2\pi/|\vec k|$. Light acquires the gravitationally induced mass $M\equiv\sqrt{(\rho-p)/(2M_{{}_{\mathrm{Pl}}}^2)}$, of the order of the curvature scale. For a light ray passing from vacuum to matter the wavenumber has to satisfy $E=|\vec k|>M$, otherwise the wave cannot enter, and the oscillating solution turns into an exponentially decaying mode. This is analogous to the phenomenon that light cannot propagate in materials whose plasma frequency exceeds the light frequency. The refractive index is $n=v^{-1}=(1-M^2/E^2)^{-1/2}$. In a region where $M\neq0$, we can write $k^\a=M w^\a$, where $w^\a$ is a unit timelike vector, and the propagation equation \re{knk} reads $w^\b\nabla_\b w_\a=-M^{-1}\hat\nabla_\a M$, where $\hat\nabla_\a$ is the derivative projected orthogonally to $w^\a$. This can be compared to the equation for an observer comoving with an ideal fluid, $w^\b\nabla_\b w_\a=-(\rho+p)^{-1}\hat\nabla_\a p$.

\paragraph{Redshift.}

The redshift is given by
\begin{align}
	1+z = \frac{E_\text{s}}{E_\text{o}} \ ,
\end{align}
where s refers to source and o to observation. To find the redshift, we proceed as in geometrical optics (see e.g.\  \cite{light_stat1, light_stat2}), and take a derivative with respect to the affine parameter $\lambda$ along the light ray,
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{E1}
  \frac{\mathrm{d} E}{\mathrm{d}\lambda} &=& -k^{\b}\nabla_{\b}\left(u^{\a}k_\a \right) 
 
  = - k^\a k^\b \nabla_\b u_\a - \frac{1}{2} u^\a \nabla_\a k^2 \ ,
\eea
where we have used \re{knk}. We decompose $\nabla_\b u_\a$ as usual (for overviews of the covariant formalism, see \cite{Ellis:1971pg, Tsagas:2007yx})
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{udec}
  \nabla_\b u_\a &=& \frac{1}{3} \theta h_{\a\b} + \sigma_{\a\b} + \omega_{\a\b} -\dot u_\a u_\b \ ,
\eea
where $h_{\a\b}\equiv g_{\a\b}+u_\a u_\b$, $\theta\equiv\nabla_\a u^\a$ is the volume expansion rate, $\sigma_{\a\b}$ is the traceless symmetric shear tensor, $\omega_{\a\b}$ is the antisymmetric vorticity tensor, $\dot u^\a\equiv u^{\b}\nabla_{\b}u^{\a}$ is the acceleration vector, and dot denotes $u^\a \nabla_\a$. The tensors $\sigma_{\a\b}$ and $\omega_{\a\b}$ and the vector $\dot u^\a$ are orthogonal to $u^\a$. Inserting \re{kdec} and \re{udec} into \re{E1}, we get
\begin{align} \label{dE}
	\frac{1}{E^2}\frac{\mathrm{d} E}{\mathrm{d}\lambda} = - \frac{1}{3} v^2 \theta - v ^2 \sigma_{\a\beta} e^{\a}e^{\b} -v e \cdot \dot u - \frac{1}{2 E^2} (k^2)\dot{} \ .
\end{align}
Integrating, we obtain
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{zgen}
	\ln\left(1+z\right) &=& \int_{\lambda_\text{s}}^{\lambda_\text{o}} \mathrm{d}\lambda E^{-1} \left[ \frac{1}{3} v^2 \theta + v ^2 \sigma_{\a\beta} e^{\a}e^{\b} + v e \cdot \dot u + \frac{1}{2 E^2}(k^2)\dot{} \right] \nonumber \\
	&=& \int_{\lambda_\text{s}}^{\lambda_\text{o}} \mathrm{d}\lambda E^{-1} \left[ v^2 \left( \frac{1}{3} \theta + \sigma_{\a\beta} e^{\a}e^{\b} \right) - \frac{v}{\rho+p} e^\a \partial_\a p - \frac{\dot\rho - \dot p}{4 M_{{}_{\mathrm{Pl}}}^2 E^2} \right] \ ,
\eea
where in the second equality we have used \re{null2} and the spatial part of the continuity equation,
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{udot}
  \dot u_\a = -\frac{1}{\rho+p} h_\a{}^\b\partial_\b p \ .
\eea
The result \re{zgen} reduces to the geometrical optics expression when $k^2=0$, i.e.\  $v=1$. In the opposite limit $|\vec k|=0$, i.e.\  $v=0$ and $E^2=M^2=(\rho-p)/(2M_{{}_{\mathrm{Pl}}}^2)$, the first two terms vanish as the curvature-induced mass becomes so large that light cannot propagate inside the matter.

\subsubsection{Size of the curvature correction}

\paragraph{Matter model.}
 
Dark matter energy density consists of widely separated sharp spikes, which we model as
\begin{align}
	\rho = m \sum_n W(r_n) \ ,
\end{align}
where $n$ labels particles, $W$ is a window function that gives the density distribution of a particle (related to its wavefunction), and $r_n\equiv|\vec r-\vec r_n|$ is the proper distance to the centre of particle $n$, located at $\vec r_n$, along a line orthogonal to $u^\a$. We take $W$ to be a Gaussian of width equal to the Compton wavelength $\lambda_\text{c}=2\pi/m$, so
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{rho}
  \rho &=& m \sum_n \frac{1}{(2\pi)^{3/2} \lambda_\text{c}^3} e^{-\frac{r_n^2}{2\lambda_\text{c}^2}} = \sqrt{\frac{2}{9\pi}} \rho_0 \sum_n e^{-\frac{r_n^2}{2\lambda_\text{c}^2}} \ ,
\eea
where $\rho_0\equiv m/(\frac{4}{3}\pi \lambda_\text{c}^3)$ is the energy density averaged over one particle. The distance between dark matter particles is so large that overlap is negligible. We model the region inside the particle using general relativity, so although the pressure is negligible, we need $p\neq0$ for the particle not to collapse, and the pressure gradient may be important. We assume $m\llM_{{}_{\mathrm{Pl}}}$, so the curvature is much smaller than the Planck scale, and we treat spacetime near and inside each particle as perturbed Minkowski space, mostly working to leading order (which often means neglecting the perturbations). We do not consider non-gravitational interactions between dark matter and light.

\paragraph{Correction to redshift.}

We have $\dot\rho=0$ everywhere. This is not in contradiction with the expansion of space and decrease of average density. The time part of the continuity equation, $\dot\rho+\theta\rho=0$, is satisfied because inside particles the expansion rate is zero and outside the density is zero. (The separation is of course blurry as the particles do not have sharp boundaries.) In the redshift integral \re{zgen} the factor $v^2$ in front of $\theta$ is essentially unity whenever $\theta$ is non-zero. The shear is zero inside the particles, because the solution is spherically symmetric and static. The acceleration term integrates to zero, because the spacetime is static inside the particle, so the light gains as much energy going in as it loses coming out. The last term is zero because $\dot\rho=\dot p=0$. There is in principle a correction to the redshift from the missing expansion and slower light travel inside the particles, but this is negligible because the particles occupy only a tiny fraction of the volume. We consider only observers whose time direction is orthogonal to the polarisation, but the redshift of different observers is related to each other by a Lorentz boost as usual. So the curvature correction to the redshift is negligible. We now turn to the angular diameter distance, where the situation is different.

\subsection{Angular diameter distance} \label{sec:distance}

\subsubsection{Expression for the angular diameter distance} \label{sec:DA}

\paragraph{Angular diameter distance and the light bundle.}

In order to determine the angular diameter distance $D_A$, we have to find the area expansion rate of a bundle of light rays orthogonal to $k^\a$ and to a four-velocity. We present the details in appendix \ref{sec:area}. In our post-geometrical approximation (unlike in geometrical optics), the area expansion rate depends on the choice of four-velocity. However, when calculating $D_A$, the velocity-dependent part is subdominant and can be neglected.

It is useful to introduce a tensor that projects orthogonally to $k^\a$ and $u^\a$,
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{htt}
  \tilde{h}_{\a\b} &\equiv& g_{\a\b} + u_\a u_\b - e_\a e_\b \nonumber \\
  &=& g_{\a\b} + \frac{k^2}{E^2+k^2} u_\a u_\b + \frac{2 E}{E^2+k^2} u_{(\a} k_{\b)} - \frac{1}{E^2+k^2} k_\a k_\b \ .
\eea
We decompose $\nabla_\b k_\a$ into parts projected in the direction of and orthogonal to $\tilde{h}_{\a\b}$ as
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{nablakdec}
  \nabla_\b k_\a = \frac{1}{2} \tilde{\theta} \tilde{h}_{\a\b} + \tilde{\sigma}_{\a\b} + P_{\a\b} \ ,
\eea
where the area expansion rate is $\tilde{\theta}\equiv \tilde{h}^{\a\b} \nabla_\b k_\a$, the light shear is $\tilde{\sigma}_{\a\b}=\tilde{h}_\a{}^\c \tilde{h}_\b{}^\d \nabla_\d k_\c-\frac{1}{2} \tilde{h}_{\a\b} \tilde{\theta}$, and $P_{\a\b}\equiv (\d_\a{}^\c \d_\b{}^\d - \tilde{h}_\a{}^\c \tilde{h}_\b{}^\d ) \nabla_\d k_\c$. Unlike in geometrical optics, $P_{\a\b}$ is neither traceless nor orthogonal to $k^\a$. As shown in appendix \ref{sec:area} we have $D_A\propto\exp\left( \frac{1}{2}\int\mathrm{d}\lambda \nabla_\a k^\a \right)=\exp\left[ \frac{1}{2}\int\mathrm{d}\lambda ( \tilde{\theta} + P^\a{}_\a ) \right]$.

\paragraph{Integrating along the light ray.}

In order to integrate $\tilde{\theta}$ over the light ray, we do the same as in geometrical optics, namely take the derivative, and then decompose and integrate it. We have
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{dtheta}
  \frac{\mathrm{d}\tilde{\theta}}{\mathrm{d}\lambda} &=& k^{\b}\nabla_{\b} \nabla_\a k^{\a} - \frac{\mathrm{d} P^\a{}_\a}{\mathrm{d}\lambda} \nonumber \\
  &=& - \frac{1}{2}\tilde{\theta}^2 - 2\tilde{\sigma}^2 - R_{\a\b} k^{\a}k^{\b} + \frac{1}{2}\Box k^2 - P_{\a\b} P^{\a\b} - \frac{\mathrm{d} P^\a{}_\a}{\mathrm{d}\lambda} \ ,
\eea
where $\tilde{\sigma}^2\equiv\frac{1}{2}\tilde{\sigma}_{\a\b}\tilde{\sigma}^{\a\b}$, and we have applied the definition of the Riemann tensor in terms of a commutator of covariant derivatives. In the geometrical optics limit, only the first three terms remain. Expanding $P_{\a\b} P^{\a\b}$ using \re{htt}, we get
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{AA}
  P_{\a\b} P^{\a\b} &=& \frac{1}{v^2 E^2} \left\{ \frac{k^2-E^2}{2 v^2 E^2} [(k^2)\dot{}]^2 + \frac{1}{v^2 E} (k^2)\dot{} (k^2)' - \frac{1}{4 v^2 E^2} [(k^2)']^2 - 2 k^2 \dot k_\a \dot k^\a \right. \nonumber \\
  && \left. - 4 E k_\a' \dot k^\a + 2 k_\a' k^\a{}' - \frac{\dot k \cdot u}{v^2 E} \left[ 2 k^2 (k^2)\dot{} + E (k^2)' \right] - \frac{(k^2)^2}{v^2 E^2} ( \dot k \cdot u )^2 \right\} \nonumber \\
  &=& \frac{1}{v^2 E^2} \left\{ - \frac{1}{16 M_{{}_{\mathrm{Pl}}}^4} [ e^\a\partial_\a (\rho-p)]^2 + \frac{\rho}{M_{{}_{\mathrm{Pl}}}^2} \dot k_\a \dot k^\a + \frac{E}{M_{{}_{\mathrm{Pl}}}^2} \dot k^\a \partial_\a (\rho-p) \right. \nonumber \\
  && + \frac{1}{8 M_{{}_{\mathrm{Pl}}}^4} \partial_\a (\rho-p) \partial^\a (\rho-p) + \frac{E^2}{2 M_{{}_{\mathrm{Pl}}}^2\rho} e^\a\partial_\a p \, e^\b \partial_\b ( \rho - p ) \nonumber \\
  && \left. - \frac{1}{4 M_{{}_{\mathrm{Pl}}}^4} ( e^\a \partial_\a p )^2 \right\} \ ,
\eea
where prime denotes $\frac{\mathrm{d}}{\mathrm{d}\lambda}$, and in the second equality we have used \re{null2}, \re{knk2} and \re{udot}, and taken into account $\dot\rho=\dot p=0$, $p\ll\rho$, and that inside the particles $\dot E=0$. Spatial derivatives come from contracting $\nabla_\b k_\a$ with $k^\a$ (entering from the definition \re{htt} of $\tilde{h}_{\a\b}$), while both $u^\a$ and $k^\a$ contribute to the time derivatives. In geometrical optics, in a spatially flat Friedmann--Lema\^itre--Robertson--Walker (FLRW) universe we have $\nabla_\a k^\a=2E(H+L^{-1})$, where $H$ is the Hubble parameter and $L$ is the proper spatial distance to the source. Because the density changes on a small length scale, the spatial derivative terms are large compared to the geometrical optics terms. However, it is not immediately obvious whether the post-geometrical correction to the angular diameter distance is large, because these terms are non-zero only in a small fraction of the volume.

Let us first show that we can neglect $\int\mathrm{d}\lambda P^\a{}_\a$, and hence the difference between $\nabla_\a k^\a$ and $\tilde{\theta}$ inside the integral. Using \re{htt}, we have
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{P}
  \int\mathrm{d}\lambda P^\a{}_\a &=& - \int\mathrm{d}\lambda \left\{ \frac{k^2}{E^2+k^2} u \cdot \dot k + \frac{1}{E^2 + k^2}\left[ E (k^2)\dot{} - \frac{1}{2} (k^2)' \right] \right\} \nonumber \\
  &=& - \int\mathrm{d}\lambda \left\{ \frac{1}{2 v^2 M_{{}_{\mathrm{Pl}}}^2 E^2}\left[ \dot u \cdot k \rho + \frac{1}{2} E e^\a \partial_\a (\rho-p) \right] \right\} \nonumber \\
  &=& - \int\mathrm{d}\lambda \frac{1}{4 v M_{{}_{\mathrm{Pl}}}^2 E} e^\a \partial_\a ( \rho - 3 p ) \ ,
\eea
where on the second line we have used \re{null2}, taking into account that $\dot\rho=\dot p=0$, $p\ll\rho$, and that inside the particles $\dot E=0$ (as the metric is static, photon energy depends on position, but not explicitly on time), and on the third line we have used \re{udot}. As $\rho$ and $p$ inside each particle depend only on the distance from the centre, the directional derivatives become sums over dot products. For each particle $n$, we have $e^\a \partial_\a = \d_{ij} e^i (r^j-r_n^j)/r_n \partial_{r_n}$. The dot product gives the cosine of the angle between the light ray and a line from the centre of particle $n$. When we average over many rays within a tube of area $A$, we integrate over the entry angle, with a distribution that is flat in the cosine of the entry angle. (Summing over a large number of particles for a single ray would give the same distribution.) Assuming the change in light ray direction is small or symmetric around the centre of each particle, we get zero.

As the contribution of $P^\a{}_\a$ is negligible, we have from \re{dtheta} and \re{AA}
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{thetaint}
  \int\mathrm{d}\lambda \nabla_\a k^\a &=& \int\mathrm{d}\lambda \tilde{\theta} \nonumber \\
  &=& \int\mathrm{d}\lambda \int\mathrm{d}\lambda' \left( - \frac{1}{2} \tilde{\theta}^2 - 2 \tilde{\sigma}^2 -R_{\a\b} k^{\a}k^{\b} + \frac{1}{2}\Box k^2 -P_{\a\b} P^{\a\b} \right) \nonumber \\
  &=& \int\mathrm{d}\lambda \int\mathrm{d}\lambda' \Bigg\{ - \frac{1}{2} \tilde{\theta}^2 - 2 \tilde{\sigma}^2 \underbrace{ - E^2 \frac{\rho + p}{M_{{}_{\mathrm{Pl}}}^2} \left[ 1 - \frac{(\rho-p)^2}{4M_{{}_{\mathrm{Pl}}}^2 E^2 (\rho+p)} \right]}_{\text{term 1}} \nonumber \\
  && \underbrace{- \frac{1}{4 M_{{}_{\mathrm{Pl}}}^2} \Box ( \rho - p )}_{\text{term 2}} \underbrace{- P_{\a\b} P^{\a\b}}_{\text{term 3}} \Bigg\} \ ,
\eea
where we have used \re{Rab} and \re{null2}. The first two terms are the same as in geometrical optics. Let us consider the three post-geometrical correction terms in turn.

\subsubsection{Size of the curvature correction} \label{sec:corr}

\paragraph{Correction term 1.}

Over one wavelength, a light ray typically passes through no particles (except for long light wavelength and small dark matter mass). However, measurements are made with finite beam size, and the effect of dark matter is smeared over the width of the beam. Let us consider an interval $\Delta\lambda$ that corresponds to a time interval $\Delta t$ that is small compared to the timescale of cosmological evolution but large compared to the distance between dark matter particles, and average over the beam width $A$, approximating spacetime in this tube as perturbed Minkowski space. The first post-geometrical correction term in \re{thetaint} involves the integral
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{rho2}
  \int_{\Delta\lambda} \mathrm{d}\lambda E^2 \frac{\rho}{M_{{}_{\mathrm{Pl}}}^2} \left( 1 - \frac{\rho}{4M_{{}_{\mathrm{Pl}}}^2 E^2} \right) &=& \int_{\Delta t}\mathrm{d} t E \frac{\rho}{M_{{}_{\mathrm{Pl}}}^2} \left(1 - \frac{\rho}{4 M_{{}_{\mathrm{Pl}}}^2 E^2} \right) \nonumber \\
  &=& \frac{1}{A} \int_{V} \mathrm{d}^3 r E \frac{\rho}{M_{{}_{\mathrm{Pl}}}^2} \left(1 - \frac{\rho}{4 M_{{}_{\mathrm{Pl}}}^2 E^2} \right) \nonumber \\
  &=& \frac{1}{A} \sum_n \int_{V} \mathrm{d}^3 r E \frac{m}{M_{{}_{\mathrm{Pl}}}^2} W(r_n) \left[ 1 - \frac{m}{4M_{{}_{\mathrm{Pl}}}^2 E^2} W(r_n) \right] \nonumber \\
  &=& \frac{\Delta t}{V} E N \frac{m}{M_{{}_{\mathrm{Pl}}}^2} \left( 1 - \frac{\rho_0}{24 \sqrt{\pi} M_{{}_{\mathrm{Pl}}}^2 E^2} \right) \nonumber \\
  &=& \int_{\Delta\lambda} \mathrm{d}\lambda E^2 \frac{\av{\rho}}{M_{{}_{\mathrm{Pl}}}^2} \left( 1 - \frac{\rho_0}{24 \sqrt{\pi} M_{{}_{\mathrm{Pl}}}^2 E^2} \right) \ ,
\eea
where we have used $\mathrm{d}\lambda E = \mathrm{d} t$, taken into account $p\ll\rho$, inserted the Gaussian density distribution \re{rho}, and neglected the small variation in $E$. Here $V=A \Delta t$ is the integration volume, $N$ is the number of particles in $V$ and $\av{\rho}=mN/V$ is the average energy density of dark matter. We have assumed that the situation is statistically homogeneous and isotropic and evolves little on the timescales it takes for light to cross the interval, so that we can replace the integral along the light ray with a spatial integral \cite{light_stat1, light_stat2, Lavinto:2013exa}. When converting back into an integral, we have smeared the contribution from the small fraction of volume inside the particles into an average contribution over the light path, again assuming that changes are slow. Note how the small average density emerges from the large spiky local density for the first term, which is the same as in geometrical optics. However, there is no such suppression for the second factor of density, which gives a correction of the order $M^2/E^2$.

\paragraph{Correction term 2.}

The second post-geometrical term in \re{thetaint} involves the integral
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{nabla2rho}
  \frac{1}{4 M_{{}_{\mathrm{Pl}}}^2} \int_{\Delta\lambda} \mathrm{d}\lambda \Box(\rho-p) &=& \frac{1}{4 M_{{}_{\mathrm{Pl}}}^2} \int\mathrm{d}\lambda \d^{ij} \partial_i \partial_j (\rho - p ) \nonumber \\
  &=& \frac{A^{-1} m}{4 M_{{}_{\mathrm{Pl}}}^2} \sum_n \int_V \mathrm{d}^3 r E^{-1} \d^{ij} \partial_i \partial_j W(r_n) \nonumber \\
  &=& \frac{A^{-1} N m}{4 M_{{}_{\mathrm{Pl}}}^2} \sum_n \int_V \mathrm{d}^3 r E^{-1} \d^{ij} \partial_i \partial_j W(r) \nonumber \\
  &=& \frac{A^{-1} N m}{4 M_{{}_{\mathrm{Pl}}}^2} \int_V \mathrm{d}^3 r E^{-2} \partial_r W(r) \partial_r E \ ,
\eea
where we have taken into account $p\ll\rho$. As spacetime inside each particle is static, $E\propto|g_{00}|^{-1/2}$. With $u^\a=|g_{00}|^{-1/2}$, we have $\dot u_r=\frac{1}{2}\partial_r\ln|g_{00}|\approx\frac{1}{2}\partial_r|g_{00}|$, which according to \re{udot} is equal to $-\rho^{-1}\partial_r p$, so $\partial_r E=E \rho^{-1} \partial_r p$. The pressure gradient is given by the Tolman--Oppenheimer--Volkov equation (dropping the subscript $n$),
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{pder}
  \frac{\partial p}{\partial r} &=& - G_{\mathrm{N}} \frac{ (\rho+p) ( \tilde m + 4\pi r^3 p) }{r (r-2G_{\mathrm{N}} \tilde m)} = - \frac{1}{8\piM_{{}_{\mathrm{Pl}}}^2} \frac{ \rho \tilde m }{r^2} \ ,
\eea
where $\tilde m(r)\equiv\int_0^r \mathrm{d} r' 4\pi r'^2 \rho$, and in the second equality we have taken into account that for a Gaussian density distribution and $m\llM_{{}_{\mathrm{Pl}}}$, the Newtonian limit is a good approximation. Therefore
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{Eder}
  \frac{\partial E}{\partial r} &=& - \frac{E}{8\piM_{{}_{\mathrm{Pl}}}^2} \frac{\tilde m}{r^2} \ ,
\eea
Applying \re{Eder}, the integral \re{nabla2rho} becomes
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{nabla2rho2}
  \int_{\Delta\lambda} \mathrm{d}\lambda \frac{1}{4 M_{{}_{\mathrm{Pl}}}^2} \Box(\rho-p) &=& - \frac{A^{-1} N m E^{-1}}{8M_{{}_{\mathrm{Pl}}}^4} \int_0^\infty \mathrm{d} r \tilde m(r) \partial_r W(r) \nonumber \\
  &=& \frac{A^{-1} N m E^{-1}}{8M_{{}_{\mathrm{Pl}}}^4} \int_0^\infty \mathrm{d} r \partial_r \tilde m(r) W(r) \nonumber \\
  &=& \int_{\Delta\lambda} \mathrm{d}\lambda \frac{\rho^2}{8M_{{}_{\mathrm{Pl}}}^2} \ ,
\eea
where we have neglected the subleading variation of $E$. This term cancels half of the correction term in \re{rho2}.

\paragraph{Correction term 3.}

The third post-geometrical term in \re{thetaint} involves an integral over \re{AA}. It contains $\dot k^\a$, which we can write as
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{kdot}
  \dot k_\a &=& u^\b \nabla_\b k_\a \nonumber \\
  &=& - \partial_\a E - E v \left( \frac{1}{3} \theta e_\a + \sigma_{\a\b} e^\b + \omega_{\a\b} e^\b - \dot u \cdot e \, u_\a \right) \nonumber \\
  &=& - \partial_\a E - E v \rho^{-1} e^\b \partial_\b p \, u_\a \ ,
\eea
where we have used $\nabla_\b k_\a=\nabla_\a k_\b$, $E=-u \cdot k$, the decomposition \re{udec}, and on the third line we have used \re{udot}, and taken into account $p\ll\rho$ and that inside the particles $\theta=0$ and $\sigma_{\a\b}=\omega_{\a\b}=0$. Using \re{kdot} in \re{AA}, the integral over $P_{\a\b}P^{\a\b}$ reads
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{AAint}
   && \int\mathrm{d}\lambda \frac{1}{M_{{}_{\mathrm{Pl}}}^4 E^2 v^2} \left[ - \frac{1}{16} ( e^\a\partial_\a \rho)^2 + \frac{1}{8} \left( 1+ \frac{4 M_{{}_{\mathrm{Pl}}}^2 E^2}{\rho} \right) e^\a\partial_\a \rho \, e^\b \partial_\b p \right. \nonumber \\
  && - \frac{1}{16} \left( 5 + 24 \frac{M_{{}_{\mathrm{Pl}}}^2 E^2}{\rho} \right) ( e^\a \partial_\a p )^2 + \frac{1}{8} \partial_\a (\rho-p) \partial^\a (\rho-p) \nonumber \\
  && \left. - M_{{}_{\mathrm{Pl}}}^2 E \partial_\a E \partial^\a (\rho-p) + M_{{}_{\mathrm{Pl}}}^2 \rho \partial_\a E \partial^\a E \right] \ .
\eea
The directional derivatives in the first three terms again become sums over dot products. For each particle $n$, we have $e^\a \partial_\a = \d_{ij} e^i (r^j-r_n^j)/r_n \partial_{r_n}$. In the sum over particles, the only pairs that give non-zero results are those where the angle is the same, so $\sum_{m,n} \d_{ij} e^i (r^j-r_m^j)/r_m \d_{kl} e^k (r^l-r_n^l)/r_n=\frac{1}{3}\sum_n$, again assuming the change in light ray direction is small or symmetric around the centre of each particle. We get 
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{AAint2}
   && \int\mathrm{d}\lambda \sum_n \frac{1}{M_{{}_{\mathrm{Pl}}}^4 E^2 v^2} \left[ \frac{5}{48} ( \partial_{r_n} \rho)^2 + \frac{1}{24} \left( - 5 + \frac{4 M_{{}_{\mathrm{Pl}}}^2 E^2}{\rho} \right) \partial_{r_n} \rho \, \partial_{r_n} p \right. \nonumber \\
  && \left. + \frac{1}{48} \left( 1 - 24 \frac{M_{{}_{\mathrm{Pl}}}^2 E^2}{\rho} \right) ( \partial_{r_n} p )^2 - M_{{}_{\mathrm{Pl}}}^2 E \partial_{r_n} E \partial_{r_n} (\rho-p) + M_{{}_{\mathrm{Pl}}}^2 \rho ( \partial_{r_n} E)^2 \right] \ .
\eea
Using $\partial_{r_n} E=E \rho^{-1} \partial_{r_n} p$, \re{AAint2} becomes
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{AAint3}
   && \int\mathrm{d}\lambda \sum_n \frac{1}{M_{{}_{\mathrm{Pl}}}^4 E^2 v^2} \left[ \frac{5}{48} ( \partial_{r_n} \rho)^2 - \frac{1}{24} \left( 5 + \frac{20 M_{{}_{\mathrm{Pl}}}^2 E^2}{\rho} \right) \partial_{r_n} \rho \, \partial_{r_n} p \right. \nonumber \\
  && \left. + \frac{1}{48} \left( 1 + \frac{72 M_{{}_{\mathrm{Pl}}}^2 E^2}{\rho} \right) ( \partial_{r_n} p )^2 \right] \ ,
\eea
where the pressure gradient is given in \re{pder}. We approximate the common prefactor $v^2$ as unity. If $v$ is appreciably different from unity, it will vary inside the particle, but this is relevant only when the light wave nearly has too little energy to cross the particle. Taking into account the variation would not change the order of magnitude of the corrections, but the contribution from the variation and the leading term would contribute to the subleading terms. Somewhat inconsistently, we do take into account the variation of $E$ in the $(\partial_{r_n} \rho)^2$ term, as it can be readily done and affects the subleading terms. We expand $E^{-2}\approx E_0^{-2}(1-2rE_0^{-1}\partial_r E)$, where $E_0$ is the energy at $r_n=0$, giving the extra term $-\frac{5}{24} \rho^{-1} ( \partial_{r_n} \rho)^2r\partial_r p$ in the integral \re{AAint3} compared to the case when the variation in $E$ is neglected. (The contribution of this term turns out to be $-2$ times the contribution of the term $-\frac{5}{24} \partial_{r_n} \rho \partial_{r_n} p$.) The terms in \re{AAint2} reduce to volume integrals involving factors of $\rho$, $\partial_r\rho$, $\tilde m$ and $r$. Using the Gaussian density distribution \re{rho} and the pressure gradient \re{pder}, we have
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{numint}
  \int\mathrm{d}^3 r (\partial_r \rho)^2 &=& \frac{m}{16\pi^{5/2}} m^2 \rho_0 \nonumber \\
  \int\mathrm{d}^3 r \partial_r \rho \partial_r p &=& \frac{m}{54\sqrt{3}\pi} \frac{\rho_0^2}{M_{{}_{\mathrm{Pl}}}^2} \nonumber \\
  \int\mathrm{d}^3 r \rho^{-1} \partial_r \rho \partial_r p &=& \frac{m}{12\sqrt{\pi}} \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2} \nonumber \\
  \int\mathrm{d}^3 r \rho^{-1} (\partial_r \rho)^2 r \partial_r p &=& - \frac{m}{27\sqrt{3}\pi} \frac{\rho_0^2}{M_{{}_{\mathrm{Pl}}}^2} \nonumber \\
  \int\mathrm{d}^3 r (\partial_r p)^2 &\approx& 3\times10^{-9} m \frac{m^4}{M_{{}_{\mathrm{Pl}}}^4} m^2 \rho_0 \nonumber \\
  \int\mathrm{d}^3 r \rho^{-1} (\partial_r p)^2 &\approx& 4\times10^{-5} m \frac{m^2}{M_{{}_{\mathrm{Pl}}}^2} \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2} \ .
\eea
The precise numerical values depend on the form of the density distribution, but the orders of magnitude are simple to understand. Over the integral, the radial derivative gives a factor of $\lambda_\text{c}^{-1}$ as that is the scale of variation of the density, so $\int\mathrm{d}^3 r (\partial_r \rho)^2\sim\lambda_\text{c}^{-2} m \rho_0\sim10^{-2} m^3 \rho_0$. The pressure gradient is suppressed by the factor $m^2/M_{{}_{\mathrm{Pl}}}^2$ compared to the density gradient, so terms proportional to $(\partial_r p)^2$ are negligible. For a steeper density profile, the integrals involving $\partial_r\rho$ would be larger; for a step function the integral over $(\partial_r\rho)^2$ diverges. In the integral over the light rays, one factor of $m$ combines with $N/V$ to become the average density $\av{\rho}$ as before. Overall, \re{AAint3} becomes
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{AAfinal}
\int\mathrm{d}\lambda E^2 \frac{\av{\rho}}{M_{{}_{\mathrm{Pl}}}^2} \left[ \frac{5}{768\pi^{5/2}} \frac{m^2}{E^2} \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2 E^2} - \frac{5}{72\sqrt{\pi}} \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2 E^2} + \frac{5}{1296\sqrt{3}\pi} \left( \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2 E^2} \right)^2 \right] \ .
\eea

\paragraph{The total correction.}

In summary, gathering the terms \re{rho2}, \re{nabla2rho2} and \re{AAfinal}, the integral \re{thetaint} is
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{thetaint2}
  \int\mathrm{d}\lambda \tilde{\theta} &=& \int\mathrm{d}\lambda \int\mathrm{d}\lambda' \left\{ - \frac{1}{2}\tilde{\theta}^2 - 2\tilde{\sigma}^2 - E^2 \frac{\av{\rho}}{M_{{}_{\mathrm{Pl}}}^2} \left[ 1 + \frac{5}{8192\pi^{13/2}} \frac{m^6}{M_{{}_{\mathrm{Pl}}}^2 E^4} \right. \right. \nonumber \\
  && \left. \left. - \frac{13}{144\sqrt{\pi}} \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2 E^2} + \frac{5}{1296\sqrt{3}\pi} \left( \frac{\rho_0}{M_{{}_{\mathrm{Pl}}}^2 E^2} \right)^2 \right] \right\} \ .
\eea
The maximum density is $\sqrt{2/(9\pi)}\rho_0$, so the maximum induced mass squared is $M^2=\rho_0/(\sqrt{18\pi}M_{{}_{\mathrm{Pl}}}^2)$. If $E^2$ is smaller than this, the wave cannot propagate through the particle. For this limiting energy, the last two terms in \re{thetaint2} reach their maximum absolute value 0.4 and 0.04, respectively. Although the numerical prefactors are sensitive to the density profile, the general picture is that these correction terms are small until the energy falls close to the curvature scale, when our approximation is not necessarily valid any more. Once the curvature-induced mass exceeds the energy, light cannot propagate inside the particles, and our plane wave approximation is certainly not valid.

These last two correction terms determined by the curvature scale are in line with the naive estimates discussed in the introduction. In contrast, the first correction term in \re{thetaint2} is given not by the curvature, but by the square of the derivative of the curvature, and is hence enhanced by the factor $m^2/E^2$. This term can give a large correction to the angular diameter distance even when the curvature is much smaller than the photon energy, and the dispersion relation of light is close to the vacuum case. The relative correction is unity at $E=0.02 (m^3/M_{{}_{\mathrm{Pl}}})^{1/2}$. (Note that in this limit, the derivative of the curvature remains small relative to $E^3$.) Expressed another way, for photons of energy $E$, this correction is significant if $m\gtrsim10 (M_{{}_{\mathrm{Pl}}} E^2)^{1/3}$. For the longest radio waves observable through the Earth's ionosphere, this gives $m\gtrsim100$ keV, and for CMB photons we get $m\gtrsim100$ MeV. If correct, this is a strong result: we can rule out heavy dark matter from the absence of a signal in the CMB. Although the precise limit on the mass changes with the density distribution inside the particle, the energy dependence does not. If we assumed that light can pass through electrons and nuclei in the same way as through dark matter, they would have negligible effect on the CMB, but not necessarily on radio waves.

\paragraph{From $\tilde{\theta}$ to $D_A$.}

In terms of the angular diameter distance $D_A$, \re{thetaint2} corresponds to (see e.g.\  \cite{light_stat1, light_stat2}), dropping those correction terms that are important only near the curvature limit, and neglecting the null shear (which is small in the real universe),
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{DA}
   D_A'' &=& - 4\piG_{\mathrm{N}} \av{\rho} \left[ 1 + \a (E/E_\text{o})^{-4} \right] D_A \ ,
\eea
where $\a\equiv\frac{5}{8192\pi^{13/2}} \frac{m^6}{M_{{}_{\mathrm{Pl}}}^2 E_\text{o}^4}$. In contrast to geometrical optics, the angular diameter distance depends on energy. Switching from the affine parameter $\lambda$ to the redshift $z$ using $\frac{\mathrm{d}}{\mathrm{d}\lambda}=-(1+z)H(z)E\frac{\mathrm{d}}{\mathrm{d} z}$, where $H(z)$ is the Hubble parameter, we obtain
\begin{eqnarray}} \newcommand{\eea}{\end{eqnarray} \label{DAz}
  H(z) \frac{\mathrm{d}}{\mathrm{d} z} \left[ (1+z)^2 H(z) \frac{\mathrm{d} D_A}{\mathrm{d} z} \right] &=& - 4\piG_{\mathrm{N}} \av{\rho} \left[ 1 + \a (1+z)^{-4} \right] D_A \ .
\eea
Because $\av{\rho}\propto (1+z)^3$, the new contribution to the source term in the angular diameter distance grows like $(1+z)^{-1}$, and becomes important at late times. The redshift dependence is the same as for dark energy with equation of state $p_\text{de}=-\frac{4}{3}\rho_\text{de}$. However, the sign $\a>0$ in \re{DAz} corresponds to $\rho_\text{de}<0$. Also, it is impossible to get a correction term like this from a matter source in a FLRW universe, because matter would also modify the mapping between $\lambda$ and $z$ given by $H(z)$, not just the source term, unlike here. This is the reverse of corrections from backreaction, which change the $\lambda-z$ mapping, but not the source term \cite{light_stat1, light_stat2, Bull:2012zx, Lavinto:2013exa}. A value $\a\gtrsim1$ on a given wavelength would significantly affect $D_A$, and measurements of distances using light with different wavelength would give discrepant results when interpreted without the correction term. As $\a>0$, the correction leads to a shorter distance, so this effect cannot mimic a cosmological constant. Nor can it explain the tension between the determination of $H_0$ (i.e.\  the inverse of the normalisation of the distance) from the CMB and local measurements of $H_0$ (which are only weakly dependent on a correction like this) \cite{Verde:2019ivm}, because the inferred value of $H_0^{-1}$ from the CMB is larger. (If $\a$ had the opposite sign, it could have such an effect, and would for the right value even change $D_A$ in a manner similar to the cosmological constant.)

We have considered dark matter only, but $\av{\rho}$ usually also includes baryonic matter. (As the full term is $\av{\rho}+\av{p}$, the cosmological constant does not contribute.) If we assume that light cannot propagate through nuclei due to electromagnetic interactions, they do not contribute to the emergence of the average density from the spiky local density. This leads to reduced density along the line of sight, and hence increased distance. In the $\Lambda$CDM model the Planck mean parameters are $\Omega_\text{m}=\Omega_\text{dm}+\Omega_\text{b}=0.315$, $\Omega_\Lambda=1-\Omega_\text{m}$, and $\Omega_\text{b}=0.0224/h^2$, where $h=H_0/(100 \text{km/s/Mpc)}$ \cite{Aghanim:2018eyx}. We can roughly estimate the effect of not propagating through baryons by dropping the baryonic contribution to the angular diameter distance integral $D_A= H_0^{-1} (1+z)^{-1} \int_0^z \frac{\mathrm{d} z'}{\sqrt{ \Omega_\Lambda + (\Omega_\text{dm}+\Omega_\text{b}) (1+z')^3 }}$. As the denominator in the integral is then smaller than in the usual $\Lambda$CDM case, the distance is longer. If we fit the $\Lambda$CDM model (with baryons included) to the CMB, the inferred value of $H_0$ will be smaller to compensate. If we adopt the local measurement $H_0=73.2\pm1.3$ km/s/Mpc \cite{Riess:2020fzl} as the real value, the inferred value is $69.0\pm1.2$ km/s/Mpc, closer to the value $67.4\pm0.5$ km/s/Mpc reported by Planck. Other determinations of $H_0$ involving the integral over the light ray would also be affected differently. In any case, this estimate is rather naive, and the issue of light propagation inside charged particles should be considered in more detail.

Possible violation of the distance duality between the angular diameter distance and the luminosity distance could also be important, as the duality is a key input in drawing conclusions about the compatibility of CMB and baryon acoustic oscillation observations and observations of type Ia supernovae. In appendix \ref{sec:duality} we show that although the deviation equation of two light rays receives a large post-geometrical correction, the reciprocity relation and hence distance duality remains unaffected.

\subsection{Discussion} \label{sec:disc}

\paragraph{Assumptions about light and matter.}

Let us consider the approximations we have made concerning light propagation and modelling of matter. We have assumed that the scales over which amplitude and polarisation of the light wave change are large compared to the wavelength, while spacetime curvature can be large, and that the electromagnetic field has the form of a local plane wave. This is only consistent if the energy flux and anisotropic stress can be neglected and the polarisation is orthogonal to the observer four-velocity. The first two conditions are reasonable, the last is more problematic, although differences in source and observer velocity are small in practice.

In the case of the amplitude, if the induced mass is close to the photon energy, $M\sim E$, the dispersion relation is strongly modified, so we could expect that the light wavefront does not remain straight, and the wave becomes choppy on small scales as it breaks on the dark matter particles. For general values of $M/E$, we can estimate the consistency of assuming that the amplitude changes slowly from \re{pn}, which gives $a'= - \frac{1}{2} a \nabla_\a k^\a$. (Note that this equation remains valid even outside our approximation as long as photon number is conserved.) So $a\propto 1/D_A$ changes slowly. As for the polarisation, contracting \re{pol2} with $k^\a$ and applying the decomposition \re{udec} gives the estimate $f_\a'\sim (E m^3)/M_{{}_{\mathrm{Pl}}}^2\sim (E M^2)/m$. The second derivative of the polarisation in the fourth term in \re{eomexp} is $\sim E^{-2} f_\a''\sim E^{-1} m f_\a'\sim M^2$. So in the observationally allowed regime where the leading correction to $D_A$ is below unity it is formally much smaller than the wavenumber $|\vec k|\approx E$, which sets the expansion scale. Of course both $k^2$ and the curvature term are smaller than $E^2$ as well (in the geometrical optics approximation, the leading term $k^2$ is zero), so this issue should be checked in detail.

In any case, these order of magnitude estimates can at best support consistency, not establish correctness. A related issue is that our approximation may not treat regimes of intermediate curvature between the curvature spikes and empty space correctly. If the polarisation and/or the amplitude also vary rapidly, this could naively be expected to push the situation even further away from the range of validity geometrical optics approximation, though correction terms could cancel. The validity of these approximations related to light propagation can in principle be straightforwardly checked with numerical solutions of the Maxwell equation in curved spacetime, given a matter distribution. This brings us to the assumptions related to the modelling of matter.

We have assumed that dark matter particle mass can be treated as being spread in a classical Gaussian distribution, which is related to spacetime curvature via the Einstein equation. While the particles are inherently quantum mechanical, the spatial distances and curvatures we consider are far below the Planck scale, so this should be a reasonable approximation. If that were not the case, then the quantum mechanical nature of the particles and/or their gravity would play a role, which would push the situation even further from the range of validity of the usual treatment. A key assumption is that the width of the Gaussian is given by the Compton wavelength, and the numerical limits we obtain depend on this choice (and on the overall shape of the mass distribution). If the width of the Gaussian is $l$ instead of $\lambda_\text{c}$, the density changes by the factor $(\lambda_\text{c}/l)^3$, and the leading correction to the angular diameter distance picks an extra factor of $(\lambda_\text{c}/l)^2$ from two derivatives. For example, if we used the de Broglie wavelength $2\pi/(m v_\text{dm})$ instead, where $v_\text{dm}$ dark matter velocity, the mass related to the curvature limit would increase by the factor $v_\text{dm}^{-3/4}\sim10^2$ (taking $v_\text{dm}\sim10^{-3}$), to $m\sim10^6$ GeV. The limiting mass for the leading correction to the angular diameter distance would increase by the factor $v_\text{dm}^{-5/6}\sim10^2$, to $\sim 10$ GeV.

The mass distribution is determined by the wave function (or more generally the density matrix) of dark matter particles. The width of a free particle wavepacket grows linearly in time, and would quickly grow to make the curvature spikes negligible. However, cosmological dark matter particles constantly interact via gravity --and many candidates also via other interactions-- with light, neutrinos, baryonic matter and each other. Dark matter particles decohere via these interactions \cite{Joos:1984uk, Tegmark:1993yn, Allali:2020ttz, Allali:2020shm, Allali:2021puy}, and should also localise, as otherwise their wavefunctions would have spread over billions of light years by today. However, if the particles remain delocalised over scales much larger than the inverse mass, the limits on the mass we have found become too weak to make any difference.

\section{Conclusions} \label{sec:conclusion}

\paragraph{Limits on dark matter from light propagation.}

In general relativity, light travel on null geodesics is an approximation valid in the geometrical optics limit when wavelength is small compared to the spacetime curvature radius and the scale over which the amplitude and the polarisation change. We have introduced a post-geometrical approximation that includes spacetime curvature, and derived the expressions for the redshift and angular diameter distance. We find that inside matter photons acquire a gravity-induced mass $M$ of the order of the curvature scale. Light paths are not null, and the gradient of the matter density (and pressure) pushes photons off the geodesic path. If $M>E$, light cannot propagate inside the particles. We approximate dark matter particles of mass $m$ as Gaussian wavepackets with width equal to the Compton wavelength $\lambda_\text{c}$. We find that the correction to the redshift is negligible. 

In the angular diameter distance $D_A$, the integral of the spiky density over the light ray reproduces the usual effect of the average density. However, we also find a post-geometrical correction that is enhanced by the factor $m^2/E^2$ over the naive expectation that the effects are suppressed by $M^2/E^2$. The reason is that not only the density but also its derivatives enter into $D_A$, and the density changes on the scale $\lambda_\text{c}\sim1/m$ inside the particle. The origin of the derivatives is the geodesic deviation vector $k^\b\nabla_\b k^\a$, which vanishes in geometrical optics, but here is proportional to the gradient of the density. The correction to the angular diameter distance is of order unity or larger for $m\gtrsim10(M_{{}_{\mathrm{Pl}}} E^2)^{1/3}$. For CMB photons $E\sim(1+z)10^{-3}$ eV, and the angular diameter distance to the CMB is strongly affected for $m\gtrsim 100$ MeV. The curvature-suppressed corrections would give the much weaker limit $m\gtrsim10\sqrt{M_{{}_{\mathrm{Pl}}} E}\sim10^4$ GeV.

As the correction term to the angular diameter distance depends on photon energy, observations made at different wavelengths could yield discrepant results if this term is not taken into account. The effect only makes $D_A$ smaller, so it cannot mimic the effect of a cosmological constant nor explain the $H_0$ tension between the CMB and other probes, and can only make the fit to observations worse.

We find that polarisation undergoes gravity-driven rotation which could provide an interesting observational signature. It would be interesting to look at the effect not only on the angular diameter distance (i.e.\  beam convergence), but also on lensing (i.e.\  beam deformation), studied for extended beams in the geometrical optics approximation in \cite{Fleury:2017owg, Fleury:2018cro, Fleury:2018odh}. In addition to considering in more detail the case $M\sim E$, when our approximation may not be valid, it would be interesting to go beyond to the case $M>E$, when light cannot propagate inside dark matter, which falls outside the local plane wave approximation.

If the results hold, dark matter masses $m\gtrsim100$ MeV are ruled out by CMB observations, and observations of $D_A$ using longer wavelengths, such as the 21 cm line, offer the prospect of tightening this bound by a few orders of magnitude. This may offer a novel way to probe the microscopic nature of dark matter via observations of light. However, a more detailed study of the validity of the assumptions we have made about light propagation and modelling of dark matter is required before drawing such conclusions.

\acknowledgments

We thank Pierre Fleury for helpful correspondence and the anonymous referee for useful criticism. Sofie Marie Koksbang is funded by the Carlsberg Foundation.

