# The effect of dark matter discreteness on light propagation

Sofie Marie Koksbang  and Syksy Räsänen

###### Abstract

Light propagation in cosmology is usually studied in the geometrical optics
approximation which requires the spacetime curvature to be much smaller than
the light wavenumber. However, for non-fuzzy particle dark matter the
curvature is concentrated in widely separated spikes at particle location. If
the particle mass is localised within a Compton wavelength, then for masses
$\gtrsim 10^{4}$ GeV the curvature is larger than the energy of CMB photons.

We consider a post-geometrical optics approximation that includes curvature.
Photons gain a gravity-induced mass when travelling through dark matter, and
light paths are not null nor geodesic. We find that the correction to the
redshift is negligible. For the angular diameter distance, we show how the
small average density emerges from the large local spikes when integrating
along the light ray. We find that there can be a large correction to the
angular diameter distance even for photon energies much larger than the
curvature. This may allow to set a strong upper limit on the mass of dark
matter particles. We discuss open issues related to the validity of our
approximations.

HIP-2021-25/TH

##  1 Introduction

##### 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 [1] (page
570), [2] (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. [3, 4, 5, 6, 7] and references therein111We focus on light; for
discussion of geometrical optics and wave optics for gravitational waves, see
[8, 9, 10, 11, 12, 13].). 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 [14, 15, 16, 17, 18]. Expanding in
the short wavelength approximation beyond the leading terms [19, 20], the
covariant derivatives in the equation of motion lead to spin-dependent light
propagation, known as the gravitational Hall effect [21, 22, 23]. 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 $\sim 10$ m and energy $E\sim 10^{-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\sim 10$ 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\sim 10^{-3}$ GeV4, which gives
$\rho/(E^{2}M_{{}_{\mathrm{Pl}}}^{2})\sim 10^{-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\sim 10^{5}$ GeV, and
the corresponding dark matter mass is $m\gtrsim 10^{12}$ GeV. Cosmic microwave
background (CMB) photons at redshift $z$ have energy $\sim 10^{-3}(1+z)$ eV,
so for $m\gtrsim 10^{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
[24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45, 46, 47, 48, 49], 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 [50]. We are concerned only
with light propagation, and start from the Maxwell equation. In section 2 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 section 3.
In appendix A we calculate the relation between the angular diameter distance
and light bundle area expansion rate, and in appendix B we show that distance
duality holds in our post-geometrical approximation.

##  2 Post-geometrical optics

###  2.1 Post-geometrical approximation

##### The equation of motion.

We start from the action
$S=-\int\mathrm{d}^{4}x\sqrt{-g}\frac{1}{4}F_{\alpha\beta}F^{\alpha\beta}$ in
curved spacetime, where the electromagnetic field strength $F_{\alpha\beta}$
is written in terms of the vector potential $A^{\alpha}$ as
$F_{\alpha\beta}=\partial_{\alpha}A_{\beta}-\partial_{\beta}A_{\alpha}$.
Varying the action with respect to $A^{\alpha}$ gives the equation of motion
(we consider an uncharged medium):

$\displaystyle\Box
A^{\alpha}-\nabla^{\alpha}\nabla_{\beta}A^{\beta}-R^{\alpha}{}_{\beta}A^{\beta}=0\
,$ (2.1)

where $\Box\equiv g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}$ and
$R_{\alpha\beta}$ is the Ricci tensor. We take the Lorenz gauge condition

$\displaystyle\nabla_{\alpha}A^{\alpha}=0\ ,$ (2.2)

so (2.1) 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 [1] (page 570), [2] (page 93)

$\displaystyle
A^{\alpha}=\sum_{n=0}^{\infty}\text{Re}(A^{\alpha}_{n}\epsilon^{n}e^{iS/\epsilon})\
,$ (2.3)

where the amplitude is expanded in a series of the constant parameter
$\epsilon\ll 1$, 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 (2.3) into the equation of motion (2.1) gives, in the
Lorenz gauge,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!0$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\text{Re}\left[\epsilon^{n}\left(-\epsilon^{-2}k^{2}A^{\alpha}_{n}+i\epsilon^{-1}\nabla_{\beta}k^{\beta}A^{\alpha}_{n}+2i\epsilon^{-1}k^{\beta}\nabla_{\beta}A^{\alpha}_{n}+\Box
A^{\alpha}_{n}-R^{\alpha}{}_{\beta}A^{\beta}_{n}\right)e^{iS/\epsilon}\right]\
,$ (2.4)

where $k_{\alpha}\equiv\partial_{\alpha}S$ and $k^{2}\equiv
k_{\alpha}k^{\alpha}$. In the geometrical optics approximation, the
derivatives of $A^{\alpha}_{n}$ are taken to be of the same order of smallness
as $A^{\alpha}_{n}$, meaning that amplitude and polarisation change slowly
compared to the wavelength. We make the same assumption. In the geometrical
optics approximation, $R^{\alpha}{}_{\beta}$ 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^{\alpha}{}_{\beta}$ is of order
$\epsilon^{-2}$. We can still formally smoothly take the limit
$R^{\alpha}{}_{\beta}\to 0$ to recover geometrical optics. However, in
between, $R^{\alpha}{}_{\beta}$ would contribute at order $\epsilon^{-1}$, so
the approximation may not cover regimes of intermediate curvature properly.
Expanding the Lorenz gauge condition (2.2) in the same way, we have

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\text{Re}\left[\epsilon^{n}\left(\epsilon^{-1}ik_{\alpha}A^{\alpha}_{n}+\nabla_{\alpha}A^{\alpha}_{n}\right)e^{iS/\epsilon}\right]\
.$ (2.5)

Considering the equation of motion (2.4) at orders $\epsilon^{-2}$ and
$\epsilon^{-1}$ gives

$\displaystyle R_{\alpha\beta}a^{\beta}+k^{2}a_{\alpha}$ $\displaystyle=$
$\displaystyle 0$ (2.6) $\displaystyle
R_{\alpha\beta}b^{\beta}+k^{2}b_{\alpha}$ $\displaystyle=$ $\displaystyle
a_{\alpha}\nabla_{\beta}k^{\beta}+2k^{\beta}\nabla_{\beta}a_{\alpha}\ ,$ (2.7)

where we have denoted $A^{\alpha}_{0}\equiv a^{\alpha}$, $A^{\alpha}_{1}\equiv
ib^{\alpha}$. Expanding the Lorenz gauge condition (2.2) to orders
$\epsilon^{-1}$ and $\epsilon^{0}$ gives (for any vectors $v^{\alpha}$ and
$w^{\alpha}$, we denote $v\cdot w\equiv g_{\alpha\beta}v^{\alpha}w^{\beta}$)

$\displaystyle a\cdot k$ $\displaystyle=$ $\displaystyle 0$ (2.8)
$\displaystyle b\cdot k$ $\displaystyle=$
$\displaystyle\nabla_{\alpha}a^{\alpha}\ .$ (2.9)

This is the same result as in geometrical optics and shows that the light ray
is transverse. Contracting (2.7) with $a^{\alpha}$ and using (2.6), we obtain

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\nabla_{\beta}(a\cdot
a\,k^{\beta})\ ,$ (2.10)

so photon number is conserved, as in geometrical optics. In geometrical optics
$k^{\alpha}$ is null, so it follows from the transversality condition (2.8)
that $a^{\alpha}$ is spacelike or null. Here that is not the case, but we
assume, as in geometrical optics, that $a^{\alpha}$ is spacelike and decompose
it as $a^{\alpha}=af^{\alpha}$, where $a>0$ and $f\cdot f=1$. Using this
decomposition in (2.7) and applying (2.10), we find

$\displaystyle k^{\beta}\nabla_{\beta}f_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{1}{2a}(R_{\alpha\beta}b^{\beta}+k^{2}b_{\alpha})\ .$
(2.11)

This shows that the polarisation vector is not necessarily parallel
transported along the light ray. Using the decomposition of $a^{\alpha}$ in
(2.6) and contracting with $f^{\beta}$, we obtain

$\displaystyle k^{2}$ $\displaystyle=$ $\displaystyle-
R_{\alpha\beta}f^{\alpha}f^{\beta}\ ,$ (2.12)

so in general the light rays are not null. Taking a covariant derivative of
(2.12) and using the fact that
$\nabla_{\beta}k_{\alpha}=\nabla_{\alpha}k_{\beta}$ (since
$k_{\alpha}=\partial_{\alpha}S$) we get

$\displaystyle k^{\beta}\nabla_{\beta}k_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\nabla_{\alpha}k^{2}=-\frac{1}{2}\nabla_{\alpha}(R_{\alpha\beta}f^{\alpha}f^{\beta})\
,$ (2.13)

showing that in general the light rays are not geodesic.

##### The energy-momentum tensor.

Consider a unit timelike vector $u^{\alpha}$ ($u\cdot u=-1$). Without loss of
generality, we can decompose the energy-momentum tensor as

$\displaystyle
T_{\alpha\beta}=(\rho+p)u_{\alpha}u_{\beta}+pg_{\alpha\beta}+2q_{(\alpha}u_{\beta)}+\Pi_{\alpha\beta}\
,$ (2.14)

where $\rho$ is the energy density, $p$ is the pressure, $q_{\alpha}$ is the
energy flux, and $\Pi_{\alpha\beta}$ is the anisotropic stress. Both
$q_{\alpha}$ and $\Pi_{\alpha\beta}$ are orthogonal to $u^{\alpha}$, and
$\Pi_{\alpha\beta}$ is traceless and symmetric. From the Einstein equation
$R_{\alpha\beta}-\frac{1}{2}Rg_{\alpha\beta}=M_{{}_{\mathrm{Pl}}}^{-2}T_{\alpha\beta}$,
where $R\equiv g^{\alpha\beta}R_{\alpha\beta}$, we get

$\displaystyle M_{{}_{\mathrm{Pl}}}^{2}R_{\alpha\beta}$ $\displaystyle=$
$\displaystyle(\rho+p)u_{\alpha}u_{\beta}+\frac{1}{2}(\rho-p)g_{\alpha\beta}+2q_{(\alpha}u_{\beta)}+\Pi_{\alpha\beta}\
.$ (2.15)

Inserting (2.15) into (2.6) and using (2.12), we obtain

$\displaystyle[(\rho+p)u\cdot f+q\cdot f]u_{\alpha}$ $\displaystyle=$
$\displaystyle[(\rho+p)(u\cdot f)^{2}+2u\cdot fq\cdot
f+\Pi_{\gamma\delta}f^{\gamma}f^{\delta}]f_{\alpha}$ (2.16)
$\displaystyle-u\cdot fq_{\alpha}-\Pi_{\alpha\beta}f^{\beta}\ .$

As $u_{\alpha}$ 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_{\alpha}=0$. The right-hand
side condition then reduces to
$\Pi^{\alpha}{}_{\beta}f^{\beta}=\Pi_{\gamma\delta}f^{\gamma}f^{\delta}f^{\alpha}$.
By the same argument, the most reasonable solution is $\Pi_{\alpha\beta}=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^{\alpha}$ and the polarisation direction $f^{\alpha}$ are
relevant. (Although $u^{\alpha}$ is timelike, the projection of $f^{\alpha}$
orthogonal to it introduces a new spatial direction unless $u\cdot f=0$.) With
$u\cdot f=0$, $q_{\alpha}=0$, $\Pi_{\alpha\beta}=0$, (2.12) reduces to

$\displaystyle k^{2}$ $\displaystyle=$
$\displaystyle-\frac{\rho-p}{2M_{{}_{\mathrm{Pl}}}^{2}}\ .$ (2.17)

For matter satisfying $\rho\geq p$, the tangent vector $k^{\alpha}$ is
timelike or null, so the gravitational effect of matter slows down light,
analogously to refraction when passing through charged matter.

Applying (2.14) to the polarisation propagation equation (2.11) gives

$\displaystyle k^{\beta}\nabla_{\beta}f_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{1}{2M_{{}_{\mathrm{Pl}}}^{2}a}(\rho+p)u\cdot
b\,u_{\alpha}\ .$ (2.18)

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^{\alpha}$ 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 (2.13) for $k^{\alpha}$ reduces to the simple
expression

$\displaystyle k^{\beta}\nabla_{\beta}k_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{1}{4M_{{}_{\mathrm{Pl}}}^{2}}\partial_{\alpha}(\rho-p)\
.$ (2.19)

The results (2.17) and (2.19) allow us to calculate the post-geometrical
correction to redshift and angular diameter distance.

###  2.2 Redshift

####  2.2.1 Expression for the redshift

##### Dispersion relation.

We decompose the light ray tangent vector as

$\displaystyle k^{\alpha}$ $\displaystyle=$ $\displaystyle
E(u^{\alpha}+ve^{\alpha})\ ,$ (2.20)

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}{2M_{{}_{\mathrm{Pl}}}^{2}E^{2}}}$.
Correspondingly, the dispersion relation is

