# Sterile Neutrino Dark Matter from
Generalized $CPT$-Symmetric Early-Universe Cosmologies

Adam Duran<EMAIL_ADDRESS>Logan Morrison<EMAIL_ADDRESS>Stefano
Profumo<EMAIL_ADDRESS>Department of Physics and Santa Cruz Institute for
Particle Physics
University of California, Santa Cruz, CA 95064, USA

###### Abstract

We generalize gravitational particle production in a radiation-dominated
$CPT$-symmetric universe to non-standard, but also $CPT$-symmetric early
universe cosmologies. We calculate the mass of a right-handed “sterile”
neutrino needed for it to be the cosmological dark matter. Since generically
sterile neutrinos mix with the Standard Model active neutrinos, we use state-
of-the-art tools to compute the expected spectrum of gamma rays and high-
energy active neutrinos from ultra-heavy sterile neutrino dark matter decay.
We demonstrate that the sterile neutrinos are never in thermal equilibrium in
the early universe. We show that very high-energy Cherenkov telescopes might
detect a signal for sterile neutrino lifetimes up to around 1027 s, while a
signal in high-energy neutrino telescopes such as IceCube could be detectable
for lifetimes up to 1030 s, offering a better chance of detection across a
vast landscape of possible masses.

##  I Introduction

Dark matter (DM) continues to stand as one of the greatest mysteries at the
interface of particle physics and cosmology [1]. Since no particle in the
standard model (SM) can account for the bulk of the cosmological dark matter,
one must extend the particle content of the theory to include at least an
additional particle species, or a bound state thereof, that can account for a
matter species making up around five times the baryonic matter. Adding right-
handed neutrino (RHN) states is a well-motivated and oft-followed route that
can at once address the origin of active neutrino masses, explain the matter-
antimatter asymmetry, and provide a viable, arguably minimal, DM candidate
[2].

Recently, Ref. [3, 4] considered a peculiar and interesting possibility for
the generation of the cosmological dark matter: while in flat space-time the
“natural” vacuum is the one respecting the isometries of Minkowski space, in a
general, curved space-time this is no longer true. The choice of vacuum is
generally observer-dependent (for instance to inertial observers at different
points in the space-time), and different observers will define different,
inequivalent vacua. As a result, the zero-particle state according to one
observer can actually be populated by particles according to a different
observer (this is well known in the context of Hawking radiation from black
hole evaporation, or the Unruh radiation experienced by an accelerated
observer) [5, 6]. In particular, this holds generally for Friedmann-Robertson-
Walker (FRW) cosmologies. Ref. [3, 4] realized that assuming time-reversal
invariance allows to identify uniquely a $CPT$ invariant vacuum in a FRW
background, enabling a well-defined calculation of the number density of
gravitationally produced particles as seen by a late-time observer.

Gravitational particle production in an FRW background without inflation,
similar to the scenario under consideration here (albeit, here without the
additional requirement of $CPT$ invariance), has been studied extensively in
the past (for reviews, we refer the Reader to Ref. [7, 8]). The key
theoretical results are that particle production (specifically fermions, or
conformally coupled scalars) occurs at the cosmological epoch when particle
masses are comparable to the Hubble expansion rate; the particle number
density $n\sim m^{3}$, and particles are created with the equation of state of
dust; finally, the fractional relic density of these particles $X$ at the time
of radiation-matter equality was found to be $\Omega_{X}\sim(m_{X}/10^{9}\
{\rm GeV})^{5/2}$ [9]. Further studies considered similar scenarios, but in
the presence of inflation (see e.g. [10] and references therein).

If the early universe can be described by a conformally flat FRW metric, the
early-time and the late-time (or the ’in’ and ’out’) vacuum states are not
equivalent due to the lack of time-translation symmetry in the metric.
Calculating the number density of particles then involves performing a
Bogoliubov transformation on the quantum field modes and finding the
appropriate Bogoliubov coefficients. Remarkably, Ref. [3] found that the
standard in/out vacuum states do not respect $CPT$, and constructing a $CPT$
invariant vacuum state leads to a certain unique form of the relevant
Bogoliubov coefficient [3, 4]. Given the number density produced by this
expansion mechanism, if this species were to account for the total dark matter
(DM) density $\rho_{DM}$, then its mass would have to be around 4.8$\times$108
GeV. This calculation was performed in the context of an early universe
cosmology described by a radiation-dominated cosmology, with a scale factor
$a(\tau)\propto\tau$, where $\tau$ indicates the conformal time, and for the
most natural and minimal DM candidate, the right-handed neutrino, assumed to
be absolutely stable.

Here, we note that the requirement of $CPT$ symmetry is not unique to
radiation domination: as we explain in this work, the universe could have been
dominated, at early times, by some other species with a different energy
density, and yet satisfy $CPT$ symmetry. Note that the gravitational
production of particles discussed here does not touch upon the question of
whether the issues of homogeneity and isotropy are addressed or not, and only
relies on $CPT$ symmetry at positive conformal times [3, 4].

The possibility of an alternate cosmological history prior to Bug Bang
Nucleosynthesis has received considerable attention (for a recent review see
Ref. [11]), and here, as far as our motivation, we build largely on the
literature on the topic. The key observational fact for radiation domination
is the synthesis of light elements; prior to that, there is no observational
evidence that underpins radiation domination. In fact, the possibility of an
early matter-dominated phase is very well justified and has been widely
considered [12]; other causes of nonstandard early-universe expansion history
include single or multi-field post-inflation reheating , multi-step or thermal
inflation (in fact the presence of any extra stages of inflation, see e.g.
Ref. [13] for a review), heavy particles and dark sectors [14], moduli fields
[15]; the resulting structure of the cosmological expansion history impacts,
among other things, gravitational particle production, as was noted in a
recent study [16] that considered the case of quintessential cosmology in the
early universe, with a formalism not dissimilar to the one employed here.
Note, however, that the present setup is uniquely different from previous work
where gravitational production is connected to inflation.

A non-standard FRW background leads to different equations of motion for the
right-handed neutrino field, which, we note here, can be approximated using
the WKB method. Once the appropriate Bogoliubov coefficient is computed from
the WKB solutions in the asymptotic future, one can follow a similar analysis
to Ref. [3] and calculate the resulting particle number density at late times.
The final result of this exercise is a DM mass as a function of a parameter
which indicates the scaling of the energy density with the cosmological scale
factor and, as we discuss below, a second parameter that depends on when the
universe transitions from the non-standard early-times cosmology to radiation
domination, needed from cosmological observations to occur prior to the
synthesis of light elements [1].

In this work, besides generalizing the results of [3, 4] to any
$CPT$-symmetric early-universe cosmology, we relax a key assumption in the
dark matter sector: that the right-handed neutrino dark matter be absolutely
stable because of a discrete symmetry. Relaxing such assumption produces the
interesting prospect of testing experimentally this production mechanism by
searching for the decay products of massive right handed neutrinos, with a
lifetime related to the mixing angle between sterile and active neutrinos. In
turn, this calls for addressing the question of whether the gravitationally-
produced neutrinos can ever thermalize, which we address in detail in sec. IV.

The remainder of this work is structured as follows: in the following section
II we describe in detail the calculation of the dark matter number density
generated in arbitrary $CPT$ symmetric cosmologies in the early universe; the
following sec. III discusses the right-handed neutrino sector and the decay
modes thereof; sec. IV explores the issue of thermalization of right-handed
neutrinos; sec. V addresses constraints and prospects for the detection of
decaying massive right-handed neutrinos as predicted in the present scenario;
finally, sec. VI concludes.

##  II Sterile Neutrino as Dark Matter in Generalized $CPT$ Symmetric
Cosmologies

Consider a cosmology dominated by a species with equation of state $P=w\rho$,
and thus with an energy density-scale factor relation

$\rho\sim a^{-3(1+w)}.$ (1)

Recalling that, for $w\neq 1$, time $t\sim a^{3(1+w)/2}$, and that conformal
time is defined as

$\tau\equiv\int\frac{dt}{a},$ (2)

the relation between conformal time and scale factor for an equation of state
parameter $w$ is

$\tau\sim a^{\frac{1+3w}{2}},\quad{\rm or}\quad
a(\tau)\sim\tau^{\frac{2}{1+3w}}.$ (3)

The requirement for $CPT$ symmetry is that $a(\tau)=-a(-\tau)$; this, in turn,
requires

$\frac{2}{1+3w}={\rm odd}=N=2k+1,\ k=0,\pm 1,\pm 2,\dots$ (4)

As a result, we obtain the equation of state of radiation ($w=1/3$) for $N=1$,
$-1/3<w<-1/9$ equations of state for positive $N\neq 2$ (we disregard the
$w=0$ case), and $-1<w<-1/3$ for negative $N$. We do not concern ourselves
with the issue of, and the specific model for, how the transition between the
era when this species dominates in the early universe and radiation domination
occurs, since this does not yield any effect on what we concern ourselves with
here, i.e. gravitational particle production.

We hereafter generalize the particle production in a radiation-dominated
$CPT$-symmetric universe discussed in Ref. [4, 3] to the non-standard
cosmologies just described above. Ref. [4, 3] considered a heavy sterile
neutrino essentially decoupled from the thermal bath of the early universe and
constructed the corresponding quantum field coupled to the expanding
Friedmann-Robertson-Walker (FRW) metric

$ds^{2}=a^{2}(\tau)[-d\tau^{2}+dx^{2}]$ (5)

that is asymptotically flat at early ($\tau\rightarrow-\infty$) and late
($\tau\rightarrow+\infty$) times. The Lagrangian for the spinor field in the
expanding background is given by the usual

${\cal L}\supset i\bar{\psi}\not{\partial}\psi-\mu\bar{\psi}\psi,$ (6)

with equation of motion

$(i\not{\partial}-\mu)\psi=0.$ (7)

Following [3], we posit that the right handed neutrino mass term originates
from the vacuum expectation value of a field $\Phi$ whose equation of motion
at early times, prior to the electroweak phase transition, is simply

$\frac{1}{2}\left(\frac{d\Phi}{d\tau}\right)^{2}=a^{4}\rho=\sim
a^{4}a^{-3(1+w)}\sim\tau^{N(1-3w)}=\tau^{2N-2}.$ (8)

We call $\tau_{1}$ the conformal time at which the universe switches from
being dominated by a non-standard energy density $\rho\sim a^{-3(1+w)}\sim
a^{-2\frac{N+1}{N}}\sim\tau^{-2(N+1)}$ to radiation domination. Solving the
equation above gives

$\displaystyle\derivative{\Phi}{\tau}=\sqrt{2\rho}a^{2}\sim\tau^{N-1}=\begin{cases}C_{1}\tau^{N-1}&0<\tau<\tau_{1}\\\
C_{2}&\tau>\tau_{1}\\\ \end{cases}$ (9)

and thus:

$\displaystyle\Phi(\tau)=\begin{cases}C_{1}\frac{\tau^{N}}{N},&0<\tau<\tau_{1}\\\
C_{1}\frac{\tau_{1}^{N}}{N}+C_{2}\quantity(\tau-\tau_{1}),&\tau>\tau_{1}\\\
\end{cases}$ (10)

and hence, the effective mass is:

$\displaystyle\mu(\tau)=y\Phi(\tau)=y\times\begin{cases}C_{1}\frac{\tau^{N}}{N},&0<\tau<\tau_{1}\\\
C_{1}\frac{\tau_{1}^{N}}{N}+C_{2}\quantity(\tau-\tau_{1}),&\tau>\tau_{1}\\\
\end{cases}$ (11)

To find the constants in the equation above, first ask that the asymptotic
value in the large $\tau$ limit give the “right” radiation domination answer,
i.e. $\mu(\tau)=y(2\rho_{1})^{1/2}\tau=\gamma\tau$, with $\rho_{1}$ the value
at the “beginning” of radiation domination, i.e. for us at conformal time
$\tau_{1}$, and $y=M/\hat{\mu}$: this sets $C_{2}=1$; second, demand
continuity of both $\Phi(\tau)$ (no constraints) and of its derivative, which
imposes $C_{1}\tau_{1}^{N-1}=C_{2}=1$ thus $C_{1}=1/\tau_{1}^{N-1}$. So indeed
for large $\tau<\tau_{1}$, where we hypothesize production occurs,
$\mu(\tau)=\gamma\tau^{N}/\tau_{1}^{N-1}$.