$\displaystyle E^{2}=\vec{k}^{2}+\frac{\rho-p}{2M_{{}_{\mathrm{Pl}}}^{2}}\ ,$
(2.21)

where $|\vec{k}|=vE$ 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\neq 0$, we can write
$k^{\alpha}=Mw^{\alpha}$, where $w^{\alpha}$ is a unit timelike vector, and
the propagation equation (2.13) reads
$w^{\beta}\nabla_{\beta}w_{\alpha}=-M^{-1}\hat{\nabla}_{\alpha}M$, where
$\hat{\nabla}_{\alpha}$ is the derivative projected orthogonally to
$w^{\alpha}$. This can be compared to the equation for an observer comoving
with an ideal fluid,
$w^{\beta}\nabla_{\beta}w_{\alpha}=-(\rho+p)^{-1}\hat{\nabla}_{\alpha}p$.

##### Redshift.

The redshift is given by

$\displaystyle 1+z=\frac{E_{\text{s}}}{E_{\text{o}}}\ ,$ (2.22)

where s refers to source and o to observation. To find the redshift, we
proceed as in geometrical optics (see e.g. [51, 52]), and take a derivative
with respect to the affine parameter $\lambda$ along the light ray,

$\displaystyle\frac{\mathrm{d}E}{\mathrm{d}\lambda}$ $\displaystyle=$
$\displaystyle-k^{\beta}\nabla_{\beta}\left(u^{\alpha}k_{\alpha}\right)=-k^{\alpha}k^{\beta}\nabla_{\beta}u_{\alpha}-\frac{1}{2}u^{\alpha}\nabla_{\alpha}k^{2}\
,$ (2.23)

where we have used (2.13). We decompose $\nabla_{\beta}u_{\alpha}$ as usual
(for overviews of the covariant formalism, see [53, 54])

$\displaystyle\nabla_{\beta}u_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\theta
h_{\alpha\beta}+\sigma_{\alpha\beta}+\omega_{\alpha\beta}-\dot{u}_{\alpha}u_{\beta}\
,$ (2.24)

where $h_{\alpha\beta}\equiv g_{\alpha\beta}+u_{\alpha}u_{\beta}$,
$\theta\equiv\nabla_{\alpha}u^{\alpha}$ is the volume expansion rate,
$\sigma_{\alpha\beta}$ is the traceless symmetric shear tensor,
$\omega_{\alpha\beta}$ is the antisymmetric vorticity tensor,
$\dot{u}^{\alpha}\equiv u^{\beta}\nabla_{\beta}u^{\alpha}$ is the acceleration
vector, and dot denotes $u^{\alpha}\nabla_{\alpha}$. The tensors
$\sigma_{\alpha\beta}$ and $\omega_{\alpha\beta}$ and the vector
$\dot{u}^{\alpha}$ are orthogonal to $u^{\alpha}$. Inserting (2.20) and (2.24)
into (2.23), we get

$\displaystyle\frac{1}{E^{2}}\frac{\mathrm{d}E}{\mathrm{d}\lambda}=-\frac{1}{3}v^{2}\theta-v^{2}\sigma_{\alpha\beta}e^{\alpha}e^{\beta}-ve\cdot\dot{u}-\frac{1}{2E^{2}}(k^{2})\dot{}\
.$ (2.25)

Integrating, we obtain

$\displaystyle\ln\left(1+z\right)$ $\displaystyle=$
$\displaystyle\int_{\lambda_{\text{s}}}^{\lambda_{\text{o}}}\mathrm{d}\lambda
E^{-1}\left[\frac{1}{3}v^{2}\theta+v^{2}\sigma_{\alpha\beta}e^{\alpha}e^{\beta}+ve\cdot\dot{u}+\frac{1}{2E^{2}}(k^{2})\dot{}\right]$
(2.26) $\displaystyle=$
$\displaystyle\int_{\lambda_{\text{s}}}^{\lambda_{\text{o}}}\mathrm{d}\lambda
E^{-1}\left[v^{2}\left(\frac{1}{3}\theta+\sigma_{\alpha\beta}e^{\alpha}e^{\beta}\right)-\frac{v}{\rho+p}e^{\alpha}\partial_{\alpha}p-\frac{\dot{\rho}-\dot{p}}{4M_{{}_{\mathrm{Pl}}}^{2}E^{2}}\right]\
,$

where in the second equality we have used (2.17) and the spatial part of the
continuity equation,

$\displaystyle\dot{u}_{\alpha}=-\frac{1}{\rho+p}h_{\alpha}{}^{\beta}\partial_{\beta}p\
.$ (2.27)

The result (2.26) 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.

####  2.2.2 Size of the curvature correction

##### Matter model.

Dark matter energy density consists of widely separated sharp spikes, which we
model as

$\displaystyle\rho=m\sum_{n}W(r_{n})\ ,$ (2.28)

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^{\alpha}$. We take $W$ to be a Gaussian of width equal to the Compton
wavelength $\lambda_{\text{c}}=2\pi/m$, so

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle
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}}}\
,$ (2.29)

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\neq 0$ for the particle not to collapse, and the pressure gradient may
be important. We assume $m\ll M_{{}_{\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.

##### 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 (2.26) 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.

###  2.3 Angular diameter distance

####  2.3.1 Expression for the angular diameter distance

##### 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^{\alpha}$
and to a four-velocity. We present the details in appendix A. 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^{\alpha}$
and $u^{\alpha}$,

$\displaystyle\tilde{h}_{\alpha\beta}$ $\displaystyle\equiv$ $\displaystyle
g_{\alpha\beta}+u_{\alpha}u_{\beta}-e_{\alpha}e_{\beta}$ (2.30)
$\displaystyle=$ $\displaystyle
g_{\alpha\beta}+\frac{k^{2}}{E^{2}+k^{2}}u_{\alpha}u_{\beta}+\frac{2E}{E^{2}+k^{2}}u_{(\alpha}k_{\beta)}-\frac{1}{E^{2}+k^{2}}k_{\alpha}k_{\beta}\
.$

We decompose $\nabla_{\beta}k_{\alpha}$ into parts projected in the direction
of and orthogonal to $\tilde{h}_{\alpha\beta}$ as

$\displaystyle\nabla_{\beta}k_{\alpha}=\frac{1}{2}\tilde{\theta}\tilde{h}_{\alpha\beta}+\tilde{\sigma}_{\alpha\beta}+P_{\alpha\beta}\
,$ (2.31)

where the area expansion rate is
$\tilde{\theta}\equiv\tilde{h}^{\alpha\beta}\nabla_{\beta}k_{\alpha}$, the
light shear is
$\tilde{\sigma}_{\alpha\beta}=\tilde{h}_{\alpha}{}^{\gamma}\tilde{h}_{\beta}{}^{\delta}\nabla_{\delta}k_{\gamma}-\frac{1}{2}\tilde{h}_{\alpha\beta}\tilde{\theta}$,
and
$P_{\alpha\beta}\equiv(\delta_{\alpha}{}^{\gamma}\delta_{\beta}{}^{\delta}-\tilde{h}_{\alpha}{}^{\gamma}\tilde{h}_{\beta}{}^{\delta})\nabla_{\delta}k_{\gamma}$.
Unlike in geometrical optics, $P_{\alpha\beta}$ is neither traceless nor
orthogonal to $k^{\alpha}$. As shown in appendix A we have
$D_{A}\propto\exp\left(\frac{1}{2}\int\mathrm{d}\lambda\nabla_{\alpha}k^{\alpha}\right)=\exp\left[\frac{1}{2}\int\mathrm{d}\lambda(\tilde{\theta}+P^{\alpha}{}_{\alpha})\right]$.

##### 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

$\displaystyle\frac{\mathrm{d}\tilde{\theta}}{\mathrm{d}\lambda}$
$\displaystyle=$ $\displaystyle
k^{\beta}\nabla_{\beta}\nabla_{\alpha}k^{\alpha}-\frac{\mathrm{d}P^{\alpha}{}_{\alpha}}{\mathrm{d}\lambda}$
(2.32) $\displaystyle=$
$\displaystyle-\frac{1}{2}\tilde{\theta}^{2}-2\tilde{\sigma}^{2}-R_{\alpha\beta}k^{\alpha}k^{\beta}+\frac{1}{2}\Box
k^{2}-P_{\alpha\beta}P^{\alpha\beta}-\frac{\mathrm{d}P^{\alpha}{}_{\alpha}}{\mathrm{d}\lambda}\
,$

where
$\tilde{\sigma}^{2}\equiv\frac{1}{2}\tilde{\sigma}_{\alpha\beta}\tilde{\sigma}^{\alpha\beta}$,
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_{\alpha\beta}P^{\alpha\beta}$ using
(2.30), we get

$\displaystyle P_{\alpha\beta}P^{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{v^{2}E^{2}}\left\\{\frac{k^{2}-E^{2}}{2v^{2}E^{2}}[(k^{2})\dot{}]^{2}+\frac{1}{v^{2}E}(k^{2})\dot{}(k^{2})^{\prime}-\frac{1}{4v^{2}E^{2}}[(k^{2})^{\prime}]^{2}-2k^{2}\dot{k}_{\alpha}\dot{k}^{\alpha}\right.$
(2.33)
$\displaystyle\left.-4Ek_{\alpha}^{\prime}\dot{k}^{\alpha}+2k_{\alpha}^{\prime}k^{\alpha}{}^{\prime}-\frac{\dot{k}\cdot
u}{v^{2}E}\left[2k^{2}(k^{2})\dot{}+E(k^{2})^{\prime}\right]-\frac{(k^{2})^{2}}{v^{2}E^{2}}(\dot{k}\cdot
u)^{2}\right\\}$ $\displaystyle=$
$\displaystyle\frac{1}{v^{2}E^{2}}\left\\{-\frac{1}{16M_{{}_{\mathrm{Pl}}}^{4}}[e^{\alpha}\partial_{\alpha}(\rho-p)]^{2}+\frac{\rho}{M_{{}_{\mathrm{Pl}}}^{2}}\dot{k}_{\alpha}\dot{k}^{\alpha}+\frac{E}{M_{{}_{\mathrm{Pl}}}^{2}}\dot{k}^{\alpha}\partial_{\alpha}(\rho-p)\right.$
$\displaystyle+\frac{1}{8M_{{}_{\mathrm{Pl}}}^{4}}\partial_{\alpha}(\rho-p)\partial^{\alpha}(\rho-p)+\frac{E^{2}}{2M_{{}_{\mathrm{Pl}}}^{2}\rho}e^{\alpha}\partial_{\alpha}p\,e^{\beta}\partial_{\beta}(\rho-p)$
$\displaystyle\left.-\frac{1}{4M_{{}_{\mathrm{Pl}}}^{4}}(e^{\alpha}\partial_{\alpha}p)^{2}\right\\}\
,$

where prime denotes $\frac{\mathrm{d}}{\mathrm{d}\lambda}$, and in the second
equality we have used (2.17), (2.19) and (2.27), 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_{\beta}k_{\alpha}$ with
$k^{\alpha}$ (entering from the definition (2.30) of
$\tilde{h}_{\alpha\beta}$), while both $u^{\alpha}$ and $k^{\alpha}$
contribute to the time derivatives. In geometrical optics, in a spatially flat
Friedmann–Lemaître–Robertson–Walker (FLRW) universe we have
$\nabla_{\alpha}k^{\alpha}=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^{\alpha}{}_{\alpha}$, and hence the difference between
$\nabla_{\alpha}k^{\alpha}$ and $\tilde{\theta}$ inside the integral. Using
(2.30), we have

$\displaystyle\int\mathrm{d}\lambda P^{\alpha}{}_{\alpha}$ $\displaystyle=$
$\displaystyle-\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})^{\prime}\right]\right\\}$
(2.34) $\displaystyle=$
$\displaystyle-\int\mathrm{d}\lambda\left\\{\frac{1}{2v^{2}M_{{}_{\mathrm{Pl}}}^{2}E^{2}}\left[\dot{u}\cdot
k\rho+\frac{1}{2}Ee^{\alpha}\partial_{\alpha}(\rho-p)\right]\right\\}$
$\displaystyle=$
$\displaystyle-\int\mathrm{d}\lambda\frac{1}{4vM_{{}_{\mathrm{Pl}}}^{2}E}e^{\alpha}\partial_{\alpha}(\rho-3p)\
,$