Given the symmetries associated with the FRW background, different observers
disagree on what is a positive- or a negative-frequency mode: each observer
will have their own corresponding mode functions, operators, and vacuum
states. Ref. [4, 3] label these the ingoing ‘$-$’ and outgoing ‘+’ states, for
observers in the asymptotic, conformally invariant past, and in the non-
conformally invariant late-time universe. The fields expansion for these modes
are given by

$\psi(x)=\int\frac{d^{3}p}{(2\pi)^{3}}\bigg{[}a_{\pm}(p,h)u_{\pm}(p,h,\tau)e^{ipx}+b_{\pm}^{\dagger}(p,h)\upsilon_{\pm}(p,h,\tau)e^{-ipx}\bigg{]}.$
(12)

In the equation above, the operators $a_{\pm}(p,h)$ and $b_{\pm}(p,h)$ act on
their respective vacuum states of the asymptotic observers, which are defined
by

$\displaystyle a_{-}(p,h)\ket{0_{-}}=b_{-}(p,h)\ket{0_{-}}=0,$ (13)
$\displaystyle a_{+}(p,h)\ket{0_{+}}=b_{+}(p,h)\ket{0_{+}}=0.$ (14)

In order to relate the field expansions and operators of the two asymptotic
observers, one employs the Bogoliubov transformation. The key result of this
transformation is that mode functions and operators for early-time observers
can be expressed as a linear combination of the late-time observer’s mode
functions and operators, which form a complete set of states, as

$\displaystyle
u_{-}(p,h,\tau)=\alpha(p,h)u_{+}(p,h,\tau)+\beta(p,h)\upsilon_{+}(p,h,\tau),$
(15) $\displaystyle
a_{-}(p,h)=\alpha(p,h)a_{+}(p,h)+\beta(p,h)a_{+}^{\dagger}(p,h),$ (16)

where $\alpha(p,h)$ and $\beta(p,h)$ are the so-called Bogoliubov
coefficients. If we apply the number operator of the ingoing observer
$N_{-}=a^{\dagger}_{-}a_{-}$ to the outgoing observers vacuum state, the
result is

$\braket{0_{+}}{N_{-}}{0_{+}}=|\beta(p,h)|^{2}.$ (17)

Integrating over the phase space gives the mean particle number density at
late times, one gets the desired result

$n=2\sum_{h}\int\frac{d^{3}p}{(2\pi)^{3}}\mid\beta(p,h)\mid^{2}.$ (18)

The equation above implies that particle number is observer-dependent, which
is a well-known fact in the context of the theory of quantum fields in
expanding backgrounds. However, Ref. [4, 3] found that since the ingoing and
outgoing vacuum states are not equivalent, they are not $CPT$ invariant (See
section 3.3 in [3]). To correct this, Ref. [4, 3] constructed a $CPT$
invariant vacuum state that respects the transformations of the field. A
Bogoliubov transformation matrix considering $CPT$ transformations was
calculated along with the corresponding coefficients (See section 3.4 in [3]).
The coefficient of interest $\beta_{+}(p,h)$ was calculated to be

$\beta_{+}(p,h)=\sin{\bigg{[}\frac{1}{2}\arcsin{\big{(}\beta(p,h)\big{)}}\bigg{]}},$
(19)

where this quantity gives the mean particle density given the assumed $CPT$
invariant vacuum state of the early universe. The Bogoliubov coefficient
$\beta(p)$ in the inverse sine is calculated considering the ingoing and
outgoing vacuum states.

In the context of the generalized conformally invariant cosmologies we
consider here, the Bogoliubov coefficient for the ingoing and outgoing states
is corrected by an effective mass that now reads
$\mu(\tau)=\gamma\tau^{N}/\tau_{1}^{N-1}$. Since the fields mode functions
form an orthonormal basis, one can follow the normalization condition

$\sum_{h}\big{[}u_{\pm}(p,h,\tau)u_{\pm}^{*}(p,h,\tau)+\upsilon_{\pm}(p,h,\tau)\upsilon_{\pm}^{*}(p,h,\tau)\big{]}=\mathds{I}_{\textrm{4x4}};$
(20)

the Bogoliubov coefficient $\beta(p)$ in (19) can be then calculated using the
inner product of ingoing and outgoing mode functions,

$\braket{\upsilon_{-}}{u_{+}}=\sum_{h}\big{[}\upsilon_{-}^{*}(p,h,\tau)u_{+}(p,h,\tau)\big{]}=\beta(p).$
(21)

The next task is to solve the equation of motion (7) and find the asymptotic
form. For a free Weyl invariant scalar field defined on an expanding FRW
metric with an effective mass $\mu$ the Lagrangian is given by

$\displaystyle
L=(\varphi^{\prime})^{\dagger}(\varphi^{\prime})-(\nabla\varphi)^{\dagger}(\nabla\varphi)-\mu^{2}\varphi^{\dagger}\varphi$
(22)

with the equation of motion

$(\partial_{\tau}^{2}-\nabla^{2}+\mu^{2})\varphi=0.$ (23)

Taking the solutions to momentum space: $\varphi(x)=u(p,\tau)e^{ipx}$ turns
the equation of motion into that of a harmonic oscillator equation with
frequency $\omega(\tau)$,

$\displaystyle
u^{\prime\prime}(p,\tau)+\omega^{2}(\tau)u(p,\tau)=0\quad;\quad\omega^{2}(\tau)=p^{2}+\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}.$
(24)

To see how the mode solutions vary with conformal time $\tau$, we numerically
solved the equation for various $N$ values, and compared with an approximate
solution from the WKB formalism [17]. For the ingoing and outgoing observers
the WKB solutions are given by

$u_{\pm}(p,\tau)=\frac{1}{\sqrt{\pm 2\omega(p,\tau)}}\exp[\mp
i\int_{\pm\tau_{0}}^{\tau}\omega(p,\tau^{{}^{\prime}})d\tau^{{}^{\prime}}\bigg{]}$
(25)

where $\tau_{0}$ is arbitrary and can be fixed for convenience. The integral
can be carried out analytically, and the result is given by

$\int...=\frac{\tau}{N+1}\bigg{(}\sqrt{p^{2}+\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}}+\\\
pN_{2}F_{1}\bigg{(}\frac{1}{2},\frac{1}{2N},1+\frac{1}{2N};-\frac{\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}}{p^{2}}\bigg{)}\bigg{)}\bigg{|}_{\pm\tau_{0}}^{\tau}.$
(26)

This expression simplifies considerably in the large $\tau$ limit. Here, we
assume that particle production occurs at conformal times smaller than the
time $\tau_{1}$ at which the universe becomes radiation dominated. Thus, the
large $\tau$ limit here means that $\tau_{1}yT_{1}\gg 1$, i.e.
$\tau_{1}T_{1}\gg\frac{1}{y}=\frac{M_{P}}{M}$, with $M$ the right-handed
neutrino mass.

Using

$\frac{T}{T_{\rm BBN}}=\frac{a_{\rm BBN}}{a}\longrightarrow T=\frac{T_{\rm
BBN}a_{\rm BBN}}{a},$

we have

$\tau_{1}T_{1}=\frac{T_{\rm BBN}a_{\rm
BBN}}{a_{1}}\int_{0}^{a_{1}}\frac{da}{a^{2}H(a)},$

where $H(a)\simeq H_{0}\sqrt{\Omega_{w}}a^{-(N+1)/N}$ at $a<a_{1}$ by
assumption, where $\Omega_{w}$ is the relative energy density of the exotic
component,

$\Omega_{r}a_{1}^{-4}=\Omega_{w}a^{-(2N+2)/N}\longrightarrow\Omega_{w}=\Omega_{r}a_{1}^{(2-2N)/N}.$

We thus find

$\tau_{1}T_{1}=\frac{T_{\rm BBN}a_{\rm
BBN}}{a_{1}H_{0}\sqrt{\Omega_{r}}a_{1}^{-(N+1)/N}}\int_{0}^{a_{1}}\frac{da}{a^{2}a^{-(N+1)/N}}=\frac{T_{\rm
BBN}a_{\rm BBN}}{H_{0}\sqrt{\Omega_{r}}a_{1}^{1/N}}\left(Na_{1}^{1/N}\right).$

Finally, we find that the condition does not depend on $a_{1}$, and is always
satisfied since, putting in the numbers (with $H_{0}\sim 10^{-42}$ GeV and
$\Omega_{r}\sim 5\times 10^{-5}$), we have

$\frac{M}{M_{P}}\gg\frac{H_{0}\sqrt{\Omega_{r}}}{NT_{\rm BBN}a_{\rm
BBN}N}\simeq\frac{4\times 10^{-32}}{N},$

implying $M\gg 10^{-5}$ eV. In order to be the dark matter, the RHN must be
heavier than at least a few keV, so the large-$\tau$ condition is always
satisfied.

In the large $\tau$ limit, the first term in Eq. (26) reads

$\sqrt{p^{2}+\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}}\approx\gamma\tau^{N}/\tau_{1}^{N-1}.$
(27)

Upon integration, one gets for that first term:

$\frac{\gamma}{N+1}(\tau^{N+1}/\tau_{1}^{N-1}-(\pm\tau_{0})^{N+1}/\tau_{1}^{N-1}).$
(28)

For the hypergeometric function in Eq. (26), in the same limit, one can make
the replacement

${}_{2}F_{1}\bigg{(}\frac{1}{2},\frac{1}{2N},1+\frac{1}{2N};-\frac{\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}}{p^{2}}\bigg{)}\Longrightarrow\
_{2}F_{1}\bigg{(}\frac{1}{2},\frac{1}{2},\frac{3}{2};-\frac{\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}}{p^{2}}\bigg{)}$

and

${}_{2}F_{1}\bigg{(}\frac{1}{2},\frac{1}{2},\frac{3}{2};-\frac{\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}}{p^{2}}\bigg{)}=\frac{p\tau_{1}^{N-1}}{\gamma\tau^{N}}\ln(\gamma\tau^{N}/\tau_{1}^{N-1}+\sqrt{p^{2}+\gamma^{2}\tau^{2N}/\tau_{1}^{2N-2}})$

The approximate asymptotic behavior of $u_{\pm}(p,h,\tau)$ and
$\upsilon_{\pm}(p,h,\tau)$ in the large $\tau$ limit is then found to be

$u_{\pm}(p\hat{z},h,\tau)\simeq\frac{1}{\sqrt{2}}\exp[\mp
i\chi\bigg{]}\begin{bmatrix}\pm 1\\\ 0\\\ +1\\\ 0\\\
\end{bmatrix};\quad\upsilon_{\pm}(p\hat{z},h,\tau)\simeq\frac{i}{\sqrt{2}}\exp[\pm
i\chi\bigg{]}\begin{bmatrix}-1\\\ 0\\\ \pm 1\\\ 0\\\ \end{bmatrix},$ (29)

where, in the equation above,

$\chi=\bigg{(}\frac{\gamma^{\prime}}{N+1}\\{\tau^{N+1}-\tau_{0}^{N+1}\\}+\frac{p^{2}N}{\gamma^{\prime}(N+1)}\bigg{(}\frac{1}{\tau^{N-1}}\ln{\tau^{N}}-\frac{1}{\tau_{0}^{N-1}}\ln{\pm\tau_{0}^{N}}\bigg{)}\bigg{)}.$

For convenience, above we indicate with
$\gamma^{\prime}=\frac{\gamma}{\tau_{1}^{N-1}}$. Plugging this expression into
the expressions above, we get the Bogoliubov coefficient $\beta(p)$ as

$\beta(p)=-i\exp[-\frac{N\pi
p^{2}}{\gamma(N+1)}\left(\frac{1}{\xi}\right)^{N-1}\bigg{]},$ (30)

where we introduced the quantity $\xi\equiv\tau_{0}/\tau_{1}$. Notice that
$\xi$ depends on the (in principle arbitrary, as long as preceding Big Bang
Nucleosynthesis) conformal time at which the universe transitions to being
radiation-dominated, and on the parameter $\tau_{0}$, which is chosen so that
Eq. (29) is valid (see [3] for details on this). Plugging this into (19) gives
the generalized number density