where on the second line we have used (2.17), 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 (2.27). 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^{\alpha}\partial_{\alpha}=\delta_{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^{\alpha}{}_{\alpha}$ is negligible, we have from
(2.32) and (2.33)

$\displaystyle\int\mathrm{d}\lambda\nabla_{\alpha}k^{\alpha}$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\lambda\tilde{\theta}$ (2.35) $\displaystyle=$
$\displaystyle\int\mathrm{d}\lambda\int\mathrm{d}\lambda^{\prime}\left(-\frac{1}{2}\tilde{\theta}^{2}-2\tilde{\sigma}^{2}-R_{\alpha\beta}k^{\alpha}k^{\beta}+\frac{1}{2}\Box
k^{2}-P_{\alpha\beta}P^{\alpha\beta}\right)$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\lambda\int\mathrm{d}\lambda^{\prime}\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}}$
$\displaystyle\underbrace{-\frac{1}{4M_{{}_{\mathrm{Pl}}}^{2}}\Box(\rho-p)}_{\text{term
2}}\underbrace{-P_{\alpha\beta}P^{\alpha\beta}}_{\text{term 3}}\Bigg{\\}}\ ,$

where we have used (2.15) and (2.17). The first two terms are the same as in
geometrical optics. Let us consider the three post-geometrical correction
terms in turn.

####  2.3.2 Size of the curvature correction

##### 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 (2.35) involves the integral

$\displaystyle\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)$
$\displaystyle=$ $\displaystyle\int_{\Delta
t}\mathrm{d}tE\frac{\rho}{M_{{}_{\mathrm{Pl}}}^{2}}\left(1-\frac{\rho}{4M_{{}_{\mathrm{Pl}}}^{2}E^{2}}\right)$
(2.36) $\displaystyle=$
$\displaystyle\frac{1}{A}\int_{V}\mathrm{d}^{3}rE\frac{\rho}{M_{{}_{\mathrm{Pl}}}^{2}}\left(1-\frac{\rho}{4M_{{}_{\mathrm{Pl}}}^{2}E^{2}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{A}\sum_{n}\int_{V}\mathrm{d}^{3}rE\frac{m}{M_{{}_{\mathrm{Pl}}}^{2}}W(r_{n})\left[1-\frac{m}{4M_{{}_{\mathrm{Pl}}}^{2}E^{2}}W(r_{n})\right]$
$\displaystyle=$ $\displaystyle\frac{\Delta
t}{V}EN\frac{m}{M_{{}_{\mathrm{Pl}}}^{2}}\left(1-\frac{\rho_{0}}{24\sqrt{\pi}M_{{}_{\mathrm{Pl}}}^{2}E^{2}}\right)$
$\displaystyle=$ $\displaystyle\int_{\Delta\lambda}\mathrm{d}\lambda
E^{2}\frac{\langle{\rho}\rangle}{M_{{}_{\mathrm{Pl}}}^{2}}\left(1-\frac{\rho_{0}}{24\sqrt{\pi}M_{{}_{\mathrm{Pl}}}^{2}E^{2}}\right)\
,$

where we have used $\mathrm{d}\lambda E=\mathrm{d}t$, taken into account
$p\ll\rho$, inserted the Gaussian density distribution (2.29), 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 $\langle{\rho}\rangle=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 [51, 52, 55]. 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}$.

##### Correction term 2.

The second post-geometrical term in (2.35) involves the integral

$\displaystyle\frac{1}{4M_{{}_{\mathrm{Pl}}}^{2}}\int_{\Delta\lambda}\mathrm{d}\lambda\Box(\rho-p)$
$\displaystyle=$
$\displaystyle\frac{1}{4M_{{}_{\mathrm{Pl}}}^{2}}\int\mathrm{d}\lambda\delta^{ij}\partial_{i}\partial_{j}(\rho-p)$
(2.37) $\displaystyle=$
$\displaystyle\frac{A^{-1}m}{4M_{{}_{\mathrm{Pl}}}^{2}}\sum_{n}\int_{V}\mathrm{d}^{3}rE^{-1}\delta^{ij}\partial_{i}\partial_{j}W(r_{n})$
$\displaystyle=$
$\displaystyle\frac{A^{-1}Nm}{4M_{{}_{\mathrm{Pl}}}^{2}}\sum_{n}\int_{V}\mathrm{d}^{3}rE^{-1}\delta^{ij}\partial_{i}\partial_{j}W(r)$
$\displaystyle=$
$\displaystyle\frac{A^{-1}Nm}{4M_{{}_{\mathrm{Pl}}}^{2}}\int_{V}\mathrm{d}^{3}rE^{-2}\partial_{r}W(r)\partial_{r}E\
,$

where we have taken into account $p\ll\rho$. As spacetime inside each particle
is static, $E\propto|g_{00}|^{-1/2}$. With $u^{\alpha}=|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 (2.27) 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$),

$\displaystyle\frac{\partial p}{\partial r}$ $\displaystyle=$ $\displaystyle-
G_{\mathrm{N}}\frac{(\rho+p)(\tilde{m}+4\pi
r^{3}p)}{r(r-2G_{\mathrm{N}}\tilde{m})}=-\frac{1}{8\pi
M_{{}_{\mathrm{Pl}}}^{2}}\frac{\rho\tilde{m}}{r^{2}}\ ,$ (2.38)

where $\tilde{m}(r)\equiv\int_{0}^{r}\mathrm{d}r^{\prime}4\pi r^{\prime
2}\rho$, and in the second equality we have taken into account that for a
Gaussian density distribution and $m\ll M_{{}_{\mathrm{Pl}}}$, the Newtonian
limit is a good approximation. Therefore

$\displaystyle\frac{\partial E}{\partial r}$ $\displaystyle=$
$\displaystyle-\frac{E}{8\pi M_{{}_{\mathrm{Pl}}}^{2}}\frac{\tilde{m}}{r^{2}}\
,$ (2.39)

Applying (2.39), the integral (2.37) becomes

$\displaystyle\int_{\Delta\lambda}\mathrm{d}\lambda\frac{1}{4M_{{}_{\mathrm{Pl}}}^{2}}\Box(\rho-p)$
$\displaystyle=$
$\displaystyle-\frac{A^{-1}NmE^{-1}}{8M_{{}_{\mathrm{Pl}}}^{4}}\int_{0}^{\infty}\mathrm{d}r\tilde{m}(r)\partial_{r}W(r)$
(2.40) $\displaystyle=$
$\displaystyle\frac{A^{-1}NmE^{-1}}{8M_{{}_{\mathrm{Pl}}}^{4}}\int_{0}^{\infty}\mathrm{d}r\partial_{r}\tilde{m}(r)W(r)$
$\displaystyle=$
$\displaystyle\int_{\Delta\lambda}\mathrm{d}\lambda\frac{\rho^{2}}{8M_{{}_{\mathrm{Pl}}}^{2}}\
,$

where we have neglected the subleading variation of $E$. This term cancels
half of the correction term in (2.36).

##### Correction term 3.

The third post-geometrical term in (2.35) involves an integral over (2.33). It
contains $\dot{k}^{\alpha}$, which we can write as

$\displaystyle\dot{k}_{\alpha}$ $\displaystyle=$ $\displaystyle
u^{\beta}\nabla_{\beta}k_{\alpha}$ (2.41) $\displaystyle=$
$\displaystyle-\partial_{\alpha}E-Ev\left(\frac{1}{3}\theta
e_{\alpha}+\sigma_{\alpha\beta}e^{\beta}+\omega_{\alpha\beta}e^{\beta}-\dot{u}\cdot
e\,u_{\alpha}\right)$ $\displaystyle=$
$\displaystyle-\partial_{\alpha}E-Ev\rho^{-1}e^{\beta}\partial_{\beta}p\,u_{\alpha}\
,$

where we have used $\nabla_{\beta}k_{\alpha}=\nabla_{\alpha}k_{\beta}$,
$E=-u\cdot k$, the decomposition (2.24), and on the third line we have used
(2.27), and taken into account $p\ll\rho$ and that inside the particles
$\theta=0$ and $\sigma_{\alpha\beta}=\omega_{\alpha\beta}=0$. Using (2.41) in
(2.33), the integral over $P_{\alpha\beta}P^{\alpha\beta}$ reads

$\displaystyle\int\mathrm{d}\lambda\frac{1}{M_{{}_{\mathrm{Pl}}}^{4}E^{2}v^{2}}\left[-\frac{1}{16}(e^{\alpha}\partial_{\alpha}\rho)^{2}+\frac{1}{8}\left(1+\frac{4M_{{}_{\mathrm{Pl}}}^{2}E^{2}}{\rho}\right)e^{\alpha}\partial_{\alpha}\rho\,e^{\beta}\partial_{\beta}p\right.$
$\displaystyle-\frac{1}{16}\left(5+24\frac{M_{{}_{\mathrm{Pl}}}^{2}E^{2}}{\rho}\right)(e^{\alpha}\partial_{\alpha}p)^{2}+\frac{1}{8}\partial_{\alpha}(\rho-p)\partial^{\alpha}(\rho-p)$
$\displaystyle\left.-M_{{}_{\mathrm{Pl}}}^{2}E\partial_{\alpha}E\partial^{\alpha}(\rho-p)+M_{{}_{\mathrm{Pl}}}^{2}\rho\partial_{\alpha}E\partial^{\alpha}E\right]\
.$ (2.42)