$n_{dm}=\sum_{h}\int\frac{d^{3}p}{(2\pi)^{3}}\mid\beta_{+}(p)\mid^{2}=I\bigg{(}\frac{(N+1)\gamma}{2\pi
N}\xi^{(N-1)}\bigg{)}^{3/2}.$ (31)

In the equation above, $I$ is a dimensionless constant defined as

$I=\frac{1}{2\pi^{2}}\int_{0}^{\infty}x^{2}\bigg{[}1-\sqrt{1-e^{-x^{2}}}\bigg{]}dx\approx
0.01276.$

Our result recovers the number density found in [3] for the radiation
dominated background, $N=1$. To extract a prediction for the dark matter
(right handed neutrino) mass out of this number density, we follow the same
procedure as in [3] (See their section 5.1.2). With $n_{dm}$, we define the
dark matter yield

$Y_{dm}=\frac{n_{dm}}{s},$ (32)

with $s$ the entropy density in the universe, which we assume to be conserved
at early times and in particular at the transition to radiation domination.

After the reference conformal time $\tau_{1}$, we have the usual radiation-
domination relations

$\displaystyle\rho=\frac{\pi^{2}}{30}g_{*}T^{4},\qquad
s=\frac{2\pi^{2}}{45}g_{*}T^{3}.$ (33)

From these relations, solving for temperature, one finds

$\frac{(2\rho)^{3/4}}{s}=\frac{3}{2}\bigg{(}\frac{15}{g_{*}\pi^{2}}\bigg{)}^{1/4}.$
(34)

Combining all expressions for $n_{dm}$, $\gamma$ and $s$, we can express
$Y_{dm}$ in terms of the particle mass $M_{dm}$

$Y_{dm}=\frac{3I}{2\pi^{2}}\bigg{(}\frac{15}{g_{*}}\bigg{)}^{1/4}\bigg{(}\frac{M_{\nu}}{\mu}\bigg{)}^{3/2}\bigg{(}\frac{N+1}{2N}\xi^{(N-1)}\bigg{)}^{3/2},$
(35)

where $g_{*}=106.75$ is the number of the degrees of freedom in the standard
model, and $\hat{\mu}=5.19\times 10^{8}$GeV is the same right-handed neutrino
reference mass scale chosen by [3]. Notice that the formula, remarkably, does
not depend on the conformal time $\tau_{1}$ at which radiation domination sets
in (there could have in principle been a dependence both in the dark matter
number density and via the energy density $\rho_{1}$).

We can now match the mass of the dark matter particle and the number density
to the present day dark matter energy density $\rho_{dm}^{(0)}=9.7\times
10^{-48}$GeV4 and the present day entropy density $s^{(0)}=2.3\times
10^{-38}$GeV3

The late-time dark matter density is then given by

$\rho_{dm}^{(0)}=M_{dm}n_{dm}=M_{dm}Y_{dm}s^{(0)}.$ (36)

Solving for the mass term for $M_{dm}=M_{\nu}$ gives the desired result for
the sterile neutrino DM mass as function of $N$ The result is a sterile
neutrino mass as a function of integer $N$ values.

$m_{\bar{\nu}}=4.8\times
10^{8}\bigg{(}\frac{2N}{N+1}\bigg{)}^{3/5}\xi^{-\frac{3(N-1)}{5}}\textrm{GeV}\qquad({\rm
odd}\ N).$ (37)

Figure 1: Mass of the right-handed neutrino giving the correct dark matter
density as a function of $\xi$ for various values of $N$.

Fig. 1 shows, as a function of the parameter $\xi$, the mass of the right-
handed neutrino $\bar{\nu}$ giving the right dark matter amount for a given
value of the parameter $N$.

This setup does not imply values for parameters such as the Majorana and Dirac
masses of the neutrino sectors, and thus the lifetime of the dark matter
particle is not a prediction of the model. We describe ways to detect sterile
neutrino decay products in this scenario in the following section, in a model-
independent way.

##  III Sterile Neutrino Decay Width

For simplicity, we assume the existence of a single, Weyl, sterile neutrino
$\bar{\nu}$. We take the sterile neutrino/SM interactions to stem from a
Yukawa interaction between the RH neutrino and the SM lepton doublet. In this
scenario, the Lagrangian density is given by:

$\displaystyle\mathcal{L}=\mathcal{L}_{\mathrm{SM}}+i\hat{\bar{\nu}}^{\dagger}\bar{\sigma}_{\mu}\partial^{\mu}\hat{\bar{\nu}}-\frac{1}{2}\quantity(\hat{m}\hat{\bar{\nu}}\hat{\bar{\nu}}+\mathrm{h.c.})+\quantity(\epsilon^{ab}\mathbf{y}^{i}_{\nu}\phi_{a}\hat{L}_{bi}\hat{\bar{\nu}}+\mathrm{h.c.}).$
(38)

In the expression above, $\hat{m}$ is the Majorana mass of the RH neutrino,
$\mathbf{y}_{\nu}$ is a Yukawa vector with three components and $\phi$, and
$\hat{L}$ are the Higgs and lepton doublets (note that all the above fermion
fields are written in two-component notation). The doublets are given by:

$\displaystyle\phi$ $\displaystyle=\matrixquantity(G^{+}\\\
\frac{1}{\sqrt{2}}\quantity(v+h+iG^{0})),$ $\displaystyle\hat{L}_{i}$
$\displaystyle=\matrixquantity(\hat{\nu}_{i}\\\ \hat{\ell}_{i})$ (39)

with $G^{+}$ and $G^{0}$ being the charged and neutral Goldstone bosons, $h$
being the standard model Higgs, $v\sim 246$ GeV being the Higgs vacuum
expectation value, $\hat{\nu}_{i}$ being the left-handed neutrinos and
$\hat{\ell}_{i}$ being the left-handed charged leptons. Expanding out the
neutrino Lagrangian and gathering all the mass terms, we find the following:

$\displaystyle-\mathcal{L}_{\nu\mathrm{-mass}}$
$\displaystyle=\frac{1}{2}\hat{m}\hat{\bar{\nu}}\hat{\bar{\nu}}+\frac{v}{\sqrt{2}}\mathbf{y}^{i}_{\nu}\hat{L}_{bi}\hat{\bar{\nu}}+\mathrm{h.c.}=\frac{1}{2}\mathbf{M}_{\nu}^{ij}\hat{{\bf\it\nu}}_{i}\hat{{\bf\it\nu}}_{j}$
(40)

where the neutrino vector
$\hat{{\bf\it\nu}}=\matrixquantity(\hat{\nu}_{1}&\hat{\nu}_{2}&\hat{\nu}_{3}&\hat{\bar{\nu}})^{T}$
and the $(4\times 4)$ neutrino mass matrix is given by:

Diagonalizing the neutrino mass matrix is done through Takagi diagonalization
with a unitary matrix ${\bf\it\Omega}$, satisfying
${\bf\it\Omega}^{T}\mathbf{M}_{\nu}{\bf\it\Omega}=\mathrm{diag}(m_{\nu_{e}},m_{\nu_{\mu}},m_{\nu_{\tau}},m_{\bar{\nu}})$.
While it is possible to analytically compute ${\bf\it\Omega}$ for general
$\mathbf{y}_{\nu}$, we choose to simplify our model by assuming only a single
non-zero Yukawa coupling. That is, we take $\mathbf{y}^{k}_{\nu}\equiv y$ and
$\mathbf{y}_{\nu}^{k\neq 0}=0$ with $k=1,2$ or $3$. in this case, the unitary
Takagi matrix is given by an orthogonal $\mathrm{O}(2)$ matrix times a
diagonal unitary matrix:

The angle satisfies

$\displaystyle\cos\theta$
$\displaystyle=\frac{\sqrt{\hat{m}^{2}+2v^{2}y^{2}}+\hat{m}}{\sqrt{\quantity(\hat{m}+\sqrt{\hat{m}^{2}+2v^{2}y^{2}})^{2}+2v^{2}y^{2}}}=\sqrt{\frac{m_{\bar{\nu}}}{m_{\nu}+m_{\bar{\nu}}}}$
(41) $\displaystyle\sin\theta$
$\displaystyle=\frac{\sqrt{\hat{m}^{2}+2v^{2}y^{2}}-\hat{m}}{\sqrt{\quantity(\hat{m}-\sqrt{\hat{m}^{2}+2v^{2}y^{2}})^{2}+2v^{2}y^{2}}}=\sqrt{\frac{m_{\nu}}{m_{\nu}+m_{\bar{\nu}}}}$
(42)

with the light and heavy neutrino masses $m_{\nu}$ and $m_{\bar{\nu}}$ are
given by:

$\displaystyle m_{\nu}$
$\displaystyle=\frac{1}{2}\quantity(\sqrt{\hat{m}^{2}+2v^{2}y^{2}}-\hat{m})\sim\hat{m}\quantity(\frac{v^{2}y^{2}}{2\hat{m}^{2}}+\order{\frac{vy}{\hat{m}}}^{4})$
(43) $\displaystyle m_{\bar{\nu}}$
$\displaystyle=\frac{1}{2}\quantity(\sqrt{\hat{m}^{2}+2v^{2}y^{2}}+\hat{m})\sim\hat{m}\quantity(1+\frac{v^{2}y^{2}}{2\hat{m}^{2}}+\order{\frac{vy}{\hat{m}}}^{4})$
(44)

The mass eigenstates for the LH and RH neutrinos, $\nu_{k}$ and $\bar{\nu}$
are related to the interaction eigenstates via
${\bf\it\Omega}\matrixquantity(\nu_{k}&\bar{\nu})^{T}$, or more explicitly:

$\displaystyle\hat{\nu}_{k}$
$\displaystyle=-i\cos\theta\nu_{k}+\sin\theta\bar{\nu}$ (45)
$\displaystyle\hat{\bar{\nu}}$
$\displaystyle=i\sin\theta\nu_{k}+\cos\theta\bar{\nu}$ (46)

Using these transformations, we find that the Lagrangian density containing
the interactions between the RH neutrino and the standard model relavent to
this study is given by:

$\displaystyle\mathcal{L}$
$\displaystyle=\cdots+\frac{g}{\sqrt{2}}\quantity[i\cos\theta
W^{+}_{\mu}\nu_{k}^{\dagger}\bar{\sigma}_{\mu}\ell_{k}+\sin\theta
W^{+}_{\mu}\bar{\nu}^{\dagger}\bar{\sigma}_{\mu}\ell_{k}+\mathrm{h.c.}]$ (47)
$\displaystyle\hskip
29.87547pt+\frac{e}{2s_{W}c_{W}}\quantity[i\sin\theta\cos\theta
Z_{\mu}\nu^{\dagger}_{k}\bar{\sigma}_{\mu}\bar{\nu}+\mathrm{h.c.}]$
$\displaystyle\hskip 29.87547pt+\frac{iy}{\sqrt{2}}h\nu_{k}\bar{\nu}$

where the $\cdots$ are terms not containing interactions between the RH
neutrino and the SM or they contain irrelevant vertices (for example,
interactions with the Goldstones.) The partial decay widths for $\bar{\nu}\to
W^{+}+\ell,W^{-}+\ell^{\dagger},Z+\nu$ and $h+\nu$ are given by:

$\displaystyle\Gamma(N\to\nu_{\ell}+h)$
$\displaystyle=\dfrac{\lambda^{1/2}\quantity(m_{\bar{\nu}}^{2},m_{\nu}^{2},m_{h}^{2})}{16\pi
m_{\bar{\nu}}v_{H}^{2}}\sin[2](\theta)\quantity[(m_{\bar{\nu}}-m_{\nu})^{2}-m_{h}^{2}]$
(48) $\displaystyle\Gamma(N\to W^{+}+\ell)$
$\displaystyle=\dfrac{e^{2}\sin[2](\theta)}{64\pi
m_{\bar{\nu}}^{3}M_{W}^{2}s_{W}^{2}}\lambda^{1/2}\quantity(m_{\bar{\nu}}^{2},m_{\ell}^{2},M_{W}^{2})$
(49)
$\displaystyle\qquad\times\quantity[\quantity(m_{\bar{\nu}}^{2}-m_{\ell}^{2})^{2}+M_{W}^{2}\quantity(m_{\bar{\nu}}^{2}+m_{\ell}^{2})-2M_{W}^{2}]$
$\displaystyle\Gamma(N\to Z+\nu_{\ell})$
$\displaystyle=\dfrac{e^{2}\sin[2](2\theta)}{256\pi
m_{\bar{\nu}}^{3}M_{W}^{2}s_{W}^{2}}\lambda^{1/2}\quantity(m_{\bar{\nu}}^{2},m_{\nu}^{2},M_{Z}^{2})$
(50)
$\displaystyle\qquad\times\quantity(m_{\bar{\nu}}-M_{Z}-m_{\nu})\quantity(m_{N}+M_{Z}-m_{\nu})\quantity(2M_{Z}^{2}+\quantity(m_{\bar{\nu}}+m_{\nu})^{2})$

where $\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$. Taking $m_{\bar{\nu}}\ll
m_{H}$, we find:

$\displaystyle\Gamma(N\to\nu_{\ell}+h)$
$\displaystyle\sim\dfrac{m_{\bar{\nu}}}{16\pi}\quantity(\frac{m_{\bar{\nu}}}{v_{H}})^{2}\sin[2](\theta)$
(51) $\displaystyle\Gamma(N\to W^{+}+\ell)$
$\displaystyle\sim\dfrac{e^{2}m_{\bar{\nu}}}{64\pi
s_{W}^{2}}\quantity(\frac{m_{\bar{\nu}}}{m_{W}})^{2}\sin[2](\theta)$ (52)
$\displaystyle\Gamma(N\to Z+\nu_{\ell})$
$\displaystyle\sim\dfrac{e^{2}m_{\bar{\nu}}}{256\pi
s_{W}^{2}}\quantity(\frac{m_{\bar{\nu}}}{m_{W}})^{2}\sin[2](2\theta)$ (53)

The total decay with of the RH neutrino is roughly (for $\theta\ll 1$
$m_{\bar{\nu}}\gg m_{H}$):

$\displaystyle\Gamma\sim\frac{10^{-6}}{\mathrm{GeV}^{2}}\theta^{2}m_{\bar{\nu}}^{3}$
(54)

and the corresponding lifetime

$\tau=\frac{6.6\times 10^{-28}\ {\rm s}}{\theta^{2}}\left(\frac{\rm 1000\
GeV}{m_{\bar{\nu}}}\right)^{3},$ (55)

meaning that

$\frac{\tau}{\tau_{U}}\simeq
15\left(\frac{10^{-23}}{\theta}\right)^{2}\left(\frac{\rm 1000\
GeV}{m_{\bar{\nu}}}\right)^{3},$ (56)

where we indicate with $\tau_{U}$ the age of the universe.

##  IV Sterile Neutrino Thermalization

A necessary condition for the validity of the dark matter mass predictions
discussed above is that the right handed neutrinos not thermalize: should they
reach chemical equilibrium with the Standard Model thermal bath, the
corresponding thermal relic density would be unrelated to the abundance from
gravitational production. We must thus ensure that the parameter space under
consideration be consistent with the absence of thermal equilibrium for the
neutrino.

At $T>T_{1}$, where $T_{1}$ is the temperature at which the energy density of
radiation equals the energy density that dominates in the early universe
during gravitational dark matter production,

$\rho\sim T^{2\frac{N+1}{N}}\qquad(T>T_{1}).$

Matching $\rho_{\rm rad}(T_{1})$ with the expression above gives

$\rho\simeq\left(\frac{T}{T_{1}}\right)^{2\frac{N+1}{N}}T_{1}^{4}.\qquad(T>T_{1}).$

Right handed neutrinos can be in thermal equilibrium only when relativistic,
which allows us to express the corresponding number density as
$n_{\bar{\nu}}\sim T^{3}$ ($T>m_{\bar{\nu}})$. At large-enough temperatures,
larger than both $m_{\bar{\nu}}$ and the electroweak scale, the cross section
responsible for neutrino thermalization is approximately
$\sigma_{\bar{\nu}-{\rm SM}}\sim\theta^{2}/T^{2}$, from decay and inverse
decay into e.g. a gauge boson and a lepton.

In the radiation domination phase ($T<T_{1}$), when the Hubble rate is
approximately $H\sim T^{2}/M_{P}$, thermalization happens at a temperature
such that

$T^{3}\cdot\frac{\theta^{2}}{T^{2}}\sim\frac{T^{2}}{M_{P}}\quad\rightarrow\quad
T\sim M_{P}\theta^{2}.$ (57)

Since by assumption $T<m_{\bar{\nu}}$, a sufficient condition to prevent
thermalization is that $m_{\bar{\nu}}<M_{P}\theta^{2}$. From the discussion
above, Eq. (56), for TeV and heavier right handed neutrinos,
$M_{P}\theta^{2}\gg 10^{-28}$ GeV, so thermalization can never occur during
radiation domination.

In the non-standard phase, $T>T_{1}$, the Hubble rate is approximately

$H(T)\sim\frac{T_{1}^{2}}{M_{P}}\left(\frac{T}{T_{1}}\right)^{\frac{N+1}{N}}\geq\frac{T_{1}^{2}}{M_{P}}\qquad(T>T_{1}),$
(58)

thus to ensure that thermal equilibrium never be attained we simply need to
require that

$n_{\bar{\nu}}\sigma_{\bar{\nu}-{\rm
SM}}\sim{T\theta^{2}}\ll\frac{T_{1}^{2}}{M_{P}}\leq H(T)$ (59)