The directional derivatives in the first three terms again become sums over
dot products. For each particle $n$, we have
$e^{\alpha}\partial_{\alpha}=\delta_{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}\delta_{ij}e^{i}(r^{j}-r_{m}^{j})/r_{m}\delta_{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

$\displaystyle\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{4M_{{}_{\mathrm{Pl}}}^{2}E^{2}}{\rho}\right)\partial_{r_{n}}\rho\,\partial_{r_{n}}p\right.$
$\displaystyle\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]\
.$ (2.43)

Using $\partial_{r_{n}}E=E\rho^{-1}\partial_{r_{n}}p$, (2.3.2) becomes

$\displaystyle\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{20M_{{}_{\mathrm{Pl}}}^{2}E^{2}}{\rho}\right)\partial_{r_{n}}\rho\,\partial_{r_{n}}p\right.$
$\displaystyle\left.+\frac{1}{48}\left(1+\frac{72M_{{}_{\mathrm{Pl}}}^{2}E^{2}}{\rho}\right)(\partial_{r_{n}}p)^{2}\right]\
,$ (2.44)

where the pressure gradient is given in (2.38). 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)^{2}r\partial_{r}p$ in the
integral (2.3.2) 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
(2.3.2) reduce to volume integrals involving factors of $\rho$,
$\partial_{r}\rho$, $\tilde{m}$ and $r$. Using the Gaussian density
distribution (2.29) and the pressure gradient (2.38), we have

$\displaystyle\int\mathrm{d}^{3}r(\partial_{r}\rho)^{2}$ $\displaystyle=$
$\displaystyle\frac{m}{16\pi^{5/2}}m^{2}\rho_{0}$
$\displaystyle\int\mathrm{d}^{3}r\partial_{r}\rho\partial_{r}p$
$\displaystyle=$
$\displaystyle\frac{m}{54\sqrt{3}\pi}\frac{\rho_{0}^{2}}{M_{{}_{\mathrm{Pl}}}^{2}}$
$\displaystyle\int\mathrm{d}^{3}r\rho^{-1}\partial_{r}\rho\partial_{r}p$
$\displaystyle=$
$\displaystyle\frac{m}{12\sqrt{\pi}}\frac{\rho_{0}}{M_{{}_{\mathrm{Pl}}}^{2}}$
$\displaystyle\int\mathrm{d}^{3}r\rho^{-1}(\partial_{r}\rho)^{2}r\partial_{r}p$
$\displaystyle=$
$\displaystyle-\frac{m}{27\sqrt{3}\pi}\frac{\rho_{0}^{2}}{M_{{}_{\mathrm{Pl}}}^{2}}$
$\displaystyle\int\mathrm{d}^{3}r(\partial_{r}p)^{2}$ $\displaystyle\approx$
$\displaystyle 3\times
10^{-9}m\frac{m^{4}}{M_{{}_{\mathrm{Pl}}}^{4}}m^{2}\rho_{0}$
$\displaystyle\int\mathrm{d}^{3}r\rho^{-1}(\partial_{r}p)^{2}$
$\displaystyle\approx$ $\displaystyle 4\times
10^{-5}m\frac{m^{2}}{M_{{}_{\mathrm{Pl}}}^{2}}\frac{\rho_{0}}{M_{{}_{\mathrm{Pl}}}^{2}}\
.$ (2.45)

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}\sim
10^{-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 $\langle{\rho}\rangle$ as before. Overall, (2.3.2) becomes

$\displaystyle\int\mathrm{d}\lambda
E^{2}\frac{\langle{\rho}\rangle}{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]\
.$ (2.46)

##### The total correction.

In summary, gathering the terms (2.36), (2.40) and (2.46), the integral (2.35)
is

$\displaystyle\int\mathrm{d}\lambda\tilde{\theta}$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\lambda\int\mathrm{d}\lambda^{\prime}\left\\{-\frac{1}{2}\tilde{\theta}^{2}-2\tilde{\sigma}^{2}-E^{2}\frac{\langle{\rho}\rangle}{M_{{}_{\mathrm{Pl}}}^{2}}\left[1+\frac{5}{8192\pi^{13/2}}\frac{m^{6}}{M_{{}_{\mathrm{Pl}}}^{2}E^{4}}\right.\right.$
(2.47)
$\displaystyle\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\\}\
.$

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 (2.47) 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 (2.47) 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\gtrsim 10(M_{{}_{\mathrm{Pl}}}E^{2})^{1/3}$. For the longest
radio waves observable through the Earth’s ionosphere, this gives $m\gtrsim
100$ keV, and for CMB photons we get $m\gtrsim 100$ 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.

##### From $\tilde{\theta}$ to $D_{A}$.

In terms of the angular diameter distance $D_{A}$, (2.47) corresponds to (see
e.g. [51, 52]), dropping those correction terms that are important only near
the curvature limit, and neglecting the null shear (which is small in the real
universe),

$\displaystyle D_{A}^{\prime\prime}$ $\displaystyle=$ $\displaystyle-4\pi
G_{\mathrm{N}}\langle{\rho}\rangle\left[1+\alpha(E/E_{\text{o}})^{-4}\right]D_{A}\
,$ (2.48)

where
$\alpha\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

$\displaystyle
H(z)\frac{\mathrm{d}}{\mathrm{d}z}\left[(1+z)^{2}H(z)\frac{\mathrm{d}D_{A}}{\mathrm{d}z}\right]$
$\displaystyle=$ $\displaystyle-4\pi
G_{\mathrm{N}}\langle{\rho}\rangle\left[1+\alpha(1+z)^{-4}\right]D_{A}\ .$
(2.49)

Because $\langle{\rho}\rangle\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 $\alpha>0$ in
(2.49) 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 [51, 52, 56, 55]. A value $\alpha\gtrsim 1$ 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 $\alpha>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) [57], because the inferred value of $H_{0}^{-1}$ from the CMB is larger.
(If $\alpha$ 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 $\langle{\rho}\rangle$ usually also
includes baryonic matter. (As the full term is
$\langle{\rho}\rangle+\langle{p}\rangle$, 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)}$ [58]. 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^{\prime}}{\sqrt{\Omega_{\Lambda}+(\Omega_{\text{dm}}+\Omega_{\text{b}})(1+z^{\prime})^{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\pm 1.3$
km/s/Mpc [59] as the real value, the inferred value is $69.0\pm 1.2$ km/s/Mpc,
closer to the value $67.4\pm 0.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 B 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.

###  2.4 Discussion

##### 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 (2.10), which gives
$a^{\prime}=-\frac{1}{2}a\nabla_{\alpha}k^{\alpha}$. (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 (2.18) with $k^{\alpha}$ and applying the decomposition (2.24)
gives the estimate
$f_{\alpha}^{\prime}\sim(Em^{3})/M_{{}_{\mathrm{Pl}}}^{2}\sim(EM^{2})/m$. The
second derivative of the polarisation in the fourth term in (2.4) is $\sim
E^{-2}f_{\alpha}^{\prime\prime}\sim E^{-1}mf_{\alpha}^{\prime}\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/(mv_{\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}\sim 10^{2}$ (taking $v_{\text{dm}}\sim 10^{-3}$), to
$m\sim 10^{6}$ GeV. The limiting mass for the leading correction to the
angular diameter distance would increase by the factor
$v_{\text{dm}}^{-5/6}\sim 10^{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 [60, 61, 62, 63, 64], 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.

##  3 Conclusions

##### 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}}\sim 1/m$ inside the particle. The origin of the
derivatives is the geodesic deviation vector
$k^{\beta}\nabla_{\beta}k^{\alpha}$, 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\gtrsim
10(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\gtrsim 10\sqrt{M_{{}_{\mathrm{Pl}}}E}\sim 10^{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 [47, 48, 49]. 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\gtrsim 100$ 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.

##  Appendix A Calculating $D_{A}$

In this appendix we relate the beam area expansion rate to the angular
diameter distance $D_{A}$, generalising the geometrical optics calculation
(see e.g. [2] pages 104-115). Let $\delta x^{\alpha}$ be a separation vector
connecting two rays in an infinitesimal light beam. Note that in order for the
difference in the position of two rays to be a vector, the rays have to belong
to the same tangent space, i.e. the curvature has be small. When $E\sim M$,
the curvature is large, and the validity of the treatment is questionable,
whereas for $E\gg M$ the curvature is small. As seen in section 2.3.2, in
order to reproduce the average cosmological density of matter in agreement
with observations we have to consider beams that are larger than the
separation of dark matter particles. We simply assume that we can sum up the
effects of nearby light rays, but it would be interesting to consider
rigorously how the infinitesimal beams can be added, and how the tangent
spaces orthogonal to the light rays mesh together.

By definition, $\delta x^{\alpha}$ is dragged along the light rays, so its Lie
derivative with respect to the tangent vector $k^{\alpha}$ vanishes,

$\displaystyle k^{\beta}\nabla_{\beta}\delta x^{\alpha}-\delta
x^{\beta}\nabla_{\beta}k^{\alpha}=0\ .$ (A.1)

Let $\tilde{U}^{\alpha}$ be a timelike vector dragged along the light rays, so
that its Lie derivative with respect to $k^{\alpha}$ also vanishes,

$\displaystyle
k^{\beta}\nabla_{\beta}\tilde{U}^{\alpha}-\tilde{U}^{\beta}\nabla_{\beta}k^{\alpha}=0\
.$ (A.2)

The unit normalised four-velocity parallel to $\tilde{U}^{\alpha}$ is
$U^{\alpha}=N\tilde{U}^{\alpha}$ (here $N\equiv
1/\sqrt{-\tilde{U}\cdot\tilde{U}}$, so that $U\cdot U=-1$):

$\displaystyle k^{\alpha}$ $\displaystyle=$ $\displaystyle
E_{(U)}(U^{\alpha}+v_{(U)}e_{(U)}^{\alpha})\ ,$ (A.3)

where $E_{(U)}=-U\cdot k$, $U\cdot e_{(U)}=0$, $e_{(U)}\cdot e_{(U)}=1$, and
$v_{(U)}\equiv\sqrt{1+k^{2}/E_{(U)}^{2}}$. Along the light ray, $\delta
x^{\alpha}$ will not stay orthogonal to $k^{\alpha}$ even if it was initially
so. We therefore introduce the projection tensor orthogonal both to
$k^{\alpha}$ and $U^{\alpha}$, analogous to (2.30),

$\displaystyle\tilde{h}_{\alpha\beta}^{(U)}$ $\displaystyle\equiv$
$\displaystyle g_{\alpha\beta}+U_{\alpha}U_{\beta}-e_{(U)\alpha}e_{(U)\beta}$
(A.4) $\displaystyle=$ $\displaystyle
g_{\alpha\beta}+\frac{k^{2}}{E_{(U)}^{2}+k^{2}}U_{\alpha}U_{\beta}+\frac{2E_{(U)}}{E_{(U)}^{2}+k^{2}}U_{(\alpha}k_{\beta)}-\frac{1}{E_{(U)}^{2}+k^{2}}k_{\alpha}k_{\beta}\
.$

The orthogonally projected separation vector is
$h_{\alpha}\equiv\tilde{h}_{\alpha\beta}\delta x^{\beta}\equiv
h\hat{h}_{\alpha}$, where $\hat{h}^{\alpha}\hat{h}_{\alpha}=1$. The norm $h>0$
gives the distance between light rays in the bundle. Let us consider its
evolution along the light ray,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\lambda}h^{2}$ $\displaystyle=$
$\displaystyle
2h^{\alpha}k^{\gamma}\nabla_{\gamma}(\tilde{h}_{\alpha\beta}^{(U)}\delta
x^{\beta})$ (A.5) $\displaystyle=$ $\displaystyle
2h_{\beta}k^{\gamma}\nabla_{\gamma}\delta x^{\beta}+2h^{\alpha}\delta
x^{\beta}k^{\gamma}\nabla_{\gamma}\tilde{h}_{\alpha\beta}$ $\displaystyle=$
$\displaystyle 2h_{\beta}\delta
x^{\gamma}\nabla_{\gamma}k^{\beta}+2h^{\alpha}\delta
x^{\beta}k^{\gamma}\nabla_{\gamma}\tilde{h}_{\alpha\beta}$ $\displaystyle=$
$\displaystyle 2h^{\alpha}h^{\beta}\nabla_{\beta}k_{\alpha}+2\frac{k^{2}\delta
x\cdot u+E_{(U)}\delta x\cdot
k}{E_{(U)}^{2}+k^{2}}h^{\alpha}(k^{\beta}\nabla_{\beta}U_{\alpha}-U^{\beta}\nabla_{\beta}k_{\alpha})$
$\displaystyle=$ $\displaystyle
2h^{\alpha}h^{\beta}\nabla_{\beta}k_{\alpha}+2\frac{k^{2}\delta x\cdot
u+E_{(U)}\delta x\cdot
k}{E_{(U)}^{2}+k^{2}}h^{\alpha}\tilde{U}_{\alpha}k^{\beta}\nabla_{\beta}N$
$\displaystyle=$ $\displaystyle 2h^{\alpha}h^{\beta}\nabla_{\beta}k_{\alpha}\
,$

where we have on the third line used (A.1), on the fourth line written $\delta
x^{\gamma}=h^{\gamma}+(\delta^{\gamma}{}_{\alpha}-\tilde{h}^{\gamma}{}_{\alpha})\delta
x^{\alpha}$ and used the definition (A.4), on the next to last line written
$U^{\alpha}=N\tilde{U}^{\alpha}$ and used the vanishing of the Lie derivative
(A.2), and on the last line taken into account that $h^{\alpha}$ and
$U^{\alpha}$ are orthogonal.

Analogously to (2.31), we decompose $\nabla_{\beta}k_{\alpha}$ into parts
projected in the direction of and orthogonal to
$\tilde{h}_{\alpha\beta}^{(U)}$ as

$\displaystyle\nabla_{\beta}k_{\alpha}=\frac{1}{2}\tilde{\theta}^{(U)}\tilde{h}_{\alpha\beta}^{(U)}+\tilde{\sigma}_{\alpha\beta}^{(U)}+P_{\alpha\beta}^{(U)}\
,$ (A.6)

where the area expansion rate is
$\tilde{\theta}^{(U)}\equiv\tilde{h}^{(U)\alpha\beta}\nabla_{\beta}k_{\alpha}$,
the light shear is
$\tilde{\sigma}_{\alpha\beta}^{(U)}=\tilde{h}^{(U)}_{\alpha\gamma}\tilde{h}_{\beta\delta}\nabla^{\delta}k^{\gamma}-\frac{1}{2}\tilde{h}_{\alpha\beta}^{(U)}\tilde{\theta}^{(U)}$,
and
$P_{\alpha\beta}^{(U)}\equiv(g_{\alpha\gamma}g_{\beta\delta}-\tilde{h}_{\alpha\gamma}^{(U)}\tilde{h}_{\beta\delta}^{(U)})\nabla^{\delta}k^{\gamma}$.
In geometrical optics, $P_{\alpha\beta}^{(U)}$ is orthogonal to $k^{\alpha}$
and can be written as $P_{\alpha\beta}^{(U)}=2k_{(\alpha}P_{\beta)}^{(U)}$,
with $P^{(U)}\cdot k=0$. That is not the case here.

Inserting (A.6) into (A.5), we obtain

$\displaystyle\frac{1}{h}\frac{\mathrm{d}h}{\mathrm{d}\lambda}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\tilde{\theta}^{(U)}+\tilde{\sigma}_{\alpha\beta}^{(U)}\hat{h}^{\alpha}\hat{h}^{\beta}\
,$ (A.7)

from which it follows that
$\tilde{\theta}=\frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}\lambda}$ describes
the relative rate of change of the beam cross section area $A$. Since the
angular diameter distance is defined as $D_{A}^{2}=A/\Omega_{\text{o}}$, where
$\Omega_{\text{o}}$ is the solid angle at observation, we obtain the usual
relation

$\displaystyle
D_{A}\propto\exp\left(\frac{1}{2}\int\mathrm{d}\lambda\tilde{\theta}^{(U)}\right)\
.$ (A.8)

Taking the trace of (A.6), we can split the area expansion rate as

$\displaystyle\tilde{\theta}^{(U)}$ $\displaystyle=$
$\displaystyle\nabla_{\alpha}k^{\alpha}-P^{(U)\alpha}{}_{\alpha}\ .$ (A.9)

Velocity-dependence of the area expansion rate is contained in the second term
(which is zero in geometrical optics). Using (A.4), we can write it as

$\displaystyle-P^{(U)\alpha}{}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{E_{(U)}^{2}+k^{2}}\left(k^{2}U^{\alpha}U^{\beta}\nabla_{\beta}k_{\alpha}+E_{(U)}U^{\alpha}\nabla_{\alpha}k^{2}-\frac{1}{2}k^{\alpha}\nabla_{\alpha}k^{2}\right)\
.$ (A.10)

We can rewrite the first term using the identity (which it is straightforward
to verify)

$\displaystyle E_{(U)}U^{\alpha}U^{\beta}\nabla_{\beta}k_{\alpha}$
$\displaystyle=$ $\displaystyle
k^{\alpha}\nabla_{\alpha}E_{(U)}+U^{\alpha}\nabla_{\alpha}k^{2}\ .$ (A.11)

Inserting (A.11) into (A.10), we obtain

$\displaystyle-P^{(U)\alpha}{}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{E_{(U)}(E_{(U)}^{2}+k^{2})}k^{\alpha}\nabla_{\alpha}E_{(U)}+\frac{1}{E_{(U)}}U^{\alpha}\nabla_{\alpha}k^{2}-\frac{1}{2}\frac{1}{E_{(U)}^{2}+k^{2}}k^{\alpha}\nabla_{\alpha}k^{2}$
(A.12) $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}\lambda}\left(\ln\frac{E_{(U)}^{2}}{E_{(U)}^{2}+k^{2}}\right)+\frac{1}{E_{(U)}}U^{\alpha}\nabla_{\alpha}k^{2}\
.$

When integrating from the source to the observer, the first term gives zero if
we assume that they are no coincident with a dark matter particle, so that
$k^{2}=0$. We can decompose $U^{\alpha}$ as

$\displaystyle U^{\alpha}=\gamma(u^{\alpha}+s^{\alpha})\ ,$ (A.13)

where $\gamma\equiv-u\cdot U$, $u\cdot s=0$. We thus have

$\displaystyle-\int\mathrm{d}\lambda P^{(U)\alpha}{}_{\alpha}$
$\displaystyle=$
$\displaystyle-\int\mathrm{d}\lambda\frac{1}{2M_{{}_{\mathrm{Pl}}}^{2}E_{(U)}}\gamma
s^{\alpha}\partial_{\alpha}(\rho-p)\ ,$ (A.14)

where we have used (2.17) and the fact that the density and pressure are
independent of the time defined by $u^{\alpha}$ (but not of the time defined
by $U^{\alpha}$). As this term is linear in the density, it is parametrically
suppressed compared to the other terms we have found that are quadratic in
density. However, it is also suppressed by a sum over directions. As in
section 2.3.1, the directional derivative becomes a sum over dot products, in
this case $\delta_{ij}(r^{i}-r^{i}_{n})s^{j}/r_{n}$. As the turning of
$U^{\alpha}$ (and hence $s^{i}$) along the light ray has no preferred spatial
direction with respect to $u^{\alpha}$ apart from $e^{\alpha}$, and the term
is linear in the cosine of the angle, it vanishes.

Therefore the velocity-dependent term $P^{(U)\alpha}{}_{\alpha}$ can be
neglected when calculating $D_{A}$, and it is sufficient to consider the term
$\nabla_{\alpha}k^{\alpha}$, which is independent of velocity.

##  Appendix B Distance duality

In this appendix, we show that the reciprocity relation and hence distance
duality $D_{L}=(1+z)^{2}D_{A}$, where $D_{L}$ is the luminosity distance, are
unaffected by the correction terms in our post-geometrical optics
approximation. We follow the proof for the reciprocity relation in the
geometrical optics approximation outlined in e.g. [2] (pages 110-115), [65]
(pages 165-166).

Consider two beams between source and observer, one that converges at the
observer and another that converges at the source. We denote a deviation
vector in the beam that converges at the source by 1 and a deviation vector in
the beam that converges in the observer by 2. The two beams have common light
rays. Consider the following quantity along such a ray:

$\displaystyle\delta x^{\alpha}_{1}k^{\beta}\nabla_{\beta}\delta
x_{2\alpha}-\delta x^{\alpha}_{2}k^{\beta}\nabla_{\beta}\delta x_{1\alpha}$
$\displaystyle=$ $\displaystyle(\delta x^{\alpha}_{2}\delta
x^{\beta}_{1}-\delta x^{\alpha}_{1}\delta
x^{\beta}_{2})\nabla_{\beta}k_{\alpha}=0\ ,$ (B.1)

where in the first equality we have used the vanishing of the Lie derivative
(A.1), and the second equality follows from the fact that
$\nabla_{\beta}k_{\alpha}$ is symmetric. Given the equality $\delta
x^{\alpha}_{1}k^{\beta}\nabla_{\beta}\delta x_{2\alpha}=\delta
x^{\alpha}_{2}k^{\beta}\nabla_{\beta}\delta x_{1\alpha}$, the proof of the
reciprocity relation between the area and subtended angle at observer and
source, and hence distance duality, follows in the same way as in geometrical
optics. Unlike in geometrical optics, the screen area is not independent of
observer velocity where $k^{\alpha}$ is not null. However, this would only
make a difference if a significant density of matter overlapped the observer,
which is not the case for particle dark matter (at least of the non-fuzzy
variety where the density at particle location is large enough to make a
difference).

## References

  * [1] C. Misner, K. Thorne and J. Wheeler, _Gravitation_.  Princeton University Press, 1973.
  * [2] P. Schneider, J. Ehlers and E. Falco, _Gravitational Lenses_.  Astronomy and Astrophysics Library. Springer New York, 1992.
  * [3] T. T. Nakamura and S. Deguchi, _Wave Optics in Gravitational Lensing_ , _Progress of Theoretical Physics Supplement_ 133 (01, 1999) 137–153.
  * [4] Y. Nambu, _Wave Optics and Image Formation in Gravitational Lensing_ , _J. Phys. Conf. Ser._ 410 (2013) 012036, [1207.6846].
  * [5] S. G. Turyshev and V. T. Toth, _Diffraction of electromagnetic waves in the gravitational field of the Sun_ , _Phys. Rev. D_ 96 (2017) 024008, [1704.06824].
  * [6] S. G. Turyshev and V. T. Toth, _Diffraction of electromagnetic waves by an extended gravitational lens_ , _Phys. Rev. D_ 103 (2021) 064076, [2102.03891].
  * [7] S. G. Turyshev and V. T. Toth, _Gravitational lensing by an extended mass distribution_ , _Phys. Rev. D_ 104 (2021) 044013, [2106.06696].
  * [8] T. T. Nakamura, _Gravitational lensing of gravitational waves from inspiraling binaries by a point mass lens_ , _Phys. Rev. Lett._ 80 (Feb, 1998) 1138–1141.
  * [9] R. Takahashi and T. Nakamura, _Wave effects in gravitational lensing of gravitational waves from chirping binaries_ , _Astrophys. J._ 595 (2003) 1039–1051, [astro-ph/0305055].
  * [10] J.-P. Macquart, _Scattering of gravitational radiation: Second order moments of the wave amplitude_ , _Astron. Astrophys._ 422 (2004) 761–775, [astro-ph/0402661].
  * [11] R. Takahashi, _Amplitude and phase fluctuations for gravitational waves propagating through inhomogeneous mass distribution in the universe_ , _Astrophys. J._ 644 (2006) 80–85, [astro-ph/0511517].
  * [12] R. Takahashi, _Arrival time differences between gravitational waves and electromagnetic signals due to gravitational lensing_ , _Astrophys. J._ 835 (2017) 103, [1606.00458].
  * [13] G. Cusin and M. Lagos, _Gravitational wave propagation beyond geometric optics_ , _Phys. Rev. D_ 101 (2020) 044041, [1910.13326].
  * [14] B. S. DeWitt and R. W. Brehme, _Radiation damping in a gravitational field_ , _Annals Phys._ 9 (1960) 220–259.
  * [15] F. Cabral and F. S. N. Lobo, _Electrodynamics and Spacetime Geometry: Foundations_ , _Found. Phys._ 47 (2017) 208–228, [1602.01492].
  * [16] F. A. Asenjo and S. A. Hojman, _Do electromagnetic waves always propagate along null geodesics?_ , _Class. Quant. Grav._ 34 (2017) 205011, [1608.06572].
  * [17] C. J. Copi, K. Pasmatsiou and G. D. Starkman, _Scalar and vector tail radiation from the interior of the lightcone_ , _JCAP_ 01 (2021) 050, [2008.13069].
  * [18] P. D. Mannheim, _Critique of the use of geodesics in astrophysics and cosmology_ , 2105.08556.
  * [19] A. M. Anile, _Geometrical optics in general relativity: A study of the higher order corrections_ , _Journal of Mathematical Physics_ 17 (1976) 576–584.
  * [20] S. R. Dolan, _Higher-order geometrical optics for electromagnetic waves on a curved spacetime_ , 1801.02273.
  * [21] M. A. Oancea, C. F. Paganini, J. Joudioux and L. Andersson, _An overview of the gravitational spin Hall effect_ , 1904.09963.
  * [22] M. A. Oancea, J. Joudioux, I. Y. Dodin, D. E. Ruiz, C. F. Paganini and L. Andersson, _Gravitational spin Hall effect of light_ , _Phys. Rev. D_ 102 (2020) 024075, [2003.04553].
  * [23] V. P. Frolov, _Maxwell equations in a curved spacetime: Spin optics approximation_ , _Phys. Rev. D_ 102 (2020) 084013, [2007.03743].
  * [24] Y. B. Zel’dovich, _Observations in a Universe Homogeneous in the Mean_ , _Soviet Astronomy_ 8 (Aug., 1964) 13.
  * [25] V. M. Dashevskii and Y. B. Zel’dovich, _Propagation of Light in a Nonhomogeneous Nonflat Universe II_ , _Soviet Astronomy_ 8 (June, 1965) 854.
  * [26] V. M. Dashevskii and V. I. Slysh, _On the Propagation of Light in a Nonhomogeneous Universe_ , _Soviet Astronomy_ 9 (Feb., 1966) 671.
  * [27] B. Bertotti, _The Luminosity of Distant Galaxies_ , _Proceedings of the Royal Society of London Series A_ 294 (Sept., 1966) 195–207.
  * [28] J. E. Gunn, _On the Propagation of Light in Inhomogeneous Cosmologies. I. Mean Effects_ , _Astrophys. J._ 150 (Dec., 1967) 737.
  * [29] C. C. Dyer and R. C. Roeder, _Observations in Locally Inhomogeneous Cosmological Models_ , _Astrophys. J._ 189 (Apr., 1974) 167–176.
  * [30] C. C. Dyer and R. C. Roeder, _On the transition from Weyl to Ricci focusing_ , _General Relativity and Gravitation_ 13 (Dec., 1981) 1157–1160.
  * [31] M. Kasai, T. Futamase and F. Takahara, _Angular diameter distance in a clumpy universe_ , _Phys. Lett. A_ 147 (1990) 97–105.
  * [32] T. W. B. Kibble and R. Lieu, _Average magnification effect of clumping of matter_ , _Astrophys. J._ 632 (2005) 718–726, [astro-ph/0412275].
  * [33] T. Clifton and P. G. Ferreira, _Archipelagian Cosmology: Dynamics and Observables in a Universe with Discretized Matter Content_ , _Phys. Rev. D_ 80 (2009) 103503, [0907.4109].
  * [34] T. Clifton and P. G. Ferreira, _Errors in Estimating $\Omega_{\Lambda}$ due to the Fluid Approximation_, _JCAP_ 10 (2009) 026, [0908.4488].
  * [35] T. Clifton, _Cosmology Without Averaging_ , _Class. Quant. Grav._ 28 (2011) 164011, [1005.0788].
  * [36] C. Clarkson, G. F. R. Ellis, A. Faltenbacher, R. Maartens, O. Umeh and J.-P. Uzan, _(Mis-)Interpreting supernovae observations in a lumpy universe_ , _Mon. Not. Roy. Astron. Soc._ 426 (2012) 1121–1136, [1109.2484].
  * [37] T. Clifton, P. G. Ferreira and K. O’Donnell, _An Improved Treatment of Optics in the Lindquist-Wheeler Models_ , _Phys. Rev. D_ 85 (2012) 023502, [1110.3191].
  * [38] J.-P. Bruneton and J. Larena, _Dynamics of a lattice Universe: The dust approximation in cosmology_ , _Class. Quant. Grav._ 29 (2012) 155001, [1204.3433].
  * [39] J.-P. Bruneton and J. Larena, _Observables in a lattice Universe_ , _Class. Quant. Grav._ 30 (2013) 025002, [1208.1411].
  * [40] J. Larena, _The fitting problem in a lattice Universe_ , _Springer Proc. Phys._ 157 (2014) 385–392, [1210.2161].
  * [41] R. G. Liu, _Lindquist-Wheeler formulation of lattice universes_ , _Phys. Rev. D_ 92 (2015) 063529, [1501.05169].
  * [42] V. A. A. Sanghai and T. Clifton, _Post-Newtonian Cosmological Modelling_ , _Phys. Rev. D_ 91 (2015) 103532, [1503.08747].
  * [43] P. Fleury, J. Larena and J.-P. Uzan, _The theory of stochastic cosmological lensing_ , _JCAP_ 11 (2015) 022, [1508.07903].
  * [44] E. Bentivegna, M. Korzyński, I. Hinder and D. Gerlicher, _Light propagation through black-hole lattices_ , _JCAP_ 03 (2017) 014, [1611.09275].
  * [45] V. A. A. Sanghai, P. Fleury and T. Clifton, _Ray tracing and Hubble diagrams in post-Newtonian cosmology_ , _JCAP_ 07 (2017) 028, [1705.02328].
  * [46] E. Bentivegna, T. Clifton, J. Durk, M. Korzyński and K. Rosquist, _Black-Hole Lattices as Cosmological Models_ , _Class. Quant. Grav._ 35 (2018) 175004, [1801.01083].
  * [47] P. Fleury, J. Larena and J.-P. Uzan, _Weak gravitational lensing of finite beams_ , _Phys. Rev. Lett._ 119 (2017) 191101, [1706.09383].
  * [48] P. Fleury, J. Larena and J.-P. Uzan, _Cosmic convergence and shear with extended sources_ , _Phys. Rev. D_ 99 (2019) 023525, [1809.03919].
  * [49] P. Fleury, J. Larena and J.-P. Uzan, _Weak lensing distortions beyond shear_ , _Phys. Rev. D_ 99 (2019) 023526, [1809.03924].
  * [50] T. Buchert and S. Räsänen, _Backreaction in late-time cosmology_ , _Ann. Rev. Nucl. Part. Sci._ 62 (2012) 57–79, [1112.5335].
  * [51] S. Räsänen, _Light propagation in statistically homogeneous and isotropic dust universes_ , _JCAP_ 02 (2009) 011, [0812.2872].
  * [52] S. Räsänen, _Light propagation in statistically homogeneous and isotropic universes with general matter content_ , _JCAP_ 03 (2010) 018, [0912.3370].
  * [53] G. F. R. Ellis, _Relativistic cosmology_ , _Proc. Int. Sch. Phys. Fermi_ 47 (1971) 104–182 (reprinted in Gen.Rel.Grav. 41 (2009) 581–660).
  * [54] C. G. Tsagas, A. Challinor and R. Maartens, _Relativistic cosmology and large-scale structure_ , _Phys. Rept._ 465 (2008) 61–147, [0705.4397].
  * [55] M. Lavinto, S. Räsänen and S. J. Szybka, _Average expansion rate and light propagation in a cosmological Tardis spacetime_ , _JCAP_ 12 (2013) 051, [1308.6731].
  * [56] P. Bull and T. Clifton, _Local and non-local measures of acceleration in cosmology_ , _Phys. Rev. D_ 85 (2012) 103512, [1203.4479].
  * [57] L. Verde, T. Treu and A. G. Riess, _Tensions between the Early and the Late Universe_ , _Nature Astron._ 3 (7, 2019) 891, [1907.10625].
  * [58] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , 1807.06209.
  * [59] A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn et al., _Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM_, _Astrophys. J. Lett._ 908 (2021) L6, [2012.08534].
  * [60] E. Joos and H. D. Zeh, _The Emergence of classical properties through interaction with the environment_ , _Z. Phys. B_ 59 (1985) 223–243.
  * [61] M. Tegmark, _Apparent wave function collapse caused by scattering_ , _Found. Phys. Lett._ 6 (1993) 571, [gr-qc/9310032].
  * [62] I. Allali and M. P. Hertzberg, _Gravitational Decoherence of Dark Matter_ , _JCAP_ 07 (2020) 056, [2005.12287].
  * [63] I. J. Allali and M. P. Hertzberg, _Decoherence from General Relativity_ , _Phys. Rev. D_ 103 (2021) 104053, [2012.12903].
  * [64] I. J. Allali and M. P. Hertzberg, _General Relativistic Decoherence with Applications to Dark Matter Detection_ , _Phys. Rev. Lett._ 127 (2021) 031301, [2103.15892].
  * [65] G. Ellis, R. Maartens and M. MacCallum, _Relativistic Cosmology_.  Relativistic Cosmology. Cambridge University Press, 2012.