Below the Planck scale, $n_{\bar{\nu}}\sigma_{\bar{\nu}-{\rm
SM}}<{M_{P}\theta^{2}}$, implying that equilibrium can only be attained for
$T_{1}\sim\theta M_{P}$ which is much smaller than the scale of Big Bang
Nucleosynthesis $T_{\rm BBN}\sim 1$ MeV, at which the universe must be
radiation dominated (in other words, $T_{1}>T_{\rm BBN}$. Thus, thermal
equilibrium is never attained at any point below the Planck scale.

##  V Constraints on Sterile Neutrino Lifetime

Tight constraints on the sterile neutrino lifetime can be computed using
experiments such as the Fermi Large Area Telescope (LAT) [18], IceCube [19],
and HAWC [20].

The gamma-ray and neutrino differential fluxes detectable on Earth from the
decay of a sterile neutrino are given by:

$\displaystyle\derivative{\phi_{\gamma,\nu}}{E_{\gamma,\nu}}$
$\displaystyle=\frac{1}{m_{\bar{\nu}}\tau}\derivative{N_{\gamma,\nu}}{E_{\gamma,\nu}}\frac{J}{4\pi}$
(60)

where $\tau$ is the sterile neutrino life-time,
$\differential{N}_{\gamma,\nu}/\differential{E}_{\gamma,\nu}$ is the gamma-
ray/neutrino spectrum per decay and $J$ is the so-called $J$-factor, the
integral along the line of sight, and averaged on the telescope’s angular
aperture, of the dark matter density, which depends on the target being
observed. Unlike for annihilation, the calculation of the $J$ factor for decay
is prone to significantly less severe uncertainties from the slope of the
inner density of the dark matter profile. Here, we use the same $J$ factors as
in [21, 22].

Figure 2: The gamma-ray (left) and neutrino spectra from the decay of a
sterile neutrino with a mass of $4.8\times 10^{8}$ GeV (corresponding to the
radiation domination case).

Computing the gamma-ray/neutrino spectrum for a sterile neutrino mass
$m_{\bar{\nu}}\sim 10^{8}$ GeV is a non-trivial task. Since $m_{\bar{\nu}}\gg
m_{W},m_{Z},m_{H}$, the electro-weak bosons emitted from the sterile neutrino
decays are highly-relativistic and are prone to undergo large numbers of
splittings into more electro-weak bosons. For this reason, we use a newly
developed software package, HDMSpectra, specifically developed for computing
spectra from the decay of very heavy dark matter ($m_{\mathrm{DM}}\gg m_{H}$)
[23]. Fig. 2 shows the gamma-ray and neutrino spectra from the decay of a
sterile neutrino with a mass of $m_{\bar{\nu}}\sim 4.8\times 10^{8}$ GeV.

To compute the constraints from HAWC data, we use the differential fluxes
provided in Ref. [20]. In that study, the HAWC collaboration measured the
gamma-ray fluxed from the Andromeda galaxy, an ideal target to look for dark
matter decay debris. They measured fluxes of photons with energies ranging
from a TeV to 100 TeV and were able to place constraints on
$E^{2}\differential{N}_{\gamma}/\differential{E_{\gamma}}\lesssim
10^{-12}-10^{-11}\ \mathrm{TeV}\ \mathrm{cm}^{-2}\ \mathrm{s}^{-1}$. Using
these results, we find a constraint of the sterile neutrino lifetime of
roughly $\tau\gtrsim 7\times 10^{24}\ \mathrm{s}$ (see Fig. (3)).

Figure 3: Upper limits on the right handed neutrino lifetime as a function of
the particle mass, from IceCube [19] (blue line), HAWC [20] (red line), and
Fermi LAT [18] (yellow line).

To calculate constraints for IceCube, we use the Galactic center neutrino
upper limit fluxes shown in Fig. 2 of Ref. [19]; similarly, for the
constraints from the Fermi Large Area Telescope (LAT) [18], we use the
sensitivity discussed in Ref. [18].

We show the upper limits to the lifetime of a decaying right-handed neutrino
in fig. 3. By far the best limits arise from the neutrino telescope IceCube,
which in the mass range of typical interest, between $10^{3}$ and $10^{8}$
GeV, achieves a sensitivity in excess of $10^{29}$ s, two to three orders of
magnitude better than gamma-ray telescopes. This lifetime sensitivity
translates into probing mixing angles as small as $\theta\sim 10^{-59}$.

##  VI Summary, Discussion, and Conclusions

We generalized gravitational particle production in the early universe in
$CPT$ symmetric cosmologies beyond the case of radiation domination discussed
in Ref. [3, 4]; our key result is that in non-standard, $CPT$ symmetric
universes the mass of the dark matter candidate, a “sterile” right-handed
neutrino, depends on the transition from the non-standard cosmological epoch
to radiation, and on the dependence of the energy density on conformal time at
early times. As a result, unlike the case of radiation domination, the dark
matter mass in largely unconstrained and could range between the electroweak
and the Planck scales. Note that we did not consider cases where the expansion
history of the universe, while respecting $CPT$ invariance, changes at
multiple points at early times: particle production in that case would happen
during one particular phase, which would fall within the cases considered
here.

We relaxed the assumption that the dark matter be protected from decay by a
$Z_{2}$ symmetry, and discussed the possible detection of massive, right-
handed neutrinos. We verified that the neutrino is never in thermal
equilibrium in the early universe. We found that the optimal way to search for
the neutrino’s decay products is to employ high-energy neutrino telescopes
such as IceCube, which at present rules out neutrinos with lifetimes in excess
of $10^{29}$ s up to masses as large as $10^{8}$ GeV.

Future tests of this mechanism for the generation of the cosmological dark
matter will depend upon observational indications for the lack of an
inflationary period, as well as in favor of a $CPT$ symmetric early phase; the
detection of the debris from the decay of a superheavy species compatible with
that species being a right-handed neutrino would also be a smoking-gun
signature for this scenario.

###### Acknowledgements.

LM and SP are partly supported by the U.S. Department of Energy grant number
DE-SC0010107.

## References

  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * Boyarsky  _et al._ [2019] A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, Sterile neutrino Dark Matter, Prog. Part. Nucl. Phys. 104, 1 (2019), arXiv:1807.07938 [hep-ph] .
  * Boyle  _et al._ [2018a] L. Boyle, K. Finn, and N. Turok, The Big Bang, CPT, and neutrino dark matter, (2018a), arXiv:1803.08930 [hep-ph] .
  * Boyle  _et al._ [2018b] L. Boyle, K. Finn, and N. Turok, CPT-Symmetric Universe, Phys. Rev. Lett. 121, 251301 (2018b), arXiv:1803.08928 [hep-ph] .
  * Birrell and Davies [1984] N. D. Birrell and P. C. W. Davies,  _Quantum Fields in Curved Space_, Cambridge Monographs on Mathematical Physics (Cambridge Univ. Press, Cambridge, UK, 1984).
  * Mukhanov and Winitzki [2007] V. Mukhanov and S. Winitzki,  _Introduction to quantum effects in gravity_ (Cambridge University Press, 2007).
  * DeWitt [1975] B. S. DeWitt, Quantum Field Theory in Curved Space-Time, Phys. Rept. 19, 295 (1975).
  * Lawrence and Martinec [1996] A. E. Lawrence and E. J. Martinec, String field theory in curved space-time and the resolution of space - like singularities, Class. Quant. Grav. 13, 63 (1996), arXiv:hep-th/9509149 .
  * Kuzmin and Tkachev [1999] V. A. Kuzmin and I. I. Tkachev, Ultrahigh-energy cosmic rays and inflation relics, Phys. Rept. 320, 199 (1999), arXiv:hep-ph/9903542 .
  * Chung  _et al._ [2012] D. J. H. Chung, L. L. Everett, H. Yoo, and P. Zhou, Gravitational Fermion Production in Inflationary Cosmology, Phys. Lett. B 712, 147 (2012), arXiv:1109.2524 [astro-ph.CO] .
  * Allahverdi  _et al._ [2020] R. Allahverdi _et al._ , The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe 10.21105/astro.2006.16182 (2020), arXiv:2006.16182 [astro-ph.CO] .
  * Barenboim and Rasero [2014] G. Barenboim and J. Rasero, Structure Formation during an early period of matter domination, JHEP 04, 138, arXiv:1311.4034 [hep-ph] .
  * Bassett  _et al._ [2006] B. A. Bassett, S. Tsujikawa, and D. Wands, Inflation dynamics and reheating, Rev. Mod. Phys. 78, 537 (2006), arXiv:astro-ph/0507632 .
  * Chung  _et al._ [1999] D. J. H. Chung, E. W. Kolb, and A. Riotto, Production of massive particles during reheating, Phys. Rev. D 60, 063504 (1999), arXiv:hep-ph/9809453 .
  * Acharya  _et al._ [2009] B. S. Acharya, G. Kane, S. Watson, and P. Kumar, A Non-thermal WIMP Miracle, Phys. Rev. D 80, 083529 (2009), arXiv:0908.2430 [astro-ph.CO] .
  * Hashiba and Yokoyama [2019] S. Hashiba and J. Yokoyama, Dark matter and baryon-number generation in quintessential inflation via hierarchical right-handed neutrinos, Phys. Lett. B 798, 135024 (2019), arXiv:1905.12423 [hep-ph] .
  * Feynman [1948] R. P. Feynman, Space-time approach to nonrelativistic quantum mechanics, Rev. Mod. Phys. 20, 367 (1948).
  * Atwood  _et al._ [2009] W. B. Atwood _et al._ (Fermi-LAT), The Large Area Telescope on the Fermi Gamma-ray Space Telescope Mission, Astrophys. J. 697, 1071 (2009), arXiv:0902.1089 [astro-ph.IM] .
  * Aartsen  _et al._ [2017] M. Aartsen, M. Ackermann, J. Adams, J. Aguilar, M. Ahlers, M. Ahrens, I. Al Samarai, D. Altmann, K. Andeen, T. Anderson, _et al._ , Constraints on galactic neutrino emission with seven years of icecube data, The Astrophysical Journal 849, 67 (2017).
  * Albert  _et al._ [2020] A. Albert _et al._ (HAWC), Constraints on the Emission of Gamma Rays from M31 with HAWC 10.3847/1538-4357/ab7999 (2020), arXiv:2001.04065 [astro-ph.HE] .
  * Coogan  _et al._ [2020] A. Coogan, L. Morrison, and S. Profumo, Direct Detection of Hawking Radiation from Asteroid-Mass Primordial Black Holes, (2020), arXiv:2010.04797 [astro-ph.CO] .
  * Coogan  _et al._ [2021] A. Coogan, A. Moiseev, L. Morrison, and S. Profumo, Hunting for Dark Matter and New Physics with (a) GECCO, (2021), arXiv:2101.10370 [astro-ph.HE] .
  * Bauer  _et al._ [2020] C. W. Bauer, N. L. Rodd, and B. R. Webber, Dark matter spectra from the electroweak to the planck scale (2020), arXiv:2007.15001 [hep-ph] .


# Updating Street Maps using
Changes Detected in Satellite Imagery

Favyen Bastani1, Songtao He1, Satvat Jagwani1, Mohammad Alizadeh1,
Hari Balakrishnan1, Sanjay Chawla2, Sam Madden1, Mohammad Amin Sadeghi2 1MIT
CSAIL, 2Qatar Computing Research Institute, HBKU
<EMAIL_ADDRESS><EMAIL_ADDRESS>2QA

(2021)

###### Abstract.

Accurately maintaining digital street maps is labor-intensive. To address this
challenge, much work has studied automatically processing geospatial data
sources such as GPS trajectories and satellite images to reduce the cost of
maintaining digital maps. An end-to-end _map update_ system would first
process geospatial data sources to extract insights, and second leverage those
insights to update and improve the map. However, prior work largely focuses on
the first step of this pipeline: these _map extraction_ methods infer road
networks from scratch given geospatial data sources (in effect creating
entirely new maps), but do not address the second step of leveraging this
extracted information to update the existing digital map data. In this paper,
we first explain why current map extraction techniques yield low accuracy when
extended to update existing maps. We then propose a novel method that
leverages the progression of satellite imagery over time to substantially
improve accuracy. Our approach first compares satellite images captured at
different times to identify portions of the physical road network that have
visibly changed, and then updates the existing map accordingly. We show that
our change-based approach reduces map update error rates four-fold.

automatic map update, machine learning

††journalyear: 2021††copyright: rightsretained††conference: 29th International
Conference on Advances in Geographic Information Systems; November 2–5, 2021;
Beijing, China††booktitle: 29th International Conference on Advances in
Geographic Information Systems (SIGSPATIAL ’21), November 2–5, 2021, Beijing,
China††doi: 10.1145/3474717.3483651††isbn: 978-1-4503-8664-7/21/11††ccs:
Applied computing Cartography

##  1\. Introduction

Maintaining street maps is a labor-intensive process. As a result, many
techniques have been proposed to automate parts of this process by using
geospatial data sources. Current map extraction techniques (Alshehhi et al.,
2017; Bastani et al., 2018b; Cheng et al., 2017; Costea et al., 2017; Máttyus
et al., 2017; Panboonyuen et al., 2017; Vakalopoulou et al., 2015) primarily
rely on satellite imagery due to its global availability, while some
techniques use GPS trajectories.

A key problem with current techniques is that they are designed to infer road
networks from scratch — however, given that we already have existing high
quality maps that cover the vast majority of the world, these inferred road
networks are not directly useful. Instead, an end-to-end map update system
must process geospatial data sources to update and improve existing digital
maps.

###  1.1. Map Extraction Methods Perform Poorly on Map Update

We first consider extending current map extraction techniques for updating
maps. We will show that these methods perform poorly on map update, creating
many false positive updates.

Suppose that the current live digital map has a set of roads $R$. We begin by
applying a map extraction method to process the most recent satellite imagery
(spanning the world). This method produces another set of roads $T$ detected
in the satellite imagery. Simply replacing $R$ with $T$ would not be sensible
for several reasons:

  1. (1)

For roads that appear in both $R$ and $T$, given that $R$ is largely human-
curated, it captures those roads substantially more accurately than $T$.

  2. (2)

$R$ includes roads such as tunnels that cannot be detected by the map
extraction method.

  3. (3)

Roads in $R$ are labeled with rich annotations such as street names, speed
limits, etc. that would be lost in the replacement.

We could instead try to combine $R$ and $T$: if for a road segment $s$,
$s\notin R$ and $s\in T$, we add $s$ to the map. (We could also remove
segments $s$ where $s\in R$ and $s\notin T$, but this would prune roads that
are not visible in the satellite image such as tunnels and roads occluded by
buildings or trees, so we do not consider it further.)

Figure 1. False positive detections when using MAiD to update OpenStreetMap.
We show new 2017 imagery, in which the model erroneously detects roads, on the
right. We include older 2015 imagery on the left to show how change detection
can reduce false positives — all of these false positives could have been
eliminated by comparing satellite images over time, and determining that these
are not recent road network changes.

However, in practice, this approach yields a large number of false positives,
where the map extraction method erroneously outputs many road segments in
places where there are no roads. We demonstrate this issue by using a state-
of-the-art inference method, MAiD (Bastani et al., 2018a), to update
OpenStreetMap (Haklay and Weber, 2008). We select a region of Massachusetts
where OpenStreetMap has high coverage, and manually remove 204 groups of roads
from the map that correspond to new construction between 2015 and 2017. We
apply MAiD to recover these removed roads, and score its performance in terms
of precision and recall comparing recovered groups of roads to the manually
removed groups. At 80% recall, MAiD yields only 67% precision — far too low
for full automation to be a realistic option. We show examples of incorrect
detections in Figure 1.

Many of these detections, such as the Runway, Walkway, and Crop Field examples
in Figure 1, arise due to paths that are virtually indistinguishable from
roads: it is restrictions set by policy governing the use of those paths, and
not physical characteristics of the paths, that make them unsuitable for
traversal by motor vehicles. Thus, simply improving the machine learning
techniques and models used in map extraction methods is unlikely to improve
accuracy; instead, a fundamental shift in approach is required.

###  1.2. Key Tasks in Maintaining Maps

To determine how accuracy can be improved, we first take a step back and
identify the major challenges associated with maintaining maps. In particular,
digital maps have near-complete coverage in most parts of the world: a 2015
study111See https://blog.mapbox.com/how-complete-is-
openstreetmap-7c369787af6e. found that, in 154 of 233 considered countries and
territories, the length of roads in OpenStreetMap exceeded the estimated road
network length reported in the CIA World Factbook, implying that the map
provided excellent coverage in those countries. Then, given that digital maps
have good coverage, identifying pre-existing roads constructed several years
ago is not a key issue: in almost all cases, such roads already appear in the
map.

Instead, _the key challenge in maintaining maps is keeping the map up-to-date
with changes in the physical road network_. Indeed, in the US alone, an
estimated 30K km of roads are constructed each year222“Public Road Mielage”,
FHWA, https://www.fhwa.dot.gov/policyinformation/statistics/2013/vmt422c.cfm.,
and map vendors spend hundreds of millions of dollars annually to keep maps
up-to-date.

This presents the question: Can we develop techniques that directly tackle the
key challenge of identifying physical road network changes, to substantially
improve accuracy at maintaining maps over map extraction methods?

###  1.3. Map Update through Change Detection

To substantially improve accuracy, we propose leveraging a source of data
largely overlooked in prior work: the progression of satellite imagery over
time. By comparing satellite images captured at different times, we can hone
in on portions of the road network that visibly changed over the satellite
image time series. Focusing on segments that visibly changed over time enables
us to disambiguate false positive road segments from genuinely new
constructions: for example, all of the false positive detections in Figure 1
arose from curvilinear features such as walkways and crop field paths that did
not undergo any recent change; thus, by processing the satellite image time
series, we can determine that these are pre-existing features, and should not
be added to the map.

To implement the proposed solution, we must detect changed road segments
across satellite images. However, most existing change detection methods are
fully supervised. They rely on collecting annotated pairs of images where
change has occurred. Since newly constructed roads are rare relative to the
size of the map, collecting positive examples of new roads for such a dataset
is tedious and costly. Additionally, the diversity in visual appearance of
roads makes change detection especially challenging. Furthermore, in Section
4, we show that prior work in unsupervised change detection exhibits low
accuracy when applied for detecting new construction.

Instead, we develop a two-stage approach that requires no hand-labeling for
comparing satellite images over time to detect new roads. In the first stage,
we apply a novel change-seeking iterative tracing procedure to detect recently
constructed roads that are missing from the existing map. Our method uses
ground-truth road labels derived from the existing map dataset to avoid
needing new annotations, and detects new roads that appear in an up-to-date
satellite image but are not visible in an old image.

Though the first stage is effective at detecting new construction, it
nevertheless yields false positive detections when occlusion and other factors
yield visual differences between the old and up-to-date images despite no
actual change. Thus, in the second stage, we propose a novel self-supervised
change detection approach to further improve precision. We train a CNN to
classify whether windows of two aligned satellite images captured at different
times are cropped at the same window (matching) or at different windows
(mismatched). The CNN learns to match features like road markers to determine
whether two crops are matching or mismatched. To apply the model for
inference, we provide it with matching crop pairs around roads detected
through tracing, and we only retain detections that fool the model into
classifying the pair as mismatched, suggesting the presence of new
construction.

We evaluate our approach on a large-scale dataset consisting of 4800 km2 of
satellite imagery. We apply approaches to improve an existing map dataset,
OpenStreetMap, by adding newly constructed roads to the map. At 50% recall,
our approach reduces error rates over existing state-of-the-art map inference
methods four-fold, from 12% to 3%. Our code and data is available at
https://favyen.com/mapupdate.

In summary, our contributions are:

Figure 2. Our two-stage approach for detecting newly constructed roads in
satellite imagery. Here, a pair of aligned images reflecting a new
neighborhood is correctly detected.

  * •

We propose a novel approach for updating street maps that leverages the
progression of satellite imagery over time. In contrast to prior work that
focuses on inferring all roads in a satellite image, our approach tackles the
directly practical map update problem.

  * •

We develop a two-stage approach for detecting new construction in satellite
imagery. Our approach does not require any hand-labeling. In the second stage,
we propose a novel self-supervised road-masking approach, where we train a CNN
classifier in a self-supervised manner to detect change.

  * •

We evaluate our approach on a large-scale dataset consisting of 3000 km2 of
satellite imagery in the Boston area for training and 1800 km2 in northeastern
Massachusetts for testing. We apply approaches to improve an existing map
dataset, OpenStreetMap, by incorporating newly constructed roads into the
dataset. At 50% recall, our approach reduces error rates over existing state-
of-the-art map inference methods four-fold, from 12% to 3%.

##  2\. Related Work

Road Extraction. Automatically inferring roads from satellite imagery is a
well-studied problem. Recent road extraction methods generally apply
convolutional neural networks (CNNs) to segment imagery for roads (Batra et
al., 2019; Costea et al., 2017; Mosinska et al., 2019; Panboonyuen et al.,
2017; Yang et al., 2019; Zhou et al., 2018), and apply various methods to
post-process the segmentation output and derive vector road network graphs.
Cheng et al. apply binary thresholding, morphological thinning, and line
following to extract a road network from segmentation probabilities (Cheng et
al., 2017). DeepRoadMapper (Máttyus et al., 2017) proposes several additional
heuristic and learning-based refinement steps, including removing short edges
and identifying potential missed roads. Some methods propose alternatives to
treating road extraction as an image segmentation problem. Alshehhi et al.
build the road network with a region adjacent graph that forms narrow
elongated regions along roads (Alshehhi et al., 2017). RoadTracer (Bastani et
al., 2018b) and PolyMapper (Li et al., 2019) propose an iterative tracing
framework to extract road networks: they train a CNN to output the
directionality of roads at each pixel, and employ an iterative search guided
by the CNN to trace the road network. VecRoad extends the iterative tracing
approach with a flexible step size and joint learning tasks (Tan et al.,
2020), and Neural Turtle Graphics extends it with a sequential generative
model (Chu et al., 2019). Another recent technique, Sat2Graph, proposes a one-
shot road extraction process where a CNN directly predicts the positions of
road network vertices and edges (He et al., 2020).

However, broadly, these approaches are unable to reason about false positive
detections made by the CNN such as those in Figure 1, especially for paths
that appear visually similar to roads in the satellite image but are not
suitable for traversal by motor vehicles. As a result, when road extraction
methods are applied for the practical task of updating existing maps, they
incorporate many non-road paths into the map, thereby substantially
deteriorating the quality of the map data. In contrast, in most of the world
where existing maps have good coverage, our method accurately keeps maps up to
date with new construction by comparing satellite imagery over time, and only
updating the map in areas where change is detected across images.

GPS Trajectories for Updating Maps. Inferring roads from GPS trajectory data
has also been studied (Ahmed et al., 2015; Biagioni and Eriksson, 2012; Cao
and Krumm, 2009; Davies et al., 2006; Prabowo et al., 2019; He et al., 2018;
Stanojevic et al., 2018). Two works in this space, CrowdAtlas (Wang et al.,
2013) and COBWEB (Shan et al., 2015), propose _map update_ methods to
incorporate new roads into existing maps. However, since these methods do not
consider GPS time series data (comparing older trajectories to recent
trajectories), they exhibit false positive errors similar to satellite image
road extraction methods due to GPS noise. Additionally, due to the lack of a
ground truth test set, prior work have not incorporated a quantitative
evaluation of the map update portion of those methods, and instead
qualitatively show results at detecting a small number of roads missing from
an existing map. While it may be possible to accurately update existing maps
by comparing old and recent GPS trajectory data, this has not been studied in
prior work; in our approach, we focus on using satellite image time series
data, since satellite imagery is globally available.

Change Detection. Change detection in satellite imagery has previously been
studied for detecting damage from natural disasters and armed conflict.
Gueguen et al. employ a semi-supervised learning approach to identify damaged
regions by comparing images before and after a calamity(Gueguen and Hamid,
2015). However, adapting supervised and semi-supervised change detection
methods for maintaining maps is difficult: annotating examples of new roads is
highly time-consuming because the density of new construction is low.
Unsupervised change detection methods have also been proposed (de Jong and
Bosman, 2019), but we will show in Section 4 that these methods exhibit poor
accuracy when used to identify newly constructed roads.

Figure 3. Illustration of change-seeking iterative tracing. The blue circle
visualizes tracing confidences at the center pixel, with lighter colors
indicating higher confidence. The existing map is green, and traced roads are
yellow. On each step, we compare confidences in $M_{\text{new}}$ (top) against
confidences in $M_{\text{old}}$ (bottom). On the right, this comparison helps
to avoid tracing a pre-existing pedestrian path (red circle) when the model
outputs similar confidences in both images.

##  3\. Detecting Street Network Changes

Existing state-of-the-art map inference methods are designed to infer maps
from scratch rather than update existing maps. When applied for updating maps
with new roads, these methods show poor accuracy — the number of false
positives overwhelms the few cases of new construction. Figure 1 shows
examples of false positive detections that arise when we apply these methods
to add missing roads to OpenStreetMap. Oftentimes, these false positives arise
due to pre-existing paths such as air fields and cycling paths that appear
visually similar to roads (Figure 1). We propose a novel approach that, in
contrast to prior work, directly detects newly constructed roads by comparing
an up-to-date satellite image against an old image. Figure 2 summarizes our
approach.

Let $M_{\text{old}}$ and $M_{\text{new}}$ be old and up-to-date satellite
images of the same region, and let $G$ be the road network graph of that
region in the existing map. Each vertex $v\in G$ is annotated with a pixel
$(i,j)$ corresponding to its location in the images, and edges correspond to
roads. At a high level, in our approach, we detect new roads by identifying
roads that appear in $M_{\text{new}}$ but not in $M_{\text{old}}$ or $G$. In
most of the world where existing maps have good coverage, roads detected in
both $M_{\text{new}}$ and $M_{\text{old}}$ likely are not roads at all, but
instead non-road paths (e.g., the examples in Figure 1); by comparing
$M_{\text{new}}$ and $M_{\text{old}}$, our approach avoids these false
positives.

We first apply a novel change-seeking iterative tracing procedure that adapts
MAiD (Bastani et al., 2018a) to selectively trace roads in $M_{\text{new}}$
that appear in neither $M_{\text{old}}$ nor the existing map $G$, i.e., roads
that were constructed after $M_{\text{old}}$ was captured. Our method traces
roads along segments where a CNN model has high confidence in $M_{\text{new}}$
and low confidence in $M_{\text{old}}$. Although this procedure improves
precision over prior work, we find that it produces false positives when
differences in off-nadir angle and lighting or visible non-construction
activity result in a sharp increase in the CNN confidence from
$M_{\text{old}}$ to $M_{\text{new}}$ despite no new roads.

Thus, we propose a novel self-supervised change detection method to
automatically prune these remaining false positives in the second stage of our
approach. Our method selectively identifies road network changes so that
visible non-construction activity does not result in a false positive. The
final road detections have high precision, and can be used to improve real-
world maps through automatic merging or human validation.

In Section 3.1, we describe our method to obtain an initial set of candidate
roads using change-seeking iterative tracing. We then introduce our novel
self-supervised road-masking approach in Section 3.2.

###  3.1. Change-Seeking Iterative Tracing

In the first stage, we apply a MAiD (Bastani et al., 2018a) model to segment
the images for tracing confidences; each tracing confidence indicates the
likelihood that a road passes a pixel in a particular direction (angle). In
MAiD, these tracing confidences are used by an iterative tracing algorithm to
draw roads along directions with high confidence. In contrast, we develop a
change-seeking iterative tracing process that avoids many of the false
positives in Figure 1 that involve pre-existing non-road paths by comparing
confidences extracted from $M_{\text{new}}$ and from $M_{\text{old}}$ to draw
roads only along directions with substantially higher confidence in
$M_{\text{new}}$, which suggests the presence of a new road. Below, we detail
each of the components in our approach.

Model. We use the MAiD CNN model architecture from (Bastani et al., 2018a).
Given an image $M$, the model produces a three-dimensional matrix $P$, where
$P_{i,j,k}$ is the probability that a road passes the pixel at $(i,j)$ in a
direction specified by $k$. $P$ includes 64 channels, and the $k$th channel
indicates the likelihood that there is a road at an angle between
$k\frac{2\pi}{64}$ and $(k+1)\frac{2\pi}{64}$ from a pixel. At a pixel that
falls along a straight road, channels specifying opposite directions along the
road would both have high confidence. $P$ is output at one-fourth the input
resolution.

Training. During training, we construct an example input by first randomly
deciding whether to input $M_{\text{old}}$ or $M_{\text{new}}$, and then
selecting a random 2D window in the image. We create training labels using the
existing road network graph $G$. The input is a $256\times 256\times 3$ image
(where the three channels are derived from RGB satellite imagery), and the
label is a $64\times 64\times 64$ matrix of tracing confidences. We train the
model using binary cross entropy loss, averaged across pixels and the 64
channels.

Tracing Procedure. We develop a novel change-seeking iterative tracing
procedure that adapts the tracing process used in prior work (Bastani et al.,
2018b, a) to focus tracing on newly constructed roads and avoid inferring
false positives along pre-existing non-road paths. Figure 3 illustrates the
tracing process.

Iterative tracing starts from an initial pixel known to lie on the road
network, and follows directions with high confidence in $P$ at the current
pixel in a depth-first search (DFS) process. The inferred road network
consists of the paths that were followed during the search.

We first compute $P_{\text{old}}$ from $M_{\text{old}}$ and $P_{\text{new}}$
from $M_{\text{new}}$. We use the existing map $G$ as a base map: our tracing
procedure starts from pixels in the base map, and follows directions with high
confidence in $P_{\text{new}}$ and low confidence in $P_{\text{old}}$ in a
depth-first search (DFS) process. The inferred road network consists of paths
that were followed during the search. Let $G^{\prime}$ be the current road
network state, which we extend during tracing. We initialize $G^{\prime}$ by
densifying $G$ so that vertices are at most 10 m apart. We use each vertex in
$G^{\prime}$ as a starting pixel for tracing, and append all vertices to a DFS
stack.

On each tracing iteration, we consider the pixel $(i,j)$ at the top of the
search stack. We identify the highest confidence direction $k$ in
$P_{\text{new}}[i,j]$ that has an angular distance of at least $30^{\circ}$
from any existing edges in $G^{\prime}$ at $(i,j)$. This prevents re-tracing
roads that are already covered by $G^{\prime}$. Additionally, though, we only
trace from $(i,j)$ if $P_{\text{new}}[i,j,k]\geq T_{\text{new}}$ and
$P_{\text{old}}[i,j,k]<T_{\text{old}}$, i.e., if the tracing confidence in the
up-to-date imagery exceeds a threshold while the confidence for the same pixel
and direction in the old imagery is small. Thus, our tracing procedure only
follows new roads that are not reflected in $P_{\text{old}}$. If we decide to
trace, then we add a new vertex $v=(i+\cos\alpha,j+\sin\alpha)$ to
$G^{\prime}$, where $\alpha$ is the angle corresponding to $k$, and add an
edge from $(i,j)$ to $v$. We then push $v$ onto the DFS stack. Otherwise, if
we decide not to trace, we pop $(i,j)$ from the stack.

We terminate tracing once the DFS stack is empty. At this point, each
connected component in $G^{\prime}-G$ is a candidate group of roads. By
comparing $P_{\text{old}}$ and $P_{\text{new}}$ during the tracing procedure,
we are able to avoid tracing along pre-existing non-road paths.

###  3.2. Self-Supervised Learning for Selective Change Detection

Figure 4. False positives after change-seeking iterative tracing. We show road
detections in yellow, which correspond to a railway and three pedestrian
paths. Lighting and off-nadir angle differences and visible non-construction
activity lead to low tracing confidences in the old image (top) and high
confidence in the new image (bottom), which results in the false positive road
detections.

Change-seeking iterative tracing improves precision over prior work, but still
produces false positives when $P_{\text{new}}$ reflects higher confidence than
$P_{\text{old}}$ due to angle and lighting differences or visible activity
without new construction between $M_{\text{old}}$ and $M_{\text{new}}$. Figure
4 shows example false positives. On the left, a railway and a pedestrian path
are partially occluded by shadows in the old image, but are visible in the up-
to-date image. On the right, although there is visible activity (bridge
demolition and crosswalk painting), there are no new roads, and two walkways
are incorrectly detected.

In the second stage of our approach, we apply change detection to filter the
candidates generated in the first stage by pruning these false positives.
However, as we will show in the evaluation, unsupervised change detection
methods are unable to robustly distinguish false positives due to angle,
lighting, and other aforementioned differences. Supervised methods are also
impractical: paired examples of new construction are tedious to annotate due
to their low density. Instead, we develop a novel approach that applies self-
supervised learning to selectively identify road network changes.

Figure 5. Examples of matching pairs (left) and mismatched pairs (right) that
we generate during training in the second stage. (a) shows the windows
selected for each example pair over a large region, (b) shows the satellite
images cropped at those windows, and (c) shows the images after masking.

In our self-supervised learning procedure, we train a classifier that inputs a
pair of windows of old and up-to-date images. The input may either be a
matching pair, where $M_{\text{old}}$ and $M_{\text{new}}$ are cropped at the
same window, or a mismatched pair, where $M_{\text{old}}$ and $M_{\text{new}}$
are cropped at disjoint windows. We train the classifier to distinguish
matching pairs from mismatched pairs. We generate training examples by
deciding to create a matching pair or mismatched pair with equal probability.
To generate a matching pair, we randomly pick one window and crop both images
at that window. To generate a mismatched pair, we randomly pick two disjoint
windows. The classifier learns to match features between the images to
determine whether they are taken at the same window despite differences in
matching pairs such as shadows, camera angle, and non-construction activity
that make unsupervised change detection methods ineffective on this task. We
show example training pairs in Figure 5.

To apply the filter for inference, we execute the classifier on crops of
$M_{\text{old}}$ and $M_{\text{new}}$ taken at a window around connected
components of roads detected during tracing. Although this pair is matching,
where the crops are aligned, if there is substantial change in the images due
to newly constructed roads, the crops are likely to fool the classifier into
outputting a higher probability for the “mismatched” class. Thus, we prune
candidate roads if the “mismatched” probability falls below a threshold.

However, in practice, new construction is often adjacent to pre-existing
roads, buildings, and other structures. If the classifier observes the same
structures in both the old and up-to-date windows, it can determine that the
pair is matching with high confidence. Thus, we develop a masking approach
that focuses the classifier on the detected candidate road.

Figure 6. CNN model architecture for the classifier used in our self-
supervised road-masking approach. Figure 7. Classifier outputs on example
groups of candidate roads inferred through change-seeking iterative tracing.
At the bottom, we show the matching class probability.

Model. The model input consists of 7 channels: 3 channels from the crop
$I_{\text{old}}$ of $M_{\text{old}}$, 3 channels from the crop
$I_{\text{new}}$ of $M_{\text{new}}$, and 1 channel containing a mask $S$.
$S[i,j]$ is either 0 or 1, and if $S[i,j]=0$, then we zero the corresponding
values in the crop channels, i.e.,
$I_{\text{old}}[i,j]=I_{\text{new}}[i,j]=0$. Because the size of the input
during inference varies based on the candidate road, we use a fully
convolutional CNN architecture consisting of 6 encoder layers followed by 5
decoder layers. Figure 6 shows the model architecture. The model outputs a
probability at each pixel that the input example is “matching”. We train the
CNN with cross entropy loss, averaged over only pixels where $S[i,j]=1$.

Training. On each training step, we construct a matching example with 50%
probability and a mismatched example with 50% probability. In both cases, we
begin by computing the mask $S$ that we will apply to the imagery crop inputs.
During inference, $S[i,j]$ will be 1 only near a group of candidate roads
obtained through tracing. For effective training, $S$ must be similar to what
we will provide during inference, e.g., $S$ should predominantly cover roads.
Thus, we leverage $G$ to compute the mask: we randomly select a vertex $v_{0}$
in $G$, and perform a breadth-first-search from $v_{0}$ to derive a subgraph
$H$ that will determine $S$. During the search, we add each traversed edge to
$H$, and terminate the search once the length of the bounding box containing
$H$ exceeds a threshold $T_{\text{box}}$. We vary $T_{\text{box}}$ to ensure
that training examples have diverse mask sizes, since during inference,
candidate groups of roads may exhibit different sizes; specifically, we pick
$T_{\text{box}}$ uniformly between 50 m and 150 m. We set $S[i,j]=1$ if pixel
$(i,j)$ falls within 20 meters of some edge in $H$, and set $S[i,j]=0$
otherwise.

(a) Precision and recall comparing detections against hand-labeled new
construction.

Method | APLS
---|---
DeepRoadMapper | 0.52
MAiD | 0.49
UnstructChange | 0.48
CmpOnly | 0.51
Cmp+Filter (ours) | 0.57
(b) APLS.

Figure 8. Evaluation of our approach (Cmp+Filter) and the baselines on
detecting newly constructed roads. Figure 9. Example detections of new roads
produced by our Cmp+Filter method, including two incorrect outputs on the
bottom right. We show the basemap in blue and detections in yellow.

To create a mismatched example, we select two windows $W_{\text{old}}$ and
$W_{\text{mask}}$. $W_{\text{mask}}$ is always a $512\times 512$ window
centered at $v_{0}$. We randomly pick $W_{\text{old}}$ so that
$W_{\text{old}}$ and $W_{\text{mask}}$ are disjoint, but the distance from
$W_{\text{old}}$ to $v_{0}$ is at most 2500 m. Choosing nearby but disjoint
windows when creating mismatched examples is crucial as it yields more
challenging examples where the old and up-to-date crops have similar style
(e.g., both in suburban neighborhoods) but different semantic content. We crop
$M_{\text{old}}$ at $W_{\text{old}}$, and $M_{\text{new}}$ and $S$ at
$W_{\text{mask}}$.

To create a matching example, with 80% probability, we simply crop both
$M_{\text{old}}$ and $M_{\text{new}}$ at $W_{\text{mask}}$. However, in some
cases, tracing may output roads over non-road paths; if we only train the
model on matching examples where the pixels where $S[i,j]=1$ fall on roads,
the model may be ineffective on non-road inputs. Thus, with 20% probability,
we crop $M_{\text{old}}$ and $M_{\text{new}}$ at a random window
$W_{\text{rand}}$.

Figure 5 shows two training examples.

Inference. For each candidate group of roads $H$, we first derive a
corresponding mask $S$ similar to the process during training: $S[i,j]=1$ only
if pixel $(i,j)$ falls within 20 meters of a candidate road. We crop
$M_{\text{old}}$, $M_{\text{new}}$, and $S$ at a window corresponding to the
bounding box of $H$ with 20-meter padding. We then compute the average
probability $p$ that the model outputs over pixels where $S[i,j]=1$. Then,
given a filter threshold $T_{\text{filter}}$, if $p<T_{\text{filter}}$, we
prune the candidate. In Figure 7, we show the average probability on several
roads inferred through change-seeking iterative tracing.

##  4\. Evaluation

We evaluate our approach against existing state-of-the-art map inference and
change detection methods on a task involving automatically updating
OpenStreetMap with new roads. We use 60 cm/pixel resolution satellite imagery
from MassGIS from 2015 and 2017 as our old imagery $M_{\text{old}}$ and up-to-
date imagery $M_{\text{new}}$, and the OpenStreetMap dataset as our road
network graph. We select two disjoint sections of this dataset for training
(3000 $\text{km}^{2}$ in the Boston metro area) and for evaluation (1800
$\text{km}^{2}$ in northeastern Massachusetts).

Metrics. For evaluation, we hand-annotated 204 groups of roads that appear in
$M_{\text{new}}$ but not $M_{\text{old}}$. We prune these roads from
OpenStreetMap to derive a road network $G$ corresponding to a map that has not
yet been updated with the new imagery. We compare the methods in terms of the
precision and recall on recovering the pruned roads. Each approach outputs a
set of map update proposals $P=\\{H_{1},\ldots,H_{n}\\}$, where each $H_{i}$
is a connected component of inferred roads. The pruned roads form a set of
ground truth proposals $P^{*}$. We say a proposal $H_{i}\in P$ matches a
ground truth proposal $H_{j}^{*}\in P^{*}$ if the proposal bounding boxes
intersect. Then, precision and recall are defined as:

$\text{precision}=\frac{\text{\\#
matches}}{|P|}\quad\text{recall}=\frac{\text{\\# matches}}{|P^{*}|}$

Under this metric, each approach yields a precision-recall curve over varying
confidence thresholds (e.g., in our approach, varying $T_{\text{filter}}$).
Because $P^{*}$ is not comprehensive, we discard a proposal $H_{i}\in P$ if it
is a correct example of a road but does not match with any ground truth
proposal.

Although the focus of our approach is on improving the precision of inferred
road segments rather than improving the geometrical accuracy of those
segments, we also evaluate the methods on the latter in terms of Average Path
Length Similarity (APLS) (Van Etten et al., 2018).

Baselines. We evaluate our method against two baselines detailed in Section 2
that implement existing state-of-the-art road inference approaches: MAiD
(Bastani et al., 2018a) and DeepRoadMapper (Máttyus et al., 2017). We apply
these methods on $M_{\text{new}}$ to derive proposed roads, and prune
proposals that correspond to roads already mapped in OpenStreetMap by pruning
segments that fall within 40 m of an edge in $G$. We also evaluate against an
unsupervised satellite image change detection method, UnstructChange (de Jong
and Bosman, 2019), which identifies change by comparing feature maps extracted
from old and up-to-date satellite images through a VGG-19 model trained for
segmentation. To apply UnstructChange, we first obtain candidate roads through
MAiD, and then eliminate candidates where the unsupervised method detects no
change.

Finally, a fourth baseline, denoted CmpOnly, applies the first stage of our
approach only (change-seeking iterative tracing). We denote our full approach
Cmp+Filter.

Results. We show precision-recall curves on detecting new roads over varying
confidence thresholds in Figure 8, and qualitative results in Figure 9.
DeepRoadMapper is unable to achieve higher than 91% precision due to false
positives, many of which correspond to pre-existing non-road paths. MAiD
provides higher precision at lower recalls, but still yields only 88%
precision at 50% recall. At 50% recall, CmpOnly improves precision to 94%, and
Cmp+Filter further improves precision to 97%. Thus, our method effectively
prunes false positives that have lighting and angle differences or visible
activity but no new roads. UnstructChange does not improve performance over
MAiD: it outputs false positives due to non-construction changes (such as
angle and lighting differences) between the old and up-to-date images. While
our work focuses on improving the precision of road detections, Figure 8(b)
shows that our method also yields a 5% improvement in APLS, which measures the
geometrical accuracy of inferred roads.

Overall, our full approach provides near-100% precision at reasonable recall
levels. Precision is crucial because automatic integration of detections into
the street map dataset is only practical if errors are rare – otherwise, the
confusion for users caused by introducing errors may outweigh the benefit from
expanded map coverage.

###  4.1. Updating Maps with New Buildings

Street maps contain numerous annotations besides roads. In this section, we
show that a simple adaption of our approach for detecting new construction of
buildings yields high accuracy.

(a) Precision and recall curves.

(b) Example inferred buildings.

Figure 10. Evaluation on detecting newly constructed buildings. At bottom, we
show detections produced by Cmp+Filter.

Baselines. We evaluate against two baselines. BldgSeg implements segmentation-
based building extraction (Li et al., 2018), applying a deep CNN to segment
imagery and then extracting building polygons from the segmentation
probabilities. BldgSeg Cmp applies the first stage of our approach, which we
adapt for buildings below.

First Stage. Our change-seeking iterative tracing method is effective at
tracing road networks, but is not applicable for tracing building polygons.
Instead, we apply the BldgSeg baseline to segment $M_{\text{old}}$ and
$M_{\text{new}}$ for buildings, and derive candidate buildings by comparing
segmentation probabilities. Specifically, we first compute segmentation
probabilities $P_{\text{old}}[i,j]$ and $P_{\text{new}}[i,j]$ using the CNN at
each pixel in the imagery. We then compute a binary image $B_{\text{compare}}$
such that $B_{\text{compare}}=1$ only if $P_{\text{old}}[i,j]<T_{\text{old}}$
and $P_{\text{new}}[i,j]>T_{\text{new}}$. Finally, we extract buildings from
$B_{\text{compare}}$.

Second Stage. We train and apply our self-supervised model as for roads, but
compute the mask $S$ by adding a fixed padding around building polygons
instead of around roads.

Metrics. We evaluate the methods on precision and recall, as defined for
roads. We construct a ground truth set of building polygons by removing 665
buildings from OpenStreetMap that were constructed in northeastern
Massachusetts between 2015 and 2017.

Results. We show results in Figure 10. On buildings, at 30% recall, BldgSeg
yields only 85% precision. BldgSeg Cmp improves precision to 87%, and our
approach, Cmp+Filter, yields 92% precision. This corresponds to an almost two-
fold reduction in error rate ($1-\text{precision}$), from 15% to 8%.

###  4.2. Removing Roads

Figure 11. Applying Cmp+Filter to identify removed roads. We show old
satellite imagery (left) and recent satellite imagery (right), with
OpenStreetMap in blue and the detection in red. Cmp+Filter succeeds in
identifying a shifted highway (top) and bulldozed road (middle), but
incorrectly removes a road occluded by trees (bottom).

Cmp+Filter can be applied in reverse to identify removed roads, where we
identify portions of the existing map that appear in $M_{\text{old}}$ but not
in $M_{\text{new}}$. Although we are unable to identify enough examples of
removed roads to conduct a quantitative evaluation, we show three detections
of removed roads in Figure 11. Cmp+Filter succeeds in identifying a shifted
highway and bulldozed road.

##  5\. Conclusion

Maintaining street maps today is labor-intensive and costly. We find that
existing state-of-the-art street map inference systems exhibit low precision
when applied to update an existing map dataset, OpenStreetMap. By leveraging
multiple satellite images collected at different times, our two-stage approach
complements prior work by identifying roads and buildings that were newly
constructed in the most recent image. Our evaluation on 4800 $km^{2}$ of
satellite imagery shows that our approach is able to update existing maps to
capture new construction with high precision.

## References

  * (1)
  * Ahmed et al. (2015) M. Ahmed, S. Karagiorgou, D. Pfoser, and C. Wenk. 2015\.  A Comparison and Evaluation of Map Construction Algorithms using Vehicle Tracking Data.  _GeoInformatica_ 19, 3 (2015), 601–632.
  * Alshehhi et al. (2017) Rasha Alshehhi, Prashanth Reddy Marpu, Wei Lee Woon, and Mauro Dalla Mura. 2017.  Simultaneous Extraction of Roads and Buildings in Remote Sensing Imagery with Convolutional Neural Networks.  _ISPRS Journal of Photogrammetry and Remote Sensing_ 130 (2017), 139–149.
  * Bastani et al. (2018a) Favyen Bastani, Songtao He, Sofiane Abbar, Mohammad Alizadeh, Hari Balakrishnan, Sanjay Chawla, and Sam Madden. 2018a.  Machine-Assisted Map Editing. In _ACM International Conference on Advances in Geographic Information Systems (SIGSPATIAL)_. 23–32.
  * Bastani et al. (2018b) Favyen Bastani, Songtao He, Sofiane Abbar, Mohammad Alizadeh, Hari Balakrishnan, Sanjay Chawla, Sam Madden, and David DeWitt. 2018b.  RoadTracer: Automatic Extraction of Road Networks from Aerial Images. In _CVPR_.
  * Batra et al. (2019) Anil Batra, Suriya Singh, Guan Pang, Saikat Basu, CV Jawahar, and Manohar Paluri. 2019\.  Improved Road Connectivity by Joint Learning of Orientation and Segmentation. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. 10385–10393.
  * Biagioni and Eriksson (2012) James Biagioni and Jakob Eriksson. 2012.  Map Inference in the Face of Noise and Disparity. In _ACM SIGSPATIAL_.
  * Cao and Krumm (2009) Lili Cao and John Krumm. 2009.  From GPS Traces to a Routable Road Map. In _ACM SIGSPATIAL_. 3–12.
  * Cheng et al. (2017) Guangliang Cheng, Ying Wang, Shibiao Xu, Hongzhen Wang, Shiming Xiang, and Chunhong Pan. 2017\.  Automatic Road Detection and Centerline Extraction via Cascaded End-to-End Convolutional Neural Network.  _IEEE Transactions on Geoscience and Remote Sensing_ 55, 6 (2017), 3322–3337.
  * Chu et al. (2019) Hang Chu, Daiqing Li, David Acuna, Amlan Kar, Maria Shugrina, Xinkai Wei, Ming-Yu Liu, Antonio Torralba, and Sanja Fidler. 2019.  Neural Turtle Graphics for Modeling City Road Layouts. In _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_. 4522–4530.
  * Costea et al. (2017) Dragos Costea, Alina Marcu, Emil Slusanschi, and Marius Leordeanu. 2017.  Creating Roadmaps in Aerial Images With Generative Adversarial Networks and Smoothing-Based Optimization. In _ICCV_.
  * Davies et al. (2006) Jonathan J Davies, Alastair R Beresford, and Andy Hopper. 2006\.  Scalable, Distributed, Real-time Map Generation.  _IEEE Pervasive Computing_ 5, 4 (2006).
  * de Jong and Bosman (2019) Kevin Louis de Jong and Anna Sergeevna Bosman. 2019.  Unsupervised Change Detection in Satellite Images Using Convolutional Neural Networks. In _International Joint Conference on Neural Networks (IJCNN)_. IEEE, 1–8.
  * Gueguen and Hamid (2015) Lionel Gueguen and Raffay Hamid. 2015.  Large-Scale Damage Detection Using Satellite Imagery. In _CVPR_.
  * Haklay and Weber (2008) Mordechai Haklay and Patrick Weber. 2008.  OpenStreetMap: User-Generated Street Maps.  _IEEE Pervasive Computing_ 7, 4 (2008), 12–18.
  * He et al. (2018) Songtao He, Favyen Bastani, Sofiane Abbar, Mohammad Alizadeh, Hari Balakrishnan, Sanjay Chawla, and Sam Madden. 2018\.  RoadRunner: Improving the Precision of Road Network Inference from GPS Trajectories. In _ACM SIGSPATIAL_.
  * He et al. (2020) Songtao He, Favyen Bastani, Satvat Jagwani, Mohammad Alizadeh, Hari Balakrishnan, Sanjay Chawla, Mohamed M Elshrif, Samuel Madden, and Amin Sadeghi. 2020\.  Sat2Graph: Road Graph Extraction through Graph-Tensor Encoding.  _Proceedings of the European Conference on Computer Vision (ECCV)_ (2020).
  * Li et al. (2018) Weijia Li, Conghui He, Jiarui Fang, and Haohuan Fu. 2018\.  Semantic Segmentation Based Building Extraction Method Using Multi-Source GIS Map Datasets and Satellite Imagery. In _CVPR_.
  * Li et al. (2019) Zuoyue Li, Jan Dirk Wegner, and Aurélien Lucchi. 2019\.  Topological Map Extraction from Overhead Images. In _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_. 1715–1724.
  * Máttyus et al. (2017) Gellért Máttyus, Wenjie Luo, and Raquel Urtasun. 2017\.  DeepRoadMapper: Extracting Road Topology From Aerial Images. In _CVPR_.
  * Mosinska et al. (2019) Agata Mosinska, Mateusz Koziński, and Pascal Fua. 2019\.  Joint Segmentation and Path Classification of Curvilinear Structures.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 42, 6 (2019), 1515–1521.
  * Panboonyuen et al. (2017) Teerapong Panboonyuen, Kulsawasd Jitkajornwanich, Siam Lawawirojwong, Panu Srestasathiern, and Peerapon Vateekul. 2017.  Road Segmentation of Remotely-Sensed Images using Deep Convolutional Neural Networks with Landscape Metrics and Conditional Random Fields.  _Remote Sensing_ 9, 7 (2017), 680.
  * Prabowo et al. (2019) Arian Prabowo, Piotr Koniusz, Wei Shao, and Flora D Salim. 2019\.  COLTRANE: Convolutional Trajectory Network for Deep Map Inference. In _Proceedings of the ACM International Conference on Systems for Energy-Efficient Buildings, Cities, and Transportation_. 21–30.
  * Shan et al. (2015) Zhangqing Shan, Hao Wu, Weiwei Sun, and Baihua Zheng. 2015\.  COBWEB: A Robust Map Update System using GPS Trajectories. In _Proceedings of the ACM International Joint Conference on Pervasive and Ubiquitous Computing_. 927–937.
  * Stanojevic et al. (2018) Rade Stanojevic, Sofiane Abbar, Saravanan Thirumuruganathan, Sanjay Chawla, Fethi Filali, and Ahid Aleimat. 2018.  Robust Road Map Inference through Network Alignment of Trajectories. In _Proceedings of the 2018 SIAM International Conference on Data Mining_. SIAM.
  * Tan et al. (2020) Yong-Qiang Tan, Shang-Hua Gao, Xuan-Yi Li, Ming-Ming Cheng, and Bo Ren. 2020.  VecRoad: Point-based Iterative Graph Exploration for Road Graphs Extraction. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. 8910–8918.
  * Vakalopoulou et al. (2015) Maria Vakalopoulou, Konstantinos Karantzalos, Nikos Komodakis, and Nikos Paragios. 2015\.  Building Detection in Very High Resolution Multispectral Data with Deep Learning Features. In _IGARSS_.
  * Van Etten et al. (2018) Adam Van Etten, Dave Lindenbaum, and Todd M Bacastow. 2018\.  SpaceNet: A Remote Sensing Dataset and Challenge Series.  _arXiv preprint arXiv:1807.01232_ (2018).
  * Wang et al. (2013) Yin Wang, Xuemei Liu, Hong Wei, George Forman, Chao Chen, and Yanmin Zhu. 2013\.  CrowdAtlas: Self-Updating Maps for Cloud and Personal Use. In _Proceeding of the International Conference on Mobile Systems, Applications, and Services (MobiSys)_. 27–40.
  * Yang et al. (2019) Xiaofei Yang, Xutao Li, Yunming Ye, Raymond YK Lau, Xiaofeng Zhang, and Xiaohui Huang. 2019\.  Road Detection and Centerline Extraction via Deep Recurrent Convolutional Neural Network U-Net.  _IEEE Transactions on Geoscience and Remote Sensing_ 57, 9 (2019), 7209–7220.
  * Zhou et al. (2018) Lichen Zhou, Chuang Zhang, and Ming Wu. 2018.  D-LinkNet: LinkNet with Pretrained Encoder and Dilated Convolution for High Resolution Satellite Imagery Road Extraction. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_. 182–186.