# Decaying Dark Matter: Simulations and Weak-Lensing Forecast

Jonathan Hubert  Aurel Schneider  Doug Potter  and Joachim Stadel

###### Abstract

Despite evidence for the existence of dark matter (DM) from very high and low
redshifts, a moderate amount of DM particle decay remains a valid possibility.
This includes both models with very long-lived yet unstable particles or mixed
scenarios where only a small fraction of dark matter is allowed to decay. In
this paper, we investigate how DM particles decaying into radiation affect
non-linear structure formation. We look at the power spectrum and its redshift
evolution, varying both the decay lifetime ($\tau$) and the fraction of
decaying-to-total dark matter ($f$), and we propose a fitting function that
reaches sub-percent precision below $k\sim 10$ h/Mpc. Based on this fit, we
perform a forecast analysis for a Euclid-like weak lensing (WL) survey,
including both massive neutrino and baryonic feedback parameters. We find that
with WL observations alone, it is possible to rule out decay lifetimes smaller
than $\tau=75$ Gyr (at 95 percent CL) for the case that all DM is unstable.
This constraint improves to $\tau=182$ Gyr if the WL data is combined with CMB
priors from the Planck satellite and to $\tau=275$ Gyr if we further assume
baryonic feedback to be fully constrained by upcoming Sunyaev-Zeldovich or
X-ray data. The latter shows a factor of 3.2 improvement compared to
constraints from CMB data alone. Regarding the scenario of a strongly decaying
sub-component of dark matter with $\tau\sim 30$ Gyr or lower, it will be
possible to rule out a decaying-to-total fraction of $f>0.49$, $f>0.21$, and
$f>0.13$ (at the 95 percent CL) for the same three scenarios. We conclude that
the upcoming stage-IV WL surveys will allow us to significantly improve
current constraints on the stability of the dark matter sector.

##  1 Introduction

Current cosmological and astrophysical observations indicate that dark matter
(DM) is present during both the very early and the late universe. Indeed, a
substantial fraction of DM is required to explain the cosmic microwave
background (CMB) radiation originating at the epoch of recombination [1]. At
the same time, weak-lensing shear measurements [2, 3], galaxy clustering
statistics [4], individual galaxy rotation curves, or strong-lensing
observations [5] point towards the need of a comparable DM abundance at low
redshift. Combining all observational evidence therefore leads to the obvious
conclusion that a hypothetical dark matter particle needs to be either stable
or, at least, very long-lived.

Current limits on the DM particle half-life are at $\tau\sim 75$ Gyr or higher
(depending on the used data combination from e.g. CMB, weak-lensing or
cluster-count observations), which is more than five times the age of the
universe [6, 7, 8, 9, 10]. Note, however, that this constraint can be
substantially relaxed if only part of the DM is assumed to be unstable. In
such a multi-component scenario, the half-life $\tau$ can be significantly
shorter, remaining entirely unconstrained if only $\sim 20$ percent of the DM
budget is allowed to decay [11].

Recent weak-lensing observations from KiDS [2], DES [3], and HSC [12] have
found a surprisingly low clustering amplitude in 2-3 $\sigma$ tension with the
CMB. While this difference might not be large enough to cause real concern, it
is worthwhile to note that a non-stable DM sector could provide a straight-
forward explanation for a reduced clustering signal towards low redshifts [13,
14, 15]. Hence, a decrease of the this tension has indeed been reported for
both one-body [3, 16, 7, 17, 18] and two-body decay scenarios [19, 20, 21, 22,
23, 24].

A more significant tension has emerged between the Hubble parameter obtained
from supernova distance measurements [25, 26] and derived via the CMB
experiment Planck [1]. Different one- and two-body decaying DM scenarios with
both early and late-time decays have been proposed as a potential solution to
the Hubble tension [e.g. 27, 28, 29]. However, to what extend such models may
indeed reconcile CMB with the supernova data [30, 31] remains disputed. A
summary of potential solutions to the Hubble tension, including decaying DM
models, can be found in Ref. [32].

In the present paper, we investigate the possibility that part or all of the
DM consist of an unstable particle decaying into radiation. The model is fully
parametrised via the particle decay rate $\Gamma=\tau^{-1}$ and the fraction
of decaying-to-total dark matter, $f=\Omega_{\rm ddm}/\Omega_{\rm dm}$. The
radiation particle is generally assumed to be part of the dark sector as well,
since decays into photons or other particles of the standard model are
strongly constrained. Note, however, that since we only investigate
gravitational effects caused by the decay process, we do not need to specify
the nature of the radiation species.

In this paper we pursue two main goals. First, we study the effect of dark
matter decays on nonlinear structure formation. We modify the $N$-body code
Pkdgrav3 [33] to include decaying DM both at the particle and the cosmological
background level. With this at hand, we run a suite of simulations and
construct a fitting function describing the suppression effect on the matter
power spectrum. As a second step, we perform a forecast analysis for stage-IV
weak-lensing surveys based on an MCMC analysis with $\Lambda$CDM mock data and
a simulated covariance matrix. In particular, we include baryonic effects
(following the procedure outlined in [34, 35]), which cause a similar
suppression effect on the matter power spectrum and could therefore be
degenerate with the signal from decaying dark matter. With this forecast
exercise we want to assess to what extent it is possible to go beyond the
current limits that come mainly from the Planck CMB survey.

The paper is structured as follows: In Sec. 2 we introduce the decaying DM
model and its effects at the level of the cosmological background. Sec. 3 is
dedicated to the N-body simulations; we describe the modifications to the
code, discuss our simulation suite, and present the fitting function for the
matter power spectrum. In Sec. 4 we present the forecast analysis with
specific focus on the decaying DM parameters before concluding in Sec. 5.

##  2 Decaying dark matter model

In this section we describe the one-body decaying dark matter model and
specify the free parameters added to $\Lambda$CDM. The model consists of both
a stable CDM component and a DDM component that decays into a form of non-
interacting dark radiation (DR). Following the notation of Ref. [7], the DDM
and DR components evolve as

$\displaystyle\dot{\rho}_{\rm ddm}+3\mathcal{H}\rho_{\rm ddm}=-\Gamma
a\rho_{\rm ddm},\hskip 14.22636pt\dot{\rho}_{\rm dr}+4\mathcal{H}\rho_{\rm
dr}=\Gamma a\rho_{\rm ddm},$ (2.1)

where the density evolution of the CDM component is simply given by
$\dot{\rho}_{\rm cdm}+3\mathcal{H}\rho_{\rm cdm}=0$. Note that we assume
conformal time as the variable, with $\mathcal{H}$ being the conformal Hubble
parameter. The DDM decay rate $\Gamma$ is considered a free model parameter
and consists of the inverse of the decay half-life, i.e. $\tau=\Gamma^{-1}$.
The conformal Hubble parameter is obtained via the Friedmann equation

$\mathcal{H}^{2}=8\pi Ga^{2}\sum_{i}\rho_{i}/3,$ (2.2)

where $\rho_{i}$ corresponds to all cosmological components, including CDM,
DDM and DR described above. Solving the closed system of Eqs. (2.1) and (2.2)
yields the full cosmological background evolution.

In order to clarify the notations, in particular with respect to the DDM
abundance, we follow the definitions of Ref. [7] and introduce the quantity

$F_{i}(a)=\frac{\rho_{i}(a)}{\rho_{c}}a^{3(1+w_{i})}$ (2.3)

where $w_{i}$ is the equation-of-state parameter of fluid $i$ and $\rho_{c}$
is the critical density of the universe. With the functions $F_{i}(a)$ at
hand, we can now define a set of parameters $\widetilde{\Omega}_{i}\equiv
F_{i}(a=0)$ as well as $\Omega_{i}\equiv F_{i}(a=1)$. These quantities
describe the abundance of component $i$ at redshift zero. However, while the
former ignores all decay processes, the latter takes them into account. This
means that $\Omega_{\rm cdm}=\widetilde{\Omega}_{\rm cdm}$, while $\Omega_{\rm
ddm}<\widetilde{\Omega}_{\rm ddm}$, and $\Omega_{\rm
dr}>\widetilde{\Omega}_{\rm dr}\equiv 0$.

A second free parameter of our DDM model is the fraction of decaying-to-total
dark matter. It is defined as

$f=\widetilde{\Omega}_{\rm ddm}/(\Omega_{\rm cdm}+\widetilde{\Omega}_{\rm
ddm}).$ (2.4)

The case $f=1$ means that there is no CDM component, i.e. that all of the dark
matter is allowed to decay. This scenario would only agree with current data
if the lifetime of the DDM particle is significantly longer than the age of
the universe. The case of $f\ll 1$, on the other hand, describes a model where
only a subdominant dark component decays during cosmic history while the bulk
of DM remains stable. In this scenario, the decay rate remains largely
unconstrained.

##  3 Cosmological $N$-body simulations

After defining the decaying dark matter scenario, we now discuss the effects
on nonlinear structure formation and provide details about the modification of
the $N$-body code and the simulations performed. At the end of this section,
we show the relative effect of DM decays on the matter power spectrum and
provide a fitting function for the observed power suppression effect. Note
that throughout this section we assume the cosmological parameters
$\Omega_{m}=0.307$, $\Omega_{b}=0.048$, $10^{9}A_{s}=2.43$, $h_{0}=0.678$, and
$n_{s}=0.96$ consistent with the Planck CMB+BAO results from 2018 [1].

###  3.1 Implementation

We implement the DDM model in the $N$-body code Pkdgrav3 following the method
described in Ref. [7]. This means that we decrease the mass of each simulation
particle according to the relation

$m(t)=m_{i}\left[(1-r_{\rm dm}f)+r_{\rm dm}fe^{-\Gamma t}\right],\hskip
14.22636ptr_{\rm dm}=\Omega_{\rm dm}/(\Omega_{b}+\Omega_{\rm dm}),$ (3.1)

where $m_{i}$ is the initial particle mass before the decay starts. Note that
this equation coincides with Eq. (3.1) of Ref. [7], except that we allow for
the additional variation of $f$.

In Pkdgrav3, initial conditions can be generated through the so-called _back-
scaling_ of a linear transfer function given at $z=0$. This grants sub-percent
precision in the linear evolution as well as enabling second-order Zel’Dovich
initial conditions generation. In a DDM scenario the _back-scaling_ approach
is known to produce additional terms when compared to the $\Lambda$CDM
prescription [36]. In principle, it would be possible to account for these
terms and keep the _back-scaling_ approach, however we opted for modifying the
code in a way such that initial conditions are set directly from a given
linear matter power spectrum at the starting redshift of the simulation
(throughout this work $z_{i}=49$). We use CLASS [37] to compute the linear
matter power spectrum starting from an early redshift. Furthermore, we use
first-order Zel’Dovich initial conditions generation (1LPT). This approach can
be referred to as the _forward_ approach and provides a correct prescription
for using $N$-body simulations to simulate non-$\Lambda$CDM scenarios.

###  3.2 Simulation results

Based on the modified code described above, we run a suite of $N$-body DDM and
CDM simulations varying the two DDM parameters $\tau$ and $f$ assuming
$\tau=\\{31.6,\,100,\,316\\}$ Gyr for the decay lifetime, and
$f=\\{0.0,\,0.2,\,0.4,\,0.6,\,0.8,\,1.0\\}$ for the decaying-to-total DM
fraction. Note that $f=0$ corresponds to the standard $\Lambda$CDM case while
$f=1$ means that all DM is allowed to decay.

For each model, we use a box size of $L=500$ Mpc/h and a total particle number
of $N=1024^{3}$. Some selected simulations are also run at lower particle
numbers of $N=256^{3}$, $512^{3}$ in order to test the convergence of our
results. The initial conditions of the simulations are set up with Pkdgrav3 at
$z_{i}=49$, using first order Lagrangian perturbation theory (1LPT) [38]. All
runs are evolved to redshift zero, producing both particle output boxes plus
power spectrum measurements with the Pkdgrav3 internal on-the-fly power
spectrum calculator. A list of all simulations is provided in Table 1.

$N$ |  $L$ [Mpc/$h$] |  $\tau$ [Gyr] | $f$
---|---|---|---
10243 | 500 | 31.6 | 1, 0.8, 0.6, 0.4, 0.2, 0
10243 | 500 | 100 | 1, 0.8, 0.6, 0.4, 0.2, 0
10243 | 500 | 316 | 1, 0.4, 0
5123 | 500 | 100 | 1, 0.3, 0
2563 | 500 | 316 | 1, 0
Table 1: List of $N$-body simulations presented in this work. $N$ is the
particle number, $L$ is the simulation box size, $\tau=1/\Gamma$ the decay
lifetime, and $f$ the decaying-to-total DM fraction. Note that $f=1$ means
that all DM is allowed to decay while $f=0$ refers to the $\Lambda CDM$
simulation.

In Fig. 1 we plot the ratio of the DDM power spectra divided by the standard
CDM case at redshift zero. The left-hand, central, and right-hand panels
correspond to three different decay lifetimes while the coloured solid lines
represent simulation results assuming different decaying-to-total DM
fractions. The dashed lines show the linear results obtained with CLASS.

The resulting power suppression of DDM relative to the standard CDM model is
of non-trivial shape. A first smooth step-like suppression can be observed in
the linear regime at around $k\sim 10^{-2}$ h/Mpc. While this feature is
visible in the results from CLASS, we do not resolve it with our simulations
due to the relatively small box-size of 500 Mpc/h. After a visible plateau at
$k\sim 0.2-2$ h/Mpc, the ratio between DDM and CDM power spectra decreases
again. This second suppression is an effect caused by the nonlinear behaviour
of structure formation. It is therefore only visible in the simulation results
but not in the linear results from CLASS.

A second notable effect visible in Fig. 1 is that both DDM parameters ($\tau$,
$f$) have a direct influence on the overall amplitude of the suppression but
not so much on the scales of the first and second step-like suppression
feature. These steps are manifest at $k\sim 10^{-2}$ and $\sim 10^{-1}$ h/Mpc,
independently of the values assumed for $\tau$ and $f$. Similar conclusions
have been drawn by the authors of Ref.[7]. Note that this particular behaviour
makes it possible to come up with a general fitting function for the DDM power
suppression as we will discuss in Sec. 3.3.

Figure 1: $P_{DDM}(k)/P_{\Lambda CDM}(k)$ for the N-body simulations ran with
Pkdgrav3 at $z=0$ and a varying fraction $f$ (solid lines). We also show the
same ratio for CLASS runs, where the same DDM parameters have been used
(dashed lines).

We investigate the convergence properties as well as comparing our code to
published results from Ref.[39] in Appendix B. Summarizing, the analysis of
the simulation convergence together with the cross-check with the CONCEPT
$N$-body code lets us conclude that our results are reliable down to at least
$k$= 6.4 h/Mpc (and possibly up to somewhat higher values).

###  3.3 Fitting function

In this section we present an updated fitting function for the DDM to CDM
power ratio of the DDM model with free decay lifetime $\tau$ and decaying-to-
total DM fraction $f$. The function is inspired by the fit of Ref. [7],
however, with the fraction $f$ as an additional parameter. We also re-
calibrate the fit to our results. We can write the nonlinear suppression as
$P_{\rm DDM}(k)/P_{\rm\Lambda CDM}(k)=1-\varepsilon_{\rm nonlin}(k)$, with

$\frac{\varepsilon_{\rm nonlin}(k)}{\varepsilon_{\rm
lin}}=\frac{1+a(k/\text{Mpc}^{-1})^{p}}{1+b(k/\text{Mpc}^{-1})^{q}}f,$ (3.2)

where the factors $a$, $b$, $p$, and $q$ are given by

$\displaystyle a(\tau,z)$ $\displaystyle=$ $\displaystyle
0.7208+2.027\left(\frac{\text{Gyr}}{\tau}\right)+3.031\left(\frac{1}{1+1.1z}\right)-0.18,$
$\displaystyle b(\tau,z)$ $\displaystyle=$ $\displaystyle
0.0120+2.786\left(\frac{\text{Gyr}}{\tau}\right)+0.6699\left(\frac{1}{1+1.1z}\right)-0.09,$
$\displaystyle p(\tau,z)$ $\displaystyle=$ $\displaystyle
1.045+1.225\left(\frac{\text{Gyr}}{\tau}\right)+0.2207\left(\frac{1}{1+1.1z}\right)-0.099,$
$\displaystyle q(\tau,z)$ $\displaystyle=$ $\displaystyle
0.9922+1.735\left(\frac{\text{Gyr}}{\tau}\right)+0.2154\left(\frac{1}{1+1.1z}\right)-0.056.$

The function $\varepsilon_{\rm lin}$ describing the main redshift evolution of
the signal is given by

$\varepsilon_{\rm lin}(\tau,z)=\alpha\left(\frac{\rm
Gyr}{\tau}\right)^{\beta}\left(\frac{1}{(0.105z)+1}\right)^{\gamma},$ (3.3)

where $\alpha$, $\beta$, $\gamma$ are functions of $\omega_{b}$, $h_{\rm
ddm}$, and $\omega_{m}=\omega_{b}+\omega_{\rm dm}$ given by

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle(5.323-1.4644u-1.391v)+(-2.055+1.329u+0.8672v)w+$
$\displaystyle(0.2682-0.3509u)w^{2},$ $\displaystyle\beta$ $\displaystyle=$
$\displaystyle 0.9260+(0.05735-0.02690v)w+(-0.01373+0.006713v)w^{2},$
$\displaystyle\gamma$ $\displaystyle=$
$\displaystyle(9.553-0.7860v)+(0.4884+0.1754v)w+(-0.2512+0.07558v)w^{2},$

with $u=\omega_{b}/0.02216$, $v=h_{\rm ddm}/0.6776$ and
$w=\omega_{m}/0.14116$. The fitting function is able to reproduce our
simulation results by better than one percent up to $k=13$ h/Mpc and for
redshifts between 0 and 2.35. This is significantly better than the 10%
accuracy reported by Ref. [7] who only focused on the case of $f=1$.

In Fig. 2 we compare the fitting function to our simulation results, assuming
three different decay lifetimes $\tau=31.6$, 100, and 316 Gyr (left to right).
The top-panels refer to the $z=0$ case with varying decaying-to-total DM
fraction ($f$); the bottom-panels focus on the $f=1$ case showing the redshift
evolution from $z=0$ to 2.35. Note that a similar agreement is obtained for
the redshift evolution of the other models with decaying-to-total DM fraction
below one.

The fit of Eq. 3.2 has been calibrated to simulations based on the Planck
cosmology. However, it is argued in Ref. [7] that varying the cosmological
parameters should only marginally affect the validity of the fit as long as
the parameters do not significantly deviate from the fiducial ones.

Figure 2: Matter power spectrum ratio $P_{DDM}(k)/P_{\Lambda CDM}(k)$ for
N-body simulations run with Pkdgrav3. Top row: simulation results at $z=0$ and
with a varying fraction $f$ parameter for the decay lifetimes $\tau=31.6\
\text{Gyr},\ 100\ \text{Gyr},\ 316\ \text{Gyr}$ (solid lines). Bottom row:
simulation results at $z=0,\ 0.22,\ 0.54,\ 1.08,\ 2.35$ with a fixed fraction
$f=1$ parameter for the decay lifetimes $\tau=31.6\ \text{Gyr},\ 100\
\text{Gyr},\ 316\ \text{Gyr}$ (solid lines). Overlapped we show the fitting
formula described in Sec.3.3 (dashed lines).

##  4 Parameter forecast analysis

After having quantified the nonlinear effects of decaying DM on the power
spectrum, we will now investigate how well DM decays can be constrained with
future galaxy surveys. We specifically focus on the tomographic weak-lensing
shear power spectrum as a potential observable, since weak lensing probes the
medium to small cosmological scales (up to $k\sim 10$ h/Mpc) where the effect
from DDM is strongest.

Our forecast study follows the analysis presented in Refs. [34, 35]. It is
based on a realistic Euclid-like mock data set with three tomographic bins
between redshift 0.1 and 1.5 and with a simulation-based covariance matrix
from statistically independent light-cone realisations that include a
simplified survey footprint, shape-noise errors, and sample variance. The mock
data is based on a $\Lambda$CDM universe and accounts for baryonic effects via
the _baryonification_ model of Refs. [40, 41]. This is particularly important
since baryonic feedback leads to a suppression effect of the power spectrum
that is qualitatively similar to the effect from decaying DM [see e.g. Refs.
42, 43].

The prediction of the weak-lensing shear power spectrum is based on the Limber
approximation and the halofit matter power spectrum from Refs. [44, 45]. The
baryonic effects are included using the emulator of Ref. [34]. The suppression
due to DM decays is accounted for via the fitting function presented in Sec.
3.3. We thereby implicitly assume that DDM effects are separable from both
cosmology and baryonic feedback effects.

Based on the Euclid-like mock data and the prediction-pipeline discussed
above, we perform an MCMC likelihood analysis to estimate the constraining
power for one-body decaying dark matter. We thereby investigate two scenarios,
one where the neutrinos are assumed to be mass-less (Sec. 4.1) and one where
we vary (and marginalise over) neutrino masses between 0.06 and 0.8 eV (Sec.
4.2). While we know that neutrinos have non zero mass, the first scenario
remains interesting, as it allows for a direct comparison with previous
results from the literature. Next to the sum of the neutrino masses ($\Sigma
m_{\nu}$) and the two DDM parameters ($\tau$, $f$), we furthermore vary five
cosmological parameters ($\Omega_{m}$, $\Omega_{b}$, $A_{s}$, $n_{s}$,
$h_{0}$), one intrinsic alignment parameter ($A_{\text{IA}}$, based on the
model of Ref. [46]), plus three baryonic parameters ($M_{c}$, $\mu$,
$\theta_{ej}$) from the _baryonification_ model of Ref. [41]. This means that
most our MCMC runs are based on a total of 11 and 12 free parameters,
respectively.

For the intrinsic alignment and baryonic parameters, we use flat priors that
are wide enough to not affect the obtained posterior results. One notable
exception is the parameter $\mu$ that remains largely unconstrained by our
mock data. Note, however, that in any case, the baryonic priors are selected
to be wide enough so conservatively include a large variety of different
results predicted by hydrodynamical simulations [see Ref. 34, for more
details].

Regarding the cosmological parameters, we first perform a weak-lensing-only
analysis with flat, uninformative priors. This corresponds to the case where
no other information other than the one from the Euclid-type WL mock is
included. In a second step, we also include Gaussian priors from the CMB
experiment Planck assuming the latest Planck18 release [1].

In this context, it is important to notice that the CMB data is not only able
to constrain the standard cosmological parameters, but it also provides
constraints on the DDM parameters themselves. This may be surprising at first,
since recombination happens well before the start of any DM decays (at least
for the models investigated here). However, since the CMB contains information
about the late-time cosmic growth history via the Sachs-Wolfe effect, it is
also sensitive to the decaying DM parameters as they lead to evolving
$\Omega_{m}$ and $\Omega_{r}$. We obtain CMB priors for the DDM parameters
using results from Ref. [11]. The details about how these priors are extracted
are outlined in Appendix A.

Parameters | Mock value | WL priors | WL+CMB priors | WL+CMB+fixed baryons priors
---|---|---|---|---
$\Omega_{b}$ | 0.049 | 0.04 - 0.06 | 0.001 | 0.001
$\Omega_{m}$ | 0.315 | 0.25 - 0.42 | 0.0084 | 0.084
$10^{9}A_{s}$ | 2.035 | 1.4 - 3.5 | 0.034 | 0.034
$n_{s}$ | 0.966 | 0.9 - 1.0 | 0.0044 | 0.0044
$h_{0}$ | 0.673 | 0.6 - 0.9 | 0.006 | 0.006
$\Sigma m_{\nu}$ (Sec. 4.1) | 0 | 0 | 0 | 0
$\Sigma m_{\nu}$ (Sec. 4.2) | 0.1 | 0.06 - 0.8 | 0.1 | 0.1
$A_{IA}$ | 1 | 0 - 2 | 0 - 2 | 0 - 2
$\log M_{c}$ | 13.8 | 13 - 16 | 13 - 16 | 0
$\mu$ | 0.21 | 0.1 - 0.7 | 0.1-0.7 | 0
$\theta_{ej}$ | 4 | 2 - 8 | 2 - 8 | 0
$\Gamma$ | 0 |  10-6 \- 0.0316 | see below | see below
$f$ | 0 | 0.0 - 1.0 | see below | see below
$\Gamma f^{c_{3}}$ | 0 | not used | $(c_{2}+c_{1}f^{c_{3}})/2$ | $(c_{2}+c_{1}f^{c_{3}})/2$
Table 2: Priors of cosmological, neutrino, intrinsic-alignment ($A_{IA}$),
baryonic parameters, and DDM parameters for the WL, WL+CMB, and WL+CMB+fixed
baryons analysis. Whenever we assume a flat prior, we specify the range (lower
limit - upper limit). For Gaussian priors, we provide just one number
(standard deviation). For the DDM Gaussian priors we follow the approach
outlined in Appendix A, using $c_{1}=0.00721714,\ c_{2}=0.00308819,\
c_{3}=1.45599761$. Note that neutrino masses are assumed fixed and massless in
Sec. 4.1 and only varied in Sec. 4.2.

###  4.1 Forecast assuming massless neutrinos

As a first step, we perform an MCMC analysis with the sum of the neutrino
masses fixed at zero ($\Sigma m_{\nu}=0$). While this assumption is known to
be inadequate, it considerably simplifies the analysis and allows us to
directly compare to previous work on decaying DM [which did not vary neutrino
masses, see e.g. Refs. 7, 11, 8].

Figure 3: Results for our MCMC analysis assuming massless neutrinos. In black
we show results when assuming exclusively WL observables ($WL$), in red we
additionally assume CMB priors ($WL$+CMB) and in green we additionally fix the
baryonic parameters $log(M_{c}),\ \mu,\ \theta_{ej}$ to their mock value. We
report results from Ref.[11] (Poulin et al. 2016) in blue. Contours in the
left panel are reported at 95$\%$ CL In the right-hand panel we show
marginalised errors ($\sigma_{\theta_{i}}$) for all the parameters, normalized
by their mock value ($\theta_{i}$), reported on 68$\%$ CL

We start by performing an MCMC run only considering the WL mock data and
assuming flat, uninformative priors (as specified in the third column of table
2). This corresponds to the hypothetical case where Euclid-like WL data is
used as an independent probe with no added information — neither from galaxy
clustering nor from any external data source.

The resulting constraints on the DDM parameters for the WL-only case are shown
as grey contours in the left-hand panel of Fig. 3. Comparing it to the current
constraints based on the Planck data from Ref. [11] (shown as blue line), we
conclude that the Euclid-like WL data will have a similar constraining power
than Planck for the $f=1$ case, while it will lead to weaker constraints for
$f<1$.

As a further step, we run an MCMC analysis combining Euclid-like weak-lensing
data with CMB data from Planck. We do this by adding Gaussian priors from the
Planck18 data set to both the cosmological and the DDM parameters. The priors
of the DDM parameters are obtained from Ref. [11] following the method
outlined in Appendix A. Detailed information about the width of the Gaussian
priors can be found in the fourth column of table. 2. The forecast posteriors
of the combined WL+CMB analysis is shown as red contours in the left-hand
panel of Fig. 3. The contours are about 50 percent tighter than the current
constraints from Ref. [11] (blue line) for all values of $f$.

So far, we have marginalised over baryonic parameters assuming that we have no
prior knowledge of the baryonic suppression effect on the matter power
spectrum. In Ref. [35], however, it was shown that future X-ray data may be
used to strongly constrain the baryonic parameters, e.g. via the observed
fraction of gas in haloes. In order to estimate how much can be gained by
including additional information about baryons, we run another MCMC chain
fixing the baryonic parameters to their default values.

The resulting posteriors of the WL+CMB analysis with fixed baryons is shown as
green contours in the left-hand panel of Fig. 3. They are 40 percent tighter
compared to the case with free baryonic parameters and 60 percent tighter
compared to current constraints from Planck. As there is a strong degeneracy
between the DDM parameters and the baryonic feedback parameters, fixing the
latter leads to a significant improvement in constraining power. It highlights
the importance of better understanding and constraining baryonic feedback for
future WL surveys.

Regarding the special case of $f=1$ (i.e. all DM being allowed to decay), we
obtain $\tau\geq 89.08$ Gyr for the WL only case, $\tau\geq 154.19$ Gyr for
the WL+CMB case, and $\tau\geq 289.28$ Gyr for the WL+CMB case with fixed
baryonic parameters. For the limit of very large decay rates ($\tau<34$ Gyr),
on the other hand, we obtain limits on the decaying-to-total DM fraction of
$f<0.38$ for the WL-only case, $f<0.17$ for the WL+CMB case, and $f<0.12$ for
WL+CMB case with fixed baryonic parameters. Note that all these limits are at
95$\%$ CL.

Ref.[11] reports a lower bound of 3.8% for the fraction of DDM that has
decayed between recombination and today (at 95% CL). We emphasize that this is
not the same constraint that we obtain on $f$ when observing the limit of
small decay lifetimes (our smallest bin being $31.6$ Gyr $<\tau<34$ Gyr). In
our case the fraction $f$ represents the decaying-to-total DM fraction at the
beginning of the simulation. Tighter constraints on $f$ could be achieved by
exploring cases with small decay lifetimes (smaller than the age of the
Universe), a domain we do not explore in this work.

The right-hand panel of Fig. 3 summarises the marginalised constraints (now at
68$\%$ CL) of all parameters in our MCMC forecast analysis. It shows that
adding prior information from the CMB significantly decreases the
uncertainties of all cosmological and DDM parameters. Unsurprisingly, the
largest improvements are obtained for the cosmological parameters $h_{0}$ and
$\Omega_{b}$, which are known to be rather poorly constrained by WL data.
Additionally fixing baryonic parameters only has a relatively mild effect (of
order 20 percent) on the cosmological parameters [as pointed out in Refs. 34,
35]. However, it significantly affects the two DDM parameters whose 68$\%$ CL
constraints become tighter by more than a factor of 2. This strong improvement
is a result of the fact that baryonic and DDM parameters affect the matter
power spectrum in similar ways and are therefore inherently degenerate.

The results presented in this section indicate that with upcoming WL surveys
such as Euclid or LSST, it will be possible to go significantly beyond current
constraints, probing unexplored parts of the DDM parameter space. In the next
section we will investigate to what extent this conclusion still holds if we
vary the neutrino masses as an additional free model parameter.

###  4.2 Forecast with varying neutrino masses

Upcoming weak-lensing and galaxy clustering surveys will be sensitive enough
to probe and constrain the sum of the neutrino masses [see e.g. Ref. 47]. We
therefore repeat the analysis above, now varying and marginalising over the
$\Sigma m_{\nu}$ parameter. While this is a more consistent approach, it is
worth noticing that the results cannot be directly compared to previous limits
that were obtained assuming neutrinos to be mass-less [e.g. 7, 8, 11, 3].

Figure 4: Results for the MCMC analysis with varying neutrino masses. The
weak-lensing only case (WL) is shown in grey, the scenarios with additional
CMB priors (WL+CMB) and with CMB priors plus fixed baryonic parameters
(WL+CMB+fixed baryons) are shown in red and green. The CMB priors are obtained
from the Planck18 analysis [1]. The left-hand panel shows the posteriors of
the decaying DM parameters (at 95 percent CL) with all other parameters
marginalised over. The right-hand panel provides an overview over all
marginalised errors ($\sigma_{\theta_{i}}$) at 68 percent CL normalized to
their mock values ($\theta_{i}$). For the special cases of $\tau$ and $f$
(where the mock values are zero), we show the absolute errors at the limits
$f=1$ and $\tau=32.8$ Gyr, instead.

Since the true value of $\Sigma m_{\nu}$ is not known, we first have to select
a reasonable mass for our mock data. We assume $\Sigma m_{\nu}=0.1$ eV based
on a normal mass hierarchy, a setup that is well within the current neutrino
mass constraints [1, 48]. More details about our mock data for a $\Lambda$CDM
case with massive neutrinos can be found in Ref. [35].

The performed MCMC runs are based on the same setup as in the previous
section. We first assume a WL-only case before adding priors for the
cosmological, DDM, and neutrino parameters from the Planck18 CMB experiment
(WL+CMB case). Finally, we investigate a third scenario where baryonic
parameters are fixed, mimicking the ideal case of strongly constrained
baryonic feedback from e.g. X-ray observations [34, 35].

In the left-hand panel of Fig. 4 we plot the posteriors of the decaying DM
parameters including varying neutrino masses. The contours of the three
scenarios are again shown in grey (WL-only), red (WL+CMB), and green (WL+CMB
with fixed baryonic parameters) at 95 percent CL. Compared to the results with
fixed neutrino masses (see Fig. 3), the contours of the WL only case are
somewhat looser, while the WL+CMB cases (without and with fixed baryonic
parameters) remain very similar. Note that although we generally expect limits
to become weaker when adding a new model parameter, this does not have to be
the case here since we have also changed the mock data (now assuming a
$\Lambda$CDM model with $\Sigma m_{\nu}=0.1$ eV).

Regarding the special case of $f=1$, we obtain constraints on the decay
lifetime of $\tau\geq 74.88$ Gyr for the case of WL only, $\tau\geq 182.23$
Gyr when additionally introducing CMB priors and $\tau\geq 275$ Gyr when the
baryonic parameters are fixed. Considering the limit of large decay rates and
small decaying-to-total DM fraction, we obtain a maximum allowed fraction of
$f<0.49$ for the WL only case, $f<0.21$ for WL+CMB case, and $f<0.13$ for the
WL+CMB case with fixed baryonic parameters.

In the right-hand panel of Fig. 4 we show the marginalised errors (at 68
percent CL) of all cosmological, intrinsic alignment, baryonic, and DDM
parameters. As in the previous case, we notice that the errors on $h_{0}$ and
$\Omega_{b}$ improve most when combining WL with CMB priors. We also confirm
that by fixing the baryonic parameters, the constraints on the DDM parameters
improve by about another factor of 2. The constraint on the neutrino parameter
($\Sigma m_{\nu}$) on the other hand, only slightly improves when fixing the
baryonic parameters. This suggests that neutrino and baryonic parameters are
not degenerate, a fact that has been pointed out previously [49, 50, 35].

Overall, we confirm that the general conclusions obtained with the assumption
of mass-less neutrinos also hold for the more general case of massive
neutrinos with varying $\Sigma m_{\nu}$. Although the constraints on the DDM
parameters become looser for the WL only case, this difference shrinks
substantially once CMB priors are included. For the WL+CMB cases, fixing or
marginalising over the neutrino mass leads to very similar DDM constraints. We
therefore conclude that even in the realistic case of a $\Lambda$CDM universe
with massive neutrinos, the the one-body decaying DM model can be efficiently
constrained by combining upcoming WL observations with data from the CMB.

##  5 Conclusions

In this paper we have investigated the nonlinear effects on the matter power
spectrum, assuming a decaying dark matter scenario where parts or all of the
dark matter may decay into (dark) radiation. We have modified the $N$-body
code Pkdgrav3, [following an approach described in Ref. 7] and run a suite of
simulations with varying decay lifetime ($\tau$) and decaying-to-total DM
fraction ($f$). Based on these simulations, we defined a fitting function for
the relative effect of decaying DM on the matter power spectrum. The fit is
sub-percent accurate up to $k=12$ h/Mpc, making it a promising tool to probe
DM decays with future weak-lensing and galaxy-clustering surveys.

With the fit of the nonlinear power spectrum from simulations, we performed a
forecast analysis for a Euclid-like tomographic weak-lensing survey using
$\Lambda$CDM mock data and a simulation-based covariance matrix [following the
method outlined in Refs. 34, 35]. We first ran an MCMC analysis,
simultaneously varying cosmological, intrinsic alignment, baryonic, and
decaying DM parameters, while assuming neutrinos to be mass-less. With this
setup in hand, we showed that, while Euclid-like weak-lensing data with flat,
uninformative priors will lead to weaker limits on decaying DM compared to
current constraints from the CMB [11], a combination of WL with CMB data will
make it possible to significantly push into currently unexplored DDM
parameters space. This is even more the case for the optimistic scenario where
baryonic feedback effects can be fully determined using external data [e.g.
from the X-ray survey eROSITA, see Ref. 35]. In this case we expect an
improvement on the DDM parameter constraints (for $\tau$ and $f$) that are at
least a factor of 2-3 better than current constraints from the CMB [11].

As a further step, we ran the same MCMC analysis, this time assuming mock data
based on a $\Lambda$CDM universe with massive neutrinos (of $\Sigma
m_{\nu}=0.1$ eV) and adding the sum of the neutrino mass as a free model
parameter. While this is a more realistic scenario, it does not allow for a
direct comparison to previous constraints where neutrinos were assumed to be
mass-less. The forecast analysis with massive neutrinos shows similar results
compared to the mass-less case with only a $\sim 5$ percent degradation of the
DDM parameter constraints. Assuming a scenario where all DM is allowed to
decay ($f=1$), we forecast a minimum decay time of $\tau\geq 182.23$ Gyr and
$\tau\geq 275$ Gyr, respectively. The former value corresponds to the
pessimistic case where baryonic parameters are marginalised over; the latter
value shows the optimistic scenario where baryonic parameters are fixed.
Assuming the limit where only a small fraction of DM is allowed to rapidly
decay ($\tau<34$ Gyr), we show that WL+CMB data will be able to rule-out
decaying-to-total DM fractions above $f=0.21$ and $f=0.13$, respectively,
depending on whether baryonic parameters are marginalised over or fixed.

We conclude that upcoming weak-lensing surveys, such as Euclid or LSST at the
Rubin Observatory, will allow us to significantly improve upon current
constraints for the one-body decaying DM scenario. Compared to current limits
from the CMB, an improvement of about a factor of 2-3 can be expected.
However, this will require additional knowledge regarding the effects of
baryonic feedback on the clustering at cosmological scales.

## Acknowledgments

This study is supported by the Swiss National Science Foundation via the grant
PCEFP2_181157.

##  Appendix A CMB priors for the decaying dark matter parameters.

In order to assume CMB priors for the DDM parameters in our forecast analysis,
we define a prior for the parameter combination $\Gamma f^{c_{3}}$, where
$c_{3}$ is a constant. We emphasize that in the MCMC analysis, the two
parameters $\Gamma$ and $f$ are separate, but we apply a single prior on the
combination $\Gamma f^{c_{3}}$. Starting from the constraints obtained in Ref.
[11] reported here in Fig. 5 in blue, we find a fit for the 95$\%$ CL contour
with form $\Gamma=c_{1}+c_{2}/f^{c_{3}}$ (shown in Fig. 5 in solid black).
Having obtained the fit, we can check the 68$\%$ CL level by replacing
$c_{1},\ c_{2}$ with $c_{1}/2,\ c_{2}/2$ (shown in dashed black in Fig. 5). We
notice that in our construction at the 68$\%$ CL our fit is more conservative
than the results from Ref.[11]. The 68$\%$ CL line will therefore be the
standard deviation of the gaussian prior for the parameter combination $\Gamma
f^{c_{3}}$ with mean value 0. Overall, this creates a likelihood distribution
which is peaked when $f=0$ or $\Gamma\xrightarrow[]{}0$ and then falls as
shown by the 68$\%$ CL and 95$\%$ CL contours shown in Fig. 5. The values of
the fit constants are $(c_{1},\ c_{2},\ c_{3})=(0.00721714,\ 0.00308819,\
1.45599761)$.

Figure 5: Results from Poulin et al. 2016 [11] (shown in light and dark blue
at 95% CL and 68% CL respectively) and our constructed prior on $\Gamma
f^{c_{3}}$ (68% and 95% levels shown in dashed and solid lines).

##  Appendix B Convergence tests and comparison with the $N$-body code CONCEPT

In this appendix we investigate the convergence of the simulations as well as
comparing our results to another $N$-body code, CONCEPT [39]. In order to
investigate the convergence properties of our simulations, we show the DDM to
CDM power ratio for three simulations with different particle numbers in the
left-hand and central panel of Fig. 6. We thereby limit ourselves to the two
cases of $\tau=100$ (left) and $316$ Gyr (right) both with $f=1$. The power
ratios from the three different simulations are well converged agreeing at the
sub-percent level at least until the Nyquist frequency ($\pi N^{1/3}/L$). Note
that this is only true for the ratio plot but not for the absolute power
spectrum, where higher resolutions and larger box-sizes are required.

A further test of our $N$-body code and the simulation runs is presented in
the right-hand panel of Fig. 6. There we compare the DDM-to-CDM power spectrum
ratio for a model with $\tau=100$ Gyr and $f=0.3$ to the published results
from Ref. [39] based on the $N$-body code CONCEPT. We find sub-percent
agreement between the results from the two codes at wave-modes below $k=1$
h/Mpc. The small deviations beyond that scale is most likely due to limited
resolution of the simulation with CONCEPT.

Figure 6: Matter power spectrum ratio $P_{DDM}(k)/P_{\Lambda CDM}(k)$ for
N-body simulations ran with Pkdgrav3 with $N=256^{3},\ 512^{3},\ 1024^{3}$
particles for decay lifetimes $\tau=316\ \text{Gyr}$, $\tau=100\ \text{Gyr}$
and a fixed fraction $f=1$ (left and center panel). In the right panel we show
results from [39] (CONCEPT, solid red) and our simulations run in Pkdgrav3
(solid black) with $\tau$=100 Gyr and $f$=0.3, 5123 particles and a simulation
box size of 500 Mpc/h.

## References

  * [1] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, and et al.  Planck 2018 results.  Astronomy & Astrophysics, 641:A6, Sep 2020.
  * [2] Giblin, Benjamin, Heymans, Catherine, Asgari, Marika, Hildebrandt, Hendrik, Hoekstra, Henk, Joachimi, Benjamin, Kannawadi, Arun, Kuijken, Konrad, Lin, Chieh-An, Miller, Lance, Tröster, Tilman, van den Busch, Jan Luca, Wright, Angus H., Bilicki, Maciej, Blake, Chris, de Jong, Jelte, Dvornik, Andrej, Erben, Thomas, Getman, Fedor, Napolitano, Nicola R., Schneider, Peter, Shan, HuanYuan, and Valentijn, Edwin.  Kids-1000 catalogue: Weak gravitational lensing shear measurements.  A&A, 645:A105, 2021.
  * [3] T.M.C. Abbott, F.B. Abdalla, S. Avila, M. Banerji, E. Baxter, K. Bechtol, M.R. Becker, E. Bertin, J. Blazek, S.L. Bridle, and et al.  Dark energy survey year 1 results: Constraints on extended cosmological models from galaxy clustering and weak lensing.  Physical Review D, 99(12), Jun 2019.
  * [4] Kyle S. Dawson, David J. Schlegel, Christopher P. Ahn, Scott F. Anderson, Éric Aubourg, Stephen Bailey, Robert H. Barkhouser, Julian E. Bautista, Alessandra Beifiori, Andreas A. Berlind, and et al.  The baryon oscillation spectroscopic survey of sdss-iii.  The Astronomical Journal, 145(1):10, Dec 2012.
  * [5] V. Bonvin, F. Courbin, S. H. Suyu, P. J. Marshall, C. E. Rusu, D. Sluse, M. Tewes, K. C. Wong, T. Collett, C. D. Fassnacht, T. Treu, M. W. Auger, S. Hilbert, L. V. E. Koopmans, G. Meylan, N. Rumbaugh, A. Sonnenfeld, and C. Spiniello.  H0LiCOW – V. New COSMOGRAIL time delays of HE 0435-1223: H0 to 3.8 percent precision from strong lensing in a flat $\Lambda$CDM model.  Monthly Notices of the Royal Astronomical Society, 465(4):4914–4930, 11 2016.
  * [6] Kiyotomo Ichiki, Masamune Oguri, and Keitaro Takahashi.  WMAP constraints on decaying cold dark matter.  Phys. Rev. Lett., 93:071302, 2004.
  * [7] Kari Enqvist, Seshadri Nadathur, Toyokazu Sekiguchi, and Tomo Takahashi.  Decaying dark matter and the tension in $\sigma_{8}$.  Journal of Cosmology and Astroparticle Physics, 2015, 05 2015.
  * [8] Kari Enqvist, Seshadri Nadathur, Toyokazu Sekiguchi, and Tomo Takahashi.  Constraints on decaying dark matter from weak lensing and cluster counts.  Journal of Cosmology and Astroparticle Physics, 2020(04):015–015, Apr 2020.
  * [9] Angela Chen, Dragan Huterer, Sujeong Lee, Agnès Ferté, Noah Weaverdyck, Otavio Alonso Alves, C. Danielle Leonard, Niall MacCrann, Marco Raveri, Anna Porredon, Eleonora Di Valentino, Jessica Muir, Pablo Lemos, Andrew Liddle, Jonathan Blazek, Andresa Campos, Ross Cawthon, Ami Choi, Scott Dodelson, Jack Elvin-Poole, Daniel Gruen, Ashley Ross, Lucas F. Secco, Ignacio Sevilla, Erin Sheldon, Michael A. Troxel, Joe Zuntz, Tim Abbott, Michel Aguena, Sahar Allam, James Annis, Santiago Avila, Emmanuel Bertin, Sunayana Bhargava, Sarah Bridle, David Brooks, Aurelio Carnero Rosell, Matias Carrasco Kind, Jorge Carretero, Matteo Costanzi, Martin Crocce, Luiz da Costa, Maria Elidaiana da Silva Pereira, Tamara Davis, Peter Doel, Tim Eifler, Ismael Ferrero, Pablo Fosalba, Josh Frieman, Juan Garcia-Bellido, Enrique Gaztanaga, David Gerdes, Robert Gruendl, Julia Gschwend, Gaston Gutierrez, Samuel Hinton, Devon L. Hollowood, Klaus Honscheid, Ben Hoyle, David James, Mike Jarvis, Kyler Kuehn, Ofer Lahav, Marcio Maia, Jennifer Marshall, Felipe Menanteau, Ramon Miquel, Robert Morgan, Antonella Palmese, Francisco Paz-Chinchon, Andrés Plazas Malagón, Aaron Roodman, Eusebio Sanchez, Vic Scarpine, Michael Schubnell, Santiago Serrano, Mathew Smith, Eric Suchyta, Gregory Tarle, Daniel Thomas, Chun-Hao To, Tamas Norbert Varga, Jochen Weller, and Reese Wilkinson.  Constraints on decaying dark matter with des-y1 and external data, 2020\.
  * [10] Benjamin Audren, Julien Lesgourgues, Gianpiero Mangano, Pasquale Dario Serpico, and Thomas Tram.  Strongest model-independent bound on the lifetime of dark matter.  Journal of Cosmology and Astroparticle Physics, 2014(12):028–028, Dec 2014.
  * [11] Vivian Poulin, Pasquale D. Serpico, and Julien Lesgourgues.  A fresh look at linear cosmological constraints on a decaying dark matter component.  Journal of Cosmology and Astroparticle Physics, 2016(08):036–036, aug 2016.
  * [12] Chiaki Hikage, Masamune Oguri, Takashi Hamana, Surhud More, Rachel Mandelbaum, Masahiro Takada, Fabian Köhlinger, Hironao Miyatake, Atsushi J Nishizawa, Hiroaki Aihara, Robert Armstrong, James Bosch, Jean Coupon, Anne Ducout, Paul Ho, Bau-Ching Hsieh, Yutaka Komiyama, François Lanusse, Alexie Leauthaud, Robert H Lupton, Elinor Medezinski, Sogo Mineo, Shoken Miyama, Satoshi Miyazaki, Ryoma Murata, Hitoshi Murayama, Masato Shirasaki, Cristóbal Sifón, Melanie Simet, Joshua Speagle, David N Spergel, Michael A Strauss, Naoshi Sugiyama, Masayuki Tanaka, Yousuke Utsumi, Shiang-Yu Wang, and Yoshihiko Yamada.  Cosmology from cosmic shear power spectra with Subaru Hyper Suprime-Cam first-year data.  Publications of the Astronomical Society of Japan, 71(2), 03 2019\.  43\.
  * [13] Mei-Yu Wang and Andrew R. Zentner.  Weak gravitational lensing as a method to constrain unstable dark matter.  Physical Review D, 82(12), Dec 2010.
  * [14] Niall MacCrann, Joe Zuntz, Sarah Bridle, Bhuvnesh Jain, and Matthew R. Becker.  Cosmic discordance: Are planck cmb and cfhtlens weak lensing measurements out of tune?, 2015.
  * [15] Renyue Cen.  Decaying cold dark matter model and small-scale power.  The Astrophysical Journal, 546(2):L77–L80, jan 2001.
  * [16] Z. Berezhiani, A. D. Dolgov, and I. I. Tkachev.  Reconciling Planck results with low redshift astronomical measurements.  Phys. Rev. D, 92(6):061303, 2015.
  * [17] A. Chudaykin, D. Gorbunov, and I. Tkachev.  Dark matter component decaying after recombination: Lensing constraints with Planck data.  Phys. Rev. D, 94:023528, 2016.
  * [18] Maria Archidiacono, Deanna C. Hooper, Riccardo Murgia, Sebastian Bohr, Julien Lesgourgues, and Matteo Viel.  Constraining dark matter-dark radiation interactions with cmb, bao, and lyman-.  Journal of Cosmology and Astroparticle Physics, 2019(10):055–055, Oct 2019.
  * [19] Guillermo F. Abellan, Riccardo Murgia, Vivian Poulin, and Julien Lavalle.  Hints for decaying dark matter from $s_{8}$ measurements, 2020.
  * [20] Dalong Cheng, M.-C. Chu, and Jiayu Tang.  Cosmological structure formation in decaying dark matter models.  Journal of Cosmology and Astroparticle Physics, 2015(07):009–009, jul 2015.
  * [21] Shohei Aoyama, Kiyotomo Ichiki, Daisuke Nitta, and Naoshi Sugiyama.  Formulation and constraints on decaying dark matter with finite mass daughter particles.  Journal of Cosmology and Astroparticle Physics, 2011(09):025–025, Sep 2011.
  * [22] Annika H. G. Peter.  Mapping the allowed parameter space for decaying dark matter models.  Physical Review D, 81(8), Apr 2010.
  * [23] Annika H. G. Peter, Christopher E. Moody, and Marc Kamionkowski.  Dark-matter decays and self-gravitating halos.  Physical Review D, 81(10), May 2010.
  * [24] Mei-Yu Wang, Annika H. G. Peter, Louis E. Strigari, Andrew R. Zentner, Bryan Arant, Shea Garrison-Kimmel, and Miguel Rocha.  Cosmological simulations of decaying dark matter: implications for small-scale structure of dark matter haloes.  Monthly Notices of the Royal Astronomical Society, 445(1):614–629, Sep 2014.
  * [25] Adam G. Riess et al.  A 2.4% Determination of the Local Value of the Hubble Constant.  Astrophys. J., 826(1):56, 2016.
  * [26] Adam G. Riess, Stefano Casertano, Wenlong Yuan, Lucas M. Macri, and Dan Scolnic.  Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM.  Astrophys. J., 876(1):85, 2019.
  * [27] Kanhaiya L. Pandey, Tanvi Karwal, and Subinoy Das.  Alleviating the $H_{0}$ and $\sigma_{8}$ anomalies with a decaying dark matter model.  JCAP, 07:026, 2020.
  * [28] Nikita Blinov, Celeste Keith, and Dan Hooper.  Warm Decaying Dark Matter and the Hubble Tension.  JCAP, 06:005, 2020.
  * [29] Kyriakos Vattis, Savvas M. Koushiappas, and Abraham Loeb.  Dark matter decaying in the late Universe can relieve the H0 tension.  Phys. Rev. D, 99(12):121302, 2019.
  * [30] Steven J. Clark, Kyriakos Vattis, and Savvas M. Koushiappas.  Cosmological constraints on late-Universe decaying dark matter as a solution to the $H_{0}$ tension.  Phys. Rev. D, 103(4):043014, 2021.
  * [31] Balakrishna S. Haridasu and Matteo Viel.  Late-time decaying dark matter: constraints and implications for the $H_{0}$-tension.  Mon. Not. Roy. Astron. Soc., 497(2):1757–1764, 2020.
  * [32] Eleonora Di Valentino, Olga Mena, Supriya Pan, Luca Visinelli, Weiqiang Yang, Alessandro Melchiorri, David F. Mota, Adam G. Riess, and Joseph Silk.  In the Realm of the Hubble tension $-$ a Review of Solutions.  3 2021.
  * [33] Douglas Potter, Joachim Stadel, and Romain Teyssier.  Pkdgrav3: Beyond trillion particle cosmological simulations for the next era of galaxy surveys, 2016.
  * [34] Aurel Schneider, Nicola Stoira, Alexandre Refregier, Andreas J. Weiss, Mischa Knabenhans, Joachim Stadel, and Romain Teyssier.  Baryonic effects for weak lensing. part i. power spectrum and covariance matrix.  Journal of Cosmology and Astroparticle Physics, 2020(04):019–019, apr 2020.
  * [35] Aurel Schneider, Alexandre Refregier, Sebastian Grandis, Dominique Eckert, Nicola Stoira, Tomasz Kacprzak, Mischa Knabenhans, Joachim Stadel, and Romain Teyssier.  Baryonic effects for weak lensing. part II. combination with x-ray data and extended cosmologies.  Journal of Cosmology and Astroparticle Physics, 2020(04):020–020, apr 2020.
  * [36] Christian Fidler, Thomas Tram, Cornelius Rampf, Robert Crittenden, Kazuya Koyama, and David Wands.  Relativistic initial conditions for n-body simulations.  Journal of Cosmology and Astroparticle Physics, 06(043), June 2017\.  20 pages, 8 figures, 12 months embargo.
  * [37] Diego Blas, Julien Lesgourgues, and Thomas Tram.  The cosmic linear anisotropy solving system (CLASS). part II: Approximation schemes.  Journal of Cosmology and Astroparticle Physics, 2011(07):034–034, jul 2011.
  * [38] Darren S. Reed, Robert E. Smith, Doug Potter, Aurel Schneider, Joachim Stadel, and Ben Moore.  Towards an accurate mass function for precision cosmology.  Monthly Notices of the Royal Astronomical Society, 431(2):1866–1882, 03 2013.
  * [39] Jeppe Dakin, Steen Hannestad, and Thomas Tram.  Fully relativistic treatment of decaying cold dark matter in n-body simulations.  Journal of Cosmology and Astroparticle Physics, 2019(06):032–032, jun 2019.
  * [40] Aurel Schneider and Romain Teyssier.  A new method to quantify the effects of baryons on the matter power spectrum.  Journal of Cosmology and Astroparticle Physics, 2015(12):049–049, Dec 2015.
  * [41] Aurel Schneider, Romain Teyssier, Joachim Stadel, Nora Elisa Chisari, Amandine M.C. Le Brun, Adam Amara, and Alexandre Refregier.  Quantifying baryon effects on the matter power spectrum and the weak lensing shear correlation.  Journal of Cosmology and Astroparticle Physics, 2019(03):020–020, Mar 2019.
  * [42] Marcel P van Daalen, Ian G McCarthy, and Joop Schaye.  Exploring the effects of galaxy formation on matter clustering through a library of simulation power spectra.  Monthly Notices of the Royal Astronomical Society, 491(2):2424–2446, Nov 2019.
  * [43] Nora Elisa Chisari, Alexander J. Mead, Shahab Joudaki, Pedro G. Ferreira, Aurel Schneider, Joseph Mohr, Tilman Tröster, David Alonso, Ian G. McCarthy, Sergio Martin-Alvarez, and et al.  Modelling baryonic feedback for survey cosmology.  The Open Journal of Astrophysics, 2(1), Jun 2019.
  * [44] Ryuichi Takahashi, Masanori Sato, Takahiro Nishimichi, Atsushi Taruya, and Masamune Oguri.  Revisiting the halofit model for the nonlinear matter power spectrum.  The Astrophysical Journal, 761(2):152, Dec 2012.
  * [45] Simeon Bird, Matteo Viel, and Martin G. Haehnelt.  Massive neutrinos and the non-linear matter power spectrum.  Monthly Notices of the Royal Astronomical Society, 420(3):2551–2561, 02 2012.
  * [46] Christopher M. Hirata and Uro š Seljak.  Intrinsic alignment-lensing interference as a contaminant of cosmic shear.  Phys. Rev. D, 70:063526, Sep 2004.
  * [47] Luca Amendola et al.  Cosmology and fundamental physics with the Euclid satellite.  Living Rev. Rel., 21(1):2, 2018.
  * [48] Nathalie Palanque-Delabrouille et al.  Neutrino masses and cosmology with Lyman-alpha forest power spectrum.  JCAP, 11:011, 2015.
  * [49] Benjamin O. Mummery, Ian G. McCarthy, Simeon Bird, and Joop Schaye.  The separate and combined effects of baryon physics and neutrino free-streaming on large-scale structure.  Mon. Not. Roy. Astron. Soc., 471(1):227–242, 2017.
  * [50] Gabriele Parimbelli, Matteo Viel, and Emiliano Sefusatti.  On the degeneracy between baryon feedback and massive neutrinos as probed by matter clustering and weak lensing.  JCAP, 01:010, 2019.


* [16] A. Debussche, J. Vovelle;  Scalar conservation laws with stochastic forcing. J. Funct. Anal. 259(4), 1014–1042 (2010).
  * [17] A. Debussche, M. Hofmanová, J. Vovelle, Degenerate parabolic stochastic partial differential equations: quasilinear case. Ann. Probab. 44 (2016), no. 3, 1916–1955. 60H15 (35K65 35R60).
  * [18] L. Denis, A. Matoussi, and L. Stoica, (2005). $L^{p}$ estimate for the uniform norm of solutions of quasilinear SPDE’s. Probab. Theory Related Fields 133 147-463. MR2197109
  * [19] S. Dotti, J. Vovelle; Convergence of Approximations to Stochastic Scalar Conservation Laws, Archive for Rational Mechanics and Analysis, Springer Verlag, 2018, 230 (2), pp.539-591.
  * [20] J. Feng and D. Nualart; Stochastic scalar conservation laws. J. Funct. Anal., 255(2): 313-373, 2008.
  * [21] M. Hofmanová;  Degenerate parabolic stochastic partial differential equations. Stochastic Process. Appl. 123 4294-4336, MR3096355.
  * [22] M. Hofmanova, U. Koley, and U. Sarkar:  Measure-valued solutions to the stochastic compressible Euler equations and incompressible limits.  https://arxiv.org/pdf/2012.07391.pdf
  * [23] J. U. Kim; On a stochastic scalar conservation law. Indiana Univ. Math. J., 52 (1), 227-256, 2003.
  * [24] U. Koley, A. K. Majee, and G. Vallet.  A finite difference scheme for conservation laws driven by Lévy noise.  IMA J. Numer. Anal., 38(2), 998–1050, 2018.
  * [25] U. Koley, A. K. Majee, and G. Vallet; Continuous dependence estimate for a degenerate parabolic-hyperbolic equation with Levy noise. Stoch. Partial Differ. Equ. Anal. Comput., 5 (2), 145-191, 2017.
  * [26] U. Koley, N. H. Risebro, C. Schwab and F. Weber.  A multilevel Monte Carlo finite difference method for random scalar degenerate convection-diffusion equations.  J. Hyperbolic Differ. Equ., 14(3), 415-454, 2017.
  * [27] U. Koley, D. Ray, and T. Sarkar.  Multi-level Monte Carlo finite difference methods for fractional conservation laws with random data.,  SIAM/ASA J. Uncertain. Quantif., 9(1), 65–105, 2021.
  * [28] S. N. Kruzhkov; First order quasilinear equations with several independent variables. Math. Sb. (N.S.) 81(123):228–255, 1970.
  * [29] P.L. Lions, B.Perthame, E. Tadmor, A kinetic formulation of multidimenional scalar conservation laws and related equations, J. Amer. Math. Soc. 7 (1) (1994) 169-191.
  * [30] G. Lv, J. Duan, H.Gao Stochastic nonlocal conservation laws on whole space, Computers $\&$ Mathematics with Applications Volume 77, Pages 1945-1962.
  * [31] G.Da Prato, and J. Zabczyk, (1992) Stochastic Equations in Infinite Dimensions.Encyclopedia of Mathematics and Its Applications 44. Cambridge Univ. Press, Cambridge. MR12
  * [32] D. Revuz, M. Yor,  Continuous Martingales and Brownian Motion, third ed., Grundlehren Math. Wiss. (FundamentalPrinciples of Mathematical Sciences), vol. 293, Springer-Verlag, Berlin, 1999
  * [33] L. Roncal, P. R. Stinga; Fractional Laplacian on the torus. Communications in Contemporary Mathematics Vol. 18, No. 03, 1550033 (2016)
  * [34] G. Vallet. and P. Wittbold; On a stochastic first-order hyperbolic equation in a bounded domain. Infin. Dimens. Anal. Quantum Probab. Relat. Top., 12(4), 613–651, 2009.


# On Simple Reactive Neural Networks for Behaviour-Based Reinforcement
Learning

Ameya Pore1 and Gerardo Aragon-Camarasa1 1 Computer Vision and Autonomous
group, School of Computing Science, University of Glasgow, G12 8QQ, UK. Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a behaviour-based reinforcement learning approach, inspired by
Brook’s subsumption architecture, in which simple fully connected networks are
trained as reactive behaviours. Our working assumption is that a pick and
place robotic task can be simplified by leveraging domain knowledge of a
robotics developer to decompose and train reactive behaviours; namely,
approach, grasp, and retract. Then the robot autonomously learns how to
combine reactive behaviours via an Actor-Critic architecture. We use an Actor-
Critic policy to determine the activation and inhibition mechanisms of the
reactive behaviours in a particular temporal sequence. We validate our
approach in a simulated robot environment where the task is about picking a
block and taking it to a target position while orienting the gripper from a
top grasp. The latter represents an extra degree-of-freedom of which current
end-to-end reinforcement learning approaches fail to generalise. Our findings
suggest that robotic learning can be more effective if each behaviour is
learnt in isolation and then combined them to accomplish the task. That is,
our approach learns the pick and place task in 8,000 episodes, which
represents a drastic reduction in the number of training episodes required by
an end-to-end approach ( 95,000 episodes) and existing state-of-the-art
algorithms.

##  I INTRODUCTION

Robots excel at using pre-programmed routines to perform repetitive tasks. A
significant barrier in their universal adoption beyond an enclosed environment
is their fragility and lack of robustness in complex environments. To tackle
these issues, recent advances in deep learning and deep Reinforcement Learning
(RL) have enabled the development of robotic solutions for complex and diverse
scenarios that have been intractable using classic traditional control
approaches. Examples include decision making for solving games [1, 2], and
continuous control tasks such as locomotion skills, dexterous manipulation and
grasping [3, 4, 5, 6]. However, a limitation to the widespread adoption of RL
algorithms in robotics is that RL approaches dramatically overfits the
idiosyncrasies of training environments [7, 8].

Figure 1: Schematic of the modular behaviour-based reinforcement learning
architecture. The goal of picking an object is subdivided into simpler
behaviours that are trained specifically for movements in x, y, z and $\theta$
of the end-effector. These behaviours are activated and inhibited by a
reactive network parametrised by an actor-critic policy. See Section III and
Fig. 2 for details.

In this paper, we depart from an end-to-end RL approach, and we investigate
whether it is possible to train a robot to pick up a block using manually
predefined behaviours that are then choreographed using RL. Specifically, we
propose a modular behaviour-based reinforcement learning architecture (see
Fig. 1) that is inspired by the subsumption architecture [9]. That is, we
start by training neural networks with a known solution, leveraging on domain
knowledge of a robotics developer. We, therefore, guide a robot to learn
specific low-level behaviours. Once these low-level behaviours are acquired,
an RL algorithm explores how to choreograph behaviours in different temporal
combinations; effectively learning to subsume behaviours that are not relevant
to the current state input. We validate our approach in a simulated
environment – MuJoCo simulator using the FetchPickandPlace environment [10,
11]. We must note that we go beyond the basic environment structure and allow
the block to spawn in random positions and orientations to include an
additional degree of freedom (gripper rotation) while grasping the block (Fig.
1(b)). Our contributions are therefore twofold: (i) an architecture that
learns isolated modular behaviours that which (ii) drastically reduces the
number of steps required to train a robot while performing a pick and place
task, i.e. picking and placing a block in a simulated environment. The source
code for this paper is available at https://github.com/cvas-ug/simple-
reactive-nn. A video summary of this paper can be found at
https://youtu.be/z7kUW9yyka0.

##  II RELATED WORK

Brooks proposed the Subsumption Architecture (SA) [9, 12] to mimic the
evolutionary path of intelligence. This architecture consists of designing
simple behaviours to achieve robust and complex behaviours in robots by
layering them in terms of complexity and execution time. Specifically, the SA
connects perception to action for robot control systems and coordinates
defined behaviours. In SA, complex behaviours subsume a set of simple
behaviours, and a task is accomplished by activating the appropriate behaviour
given an input state. Hierarchical Reinforcement Learning (HRL) [13, 14, 15]
resembles the SA in the sense that a complex task is automatically decomposed
into sub-task sequences, that are themselves built by machine-defined simple
actions. HRL learns and operates at different levels of temporal abstraction
by using multiple layers of policies that are trained to perform decision-
making and control at the successively higher level of behavioural and
temporal abstractions. The lowest-level policy of the hierarchy (subordinate
actions) applies actions to the environment, whereas the higher-level policies
are trained over a longer time scale.

Current HRL approaches include the work by Nachum et al. [13], where the
authors propose HIRO, a 2-level HRL approach that can learn off-policies.
While Levy et al. .[14] have built a hindsight experience with similar
hierarchical architecture as HIRO to increase the sample efficiency in sparse
reward conditions. The key distinction between these HRL approaches and our
approach is that in HRL, multiple policies are learnt in parallel and end-to-
end, whereas we attempt to learn them incrementally. Concretely, HRL
algorithms autonomously decide how to segment the main task into sub-tasks.
This segmentation is task-specific, and the decomposed sub-task, once trained,
would hardly be able to generalise to a different high-level task (Section I).
Moreover, learning multiple behaviours end-to-end leads to the curse of
dimensionality as the task becomes temporally elongated [15].

Our work is closely related to the architecture proposed by Konidaris et al.
[16] in behaviour-based reinforcement learning, where a topological map is
learnt to create task-relevant state spaces, and layered reinforcement
learning takes place over this map. Multiple learning models, multiple control
processes, and a complex environment result in complex learning behaviours.
However, a significant drawback of this approach is that it is not feasible to
build topological maps in many situations. We have addressed this by using
neural networks that learn feature representations from raw sensory data.
Similarly, the work by Frans et al. [17] presents an end-to-end method to use
shared policy primitives, within a distribution of tasks, and are switched
between by task-specific policies to execute over a large number of timesteps.
A master policy is learnt, and this policy selects a sub-policy to be active.
In this paper, we use shared parameters for state feature representation to
train a policy that selects the underlying trained primitive policies, and the
primitive policies are incrementally trained.

From the above, we can observe that state of the art approaches adopts an end-
to-end strategy to optimise and learn tasks and sub-tasks. However, humans
tend to learn simple behaviours first in order to compose complex behaviours
[18]. For example, while learning tennis, we start by learning basic
behaviours separately such as bouncing the ball, hitting the ball, serve, tp
name but a few. In contrast, an end-to-end approach would attempt to optimise
all possible behaviours mat the same time. That is, an artificially
intelligent agent requires millions of trials using sophisticated model-free
RL algorithms to complete simple tasks on simulations and games, whereas
humans learn behaviours in 50-100 attempts [19]. RL agents start solving each
problem tabula rasa with no human expert used as part of the training, whereas
we come in with a wealth of prior knowledge about the world, from physics to
semantics to affordances. In this paper, we propose that a robot learns basic
and simple behaviours. After building a set of these behaviours, a robot can
then learn how to choreograph these autonomously using reinforcement learning.

##  III MATERIALS & METHODS

Figure 2: The network takes as input a 28 state vector and outputs the $x$,
$y$, $z$ movement of the end-effector. The modules of approach (a), grasp (b)
and retract (c) are trained separately using behaviour cloning and combined to
accomplish the task. First, the state vector is given as an input to the
feature extraction layer. The extracted features are relayed on to the
reactive layers.

###  III-A Rationale

Our working assumption is that human knowledge can simplify the task of a
robot solving a pick and place task. Hence, we manually break up tasks into
simple behaviours and enable the robot to sequence them according to the state
encountered. The goal of the robot is to pick and take the block to a target
location. We, therefore, decompose the main pick and place task into
approaching, grasping a block, and retracting to a target point simple
behaviours (Fig. 2). These low-level behaviours start with general and
primitive abilities that are controlled, overridden, or subsumed by specific
goal-directed behaviours. We, therefore, adopt the SA architecture as the
means of defining behaviours to enable an incremental and sequential, bottom-
up operation of the system.

The overall behaviour of the robot is thus a consequence of the responses
within the environment. Behaviours rely on the state of the world without
maintaining a global internal representation [9]. We use a supervised learning
approach to learn from demonstrations instead of learning new policies from
scratch (Section III-B [20, 21]. That is, if the robot is approaching the
block, we explicitly teach it to move in the $x$, $y$ and $z$ Cartesian space.
We further repeat this for other behaviours. Once, these low-level behaviours
are learnt, we train an actor-critic RL architecture (Section III-C). The RL
architecture determines the activation and inhibition mechanisms of low-level
behaviours in a particular temporal sequence that ultimately give rise to the
high-level behaviour of picking and placing a block.

###  III-B Low-level Behaviours

Under the Subsumption Architecture, behaviours do not own memory and are
decomposed in layers, each with a predefined goal [9]. We, therefore, train a
neural network to learn a specific behaviour using demonstrations of low-level
behaviour. For this, we use behaviour cloning (BC), which learns a policy
through supervised learning in order to mimic the demonstrations [20, 21, 22].
Expert demonstrations of successful behaviours are used to train a network
which learns to imitate the expert providing these successful trajectories
[23].

We thus use a loss function computed on the demonstration examples as follows:

$\displaystyle L_{BC}=||\pi(s_{i}|\theta_{\pi})-a_{i}||^{2}$ (1)

where $a_{i}$ refers to the intended output of the behaviour. $\pi(s_{i})$
refers to the action predicted by the robot at state $s_{i}$ under the policy
parametrised by $\theta(\pi)$. We reparameterise the policy such that a sample
from $\pi_{\theta}(\cdot|s_{i})$ is drawn by computing a deterministic
function of the state, policy parameters and independent noise. We use a
neural network transformation $a_{i}=f(\phi_{i};s_{i})$, where $\phi_{i}$ is
an input noise vector, sampled from a fixed distribution, such as a spherical
Gaussian,

$\displaystyle
a_{i}=\tanh{(\mu_{\theta}(s_{i})+\sigma_{\theta}(s_{i})\phi_{i})},\phi_{i}\sim
N(0,1)$ (2)

The raw input of the kinematic coordinates and velocity of block and gripper
are fed to a feature extraction network consisting of two fully connected
layers with 128 neurons each. The output is mapped to six neurons in the last
layer that determines the mean and standard deviation of a Gaussian
distribution from which the movement values are sampled for the end-effector
in $x$, $y$ and $z$ for approach(a), grasping (b) and retract (c). For
training (Fig. 2), we start with the approach (a) behaviour module (Fig. 2).
That is, the output of the network is used to control the end-effector while
it is approaching the block (i.e $\mid
d_{eff}-d_{block}\mid<error$)111$d_{eff}$ refers to position of end-effector,
$d_{block}$ refers to the position of the block, and $error$ is manually set
to $0.01$ for training (a) and (c), and $0.005$ for training (b) and we use
hand-engineered solutions for grasping (b) and retract (c).

The BC loss function is backpropagated after each step while approaching. Once
the success rate reaches its maximum, training is stopped, and we save the
weights. Then, we train the grasp (b) behaviour module. For this, the weights
of the feature extraction (the first two layers) are frozen, and a similar
training process using BC is carried out where the output of the network is
used to control the end-effector while it is grasping the block. Finally, for
training the retract module, we froze the weights of (a) and (b) and a similar
training strategy is carried out until it reaches maximum success rate for
(c); in this case, we minimise the distance between the end-effector and the
target point (i.e. $\mid d_{target}-d_{block}\mid<0.01$).

TABLE I: Experimental strategies for training reactive behaviours before interfacing them with the high-level choreographer. # |  Name | Training Strategy
---|---|---
1 |  Sequential |  The approach behaviour (a) with a hand-engineered solution for grasping (b) and retract (c) behaviours. Once (a) is trained (i.e. reaches a desired performance level), (b) is trained using the output from the network (a) and a hand-engineered solution for (c). Similarly, once (b) is learnt, we train on (c), using outputs from the (a) and (b) networks.
2 |  Sequential + Freezing |  We start by training (a), with a hand-engineered solution for (b) and (c). After (a) is trained, we freeze the layers that extract features from the raw input. We then train (b), using the frozen weights for feature extraction and the output of (a), and a hand-engineered solution for (c). A similar approach is applied for scaling up and training on (c).
3 |  Separate |  We start by training (a) with a hand-engineered solution for (b) and (c). Once, (a) is trained, we do not use the output of the trained network (a). For training (b), we use a hand-engineered solution for (a) and (c). Similarly, training is carried out for (c). Once, each module is trained separately; we combine all the behaviours to accomplish the task.
4 |  Separate + Freezing |  We start by training (a) with a hand-engineered solution for (b) and (c). After (a) is trained, we freeze the layers that extract features from the raw input. For training (b), the output from the network (a) is not used. Instead, we use a hand-engineered solution for (a) and (c) with frozen weights of the feature extraction layer. Similarly, training is carried out for (c). Once, each module is trained separately; we combine all the behaviours to accomplish the task.
5 |  End-to-end |  In this case, we do not decompose simple behaviours. That is, state inputs are directly mapped to actions training using behaviour cloning for all time steps in the episodes. The action space consists of only the $x$, $y$ and $z$ coordinates of the end-effector without considering the end-effector orientation, $\theta$.

###  III-C High-level Choreographer

As stated in Section III-A, the high-level sequencer learns a policy that
choreographs a set of behaviours in order to solve a robotic pick and place
task. We consider the standard Markov decision process framework for picking
optimal behaviours to maximise rewards over discrete timesteps in an
environment $E$ [24]. At every timestep $t$, the robot is in a state $s$,
executes a behaviour $u_{t}$, receives a reward $r_{t}$, and $E$ evolves to
state $s_{t+1}$. Lets now denote the return by
$R_{t}=\sum_{i=t}^{T}{\gamma}^{(i-t)}r_{i}$, where $T$ is the horizon that the
robot optimises over, and ${\gamma}$ is a discount factor for future rewards.
The robot’s objective is to maximise the expected return from the start
distribution,

$\displaystyle J=E_{r_{t},s_{t}{\sim}E,u_{t}{\sim}\pi}[R_{0}].$ (3)

For the high-level choreographer, the extracted features from the block’s
position (Section III-B) are passed through an LSTM layer with $32$ units. Two
separate fully connected layers are used to predict the value function and the
activation/inhibition from the LSTM feature representation – see Fig. 1. The
aim of using a recurrent layer is that the agent should have a local memory of
the amount of task it has accomplished.In this paper, we adopted and tailored
the Asynchronous Actor-Critic architecture (A3C) [25] to serve as our higher
level choreographer. This architecture learns to sequence the low-level
behaviours described in Section III-B, and consists of a neural network called
the actor that predicts actions, and a network called critic that learns to
predict the value of a state-behaviour pair by optimising the Q-function. We
have used generalised advantage estimation to optimise the actor-critic model
[26].

##  IV EXPERIMENTS

Figure 3: OpenAI FetchPickandPlace environment as used in this paper.

In our experiments, we use a simulated environment based on the OpenAI
FetchPickandPlace environment [10, 11]. This environment is used as a
benchmark for testing algorithms for continuous control tasks such as robotic
manipulation and grasping. The goal in FetchPickandPlace environment is to
grasp a randomly positioned block and lift it to a target position. The
environment provides kinematic values of position, velocity and orientation of
the block and the gripper. In most studies related to RL algorithms
controlling the gripper in Fetch, the orientation of the block is fixed in
order to reduce the task complexity. We, however, activate the orientation of
the gripper to grasp different block orientations, as shown in Fig. 3.

Figure 4: Comparison between end-to-end and the 4 proposed training
strategies; i) Separate + Freezing (Yellow) and end-to-end (Blue); ii)
Separate (Blue) and end-to-end (yellow); iii) Sequential + Freezing (Blue) and
end-to-end (yellow); and iv) Sequential (blue) and end-to-end (yellow). Each
training strategy is run independently for three runs (Learning $n=3$) with
different seeds.

In order to investigate which training strategy is the most optimal within our
behaviour-based approach, we design five different training strategies, as
described in Table I. We must note that the sequence of the behaviours is
controlled manually in these training strategies, and behaviours are trained
following the approach in Section III-B. For the end-to-end training strategy,
we use a similar network structure for the feature extraction network as
described in Section III-B. The output is then mapped to 6 neurons in the last
layer that determines the mean and standard deviation of a Gaussian
distribution from which the movement values are sampled for the end-effector
in $x$, $y$ and $z$.

We then select the best performing training strategy and train the high-level
choreographer, as described in Section III-C. After training the high-level
choreographer, we evaluate two scenarios:

  1. 1.

Sparse reward condition: The robot receives a reward after it picks the block
to a target position.

  2. 2.

Dense reward condition: The robot receives a reward after each successful
completion of a behaviour. For example: if the robot selects approach at the
start of the episode, it receives a positive reward, on the contrary, if the
robot selects other behaviours, it receives negative rewards.

##  V RESULTS

Results are shown in Fig. 4. The first peak in Fig. 4(i-iv) denotes the
completion of training network (a) (approach) when the success rate reaches
100%, the second peak denotes the completion of training network (b)
(grasping), and the third peak denotes the completion for the network (c)
(retract). After each behaviour is trained, the success rate drops to zero
since learning the new behaviour starts from scratch. Step (d) denotes the
manual combination of (a), (b) and (c) to complete the task. In the end-to-end
case (i.e. training strategy 5 in Table I), actions are optimised for the
entire pick-and-place task.

We can observe that it takes a similar number of training episodes for the
end-to-end, Sequential + Freezing (Fig. 4-iii) and Sequential + No Freezing
(Fig. 4-iv) approaches. Separate + No Freezing (Fig. 4-ii) shows a slight
increase in the success rate, but its trend is close to the end-to-end
approach. However, Separate + Freezing (Fig. 4-i) shows a drastic increase in
success rate. In this case, the robot learns each skill separately with
freezing layers after initial training, and our behaviour-based approach can
reach 100% success rate in 6,000 episodes of training. This result suggests
that learning can be more effective if each skill is learnt in isolation and
then combined in order to learn the high-level task of pick and place.
Secondly, once the feature extraction network learns the latent feature space,
freezing the knowledge shows more learning potential for subsequent behaviours
(e.g. Fig. 4-i). That is, once the robot knows how to perceive the state
inputs, there is no advantage to learn the feature extraction layers again.

The best performing strategy is the Separate + Freezing (Fig. 4-i), and we,
therefore, use this strategy for the higher-level choreographer to learn how
to sequence these behaviours. The choreographer uses the features from the
feature extractor and selects the low-level actions based on the feedback of
the rewards provided, as described in Section III-C. This represents a simple
RL setting where an agent has to decide the actions to maximise the total
cumulative reward. Hence, we start training the high-level choreographer after
6,000 episodes (see Fig. 5). For this, we evaluate two RL agents that receive
rewards on different time-scales, namely dense and sparse reward settings
(Section IV. As expected, the RL agent that receives dense feedback can
accomplish the task faster; however, the RL agent that receives sparse rewards
performs close to the dense reward agent. One reason for this is that the pick
and place task considered in this paper is short-time horizon and requires
only three behaviours. We speculate that the robot with a dense reward
condition would outperform the sparse reward setting agent in a log-term
horizon task. This is left for future work as described in Section VI.

Figure 5: Comparison between the end-to-end approach (green) with the complete
behaviour-based RL architecture (i.e. Separate + Freezing combined with the
high-level choreographer – blue and orange).

By further inspecting Fig. 5, we can observe that the Separate + Freezing
combined with the high-level choreographer achieves 100% success accuracy in
approx. 8,000 episodes, effectively using only 2,000 episodes for training the
high-level choreographer. This result represents a drastic reduction in the
number of training episodes required by an end-to-end approach and the
existing state-of-the-art RL algorithms such as Deep Deterministic Policy
Gradients(DDPG) + Hindsight experience reply(HER) which takes 95,000 episodes
to learn the grasping task [11]. For our end-to-end approach, the robot is
able to reach a maximum of 60% success rate after 100,000 episodes. We
speculate that the reason for this is due to the network shallowness to learn
multiple behaviours at once, and backpropagation gets stuck in a local
minimum. However, the latter is outside the scope of this paper.

Our results suggest that training can be more effective if each module is
trained separately with other modules already trained; similar to what has
been found in curriculum learning [18]. Also, freezing the feature extraction
layer after initial training and using this layer for training other modules,
shows a considerable reduction in training time. This indicates that training
of complex learning systems should be accomplished in a structured fashion,
i.e. training simple modules first and independent of the rest of the network.
It aligns with the paradigm of bottom-up learning approach where complex
behaviours could arise by generating and combining simple ones [27] of which
motivates the adoption of Brook’s Subsumption Architecture [9] in this work.

##  VI CONCLUSIONS & FUTURE WORK

In this paper, we have used a bottom-up approach for training these modular
behaviours, i.e. Subsumption Architecture [9]. We have demonstrated, in
simulation, that long-time tasks can be decomposed and can be learnt
independently. The latter can give rise to behaviours that could accomplish a
variety of tasks. We train all the simple behaviours independently and then
combine them sequentially to complete the task. We argue that these decomposed
behaviours once trained could be used to accomplish different tasks and are
task agnostic. Further, the proposed behaviour-based RL architecture is a
simple feed forward neural network that maps the positional coordinates and
kinematic state input to low-level actions and high-level behaviours via a
learned and distributed internal representation.

From the results, we can state that our approach can learn to pick up a block
in approximately 8,000 episodes, as opposed to an end-to-end learning approach
and state-of-the-art RL approaches that take 95,000 episodes [11] on
simulation. The latter suggests that finding solutions using a model-free
approach is data inefficient and requires several trial and error iterations
for an RL agent to solve a task which is impractical in robotics. Our results
also suggest that by tapping into human knowledge to decompose simple
behaviours and separately learning these simple behaviours show a drastic
reduction in training time. For future work, we will deploy our approach in a
real robotic task to investigate how knowledge acquired by in simulation can
be generalised to other types of objects. For this, we will use deep learning
solutions for object recognition and for estimating the pose of objects, e.g.
[28], and investigate the use of continuous perception [29] to maintain
temporal consistency during task execution.

## ACKNOWLEDGEMENT

We thank Paul Siebert, Ali Al-Qallaf, Piotr Ozimek, Julio Caballero and John
Williamson for valuable discussions at the earlier stages of this research.
Ameya Pore thanks to Erasmus+international credit mobility programme and IISER
Pune. Gerardo Aragon-Camarasa acknowledges financial support from the EPSRC
Grant No. EP/S019472/1. We also acknowledge the support of NVIDIA Corporation
for the donation of the Titan Xp GPU used in this research.

## References

  * [1] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, “Playing atari with deep reinforcement learning,” _arXiv preprint arXiv:1312.5602_ , 2013.
  * [2] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot _et al._ , “Mastering the game of go with deep neural networks and tree search,” _nature_ , vol. 529, no. 7587, p. 484, 2016.
  * [3] T. Haarnoja, A. Zhou, S. Ha, J. Tan, G. Tucker, and S. Levine, “Learning to walk via deep reinforcement learning,” _arXiv preprint arXiv:1812.11103_ , 2018.
  * [4] S. Gu, E. Holly, T. Lillicrap, and S. Levine, “Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 3389–3396.
  * [5] T. Haarnoja, V. Pong, A. Zhou, M. Dalal, P. Abbeel, and S. Levine, “Composable deep reinforcement learning for robotic manipulation,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 6244–6251.
  * [6] M. Asenov, M. Burke, D. Angelov, T. Davchev, K. Subr, and S. Ramamoorthy, “Vid2param: Online system identification from video for robotics applications,” _CoRR_ , vol. abs/1907.06422, 2019. [Online]. Available: http://arxiv.org/abs/1907.06422
  * [7] C. Packer, K. Gao, J. Kos, P. Krähenbühl, V. Koltun, and D. Song, “Assessing generalization in deep reinforcement learning,” _arXiv preprint arXiv:1810.12282_ , 2018.
  * [8] A. Zhang, N. Ballas, and J. Pineau, “A dissection of overfitting and generalization in continuous reinforcement learning,” _arXiv preprint arXiv:1806.07937_ , 2018.
  * [9] R. A. Brooks, “Intelligence without representation,” _Artificial intelligence_ , vol. 47, no. 1-3, pp. 139–159, 1991.
  * [10] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.
  * [11] M. Plappert, M. Andrychowicz, A. Ray, B. McGrew, B. Baker, G. Powell, J. Schneider, J. Tobin, M. Chociej, P. Welinder _et al._ , “Multi-goal reinforcement learning: Challenging robotics environments and request for research,” _arXiv preprint arXiv:1802.09464_ , 2018.
  * [12] R. A. Brooks, “Elephants don’t play chess,” _Robotics and autonomous systems_ , vol. 6, no. 1-2, pp. 3–15, 1990.
  * [13] O. Nachum, S. S. Gu, H. Lee, and S. Levine, “Data-efficient hierarchical reinforcement learning,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 3303–3313.
  * [14] A. Levy, G. Konidaris, R. Platt, and K. Saenko, “Learning multi-level hierarchies with hindsight,” in _Proceedings of International Conference on Learning Representations_ , 2019.
  * [15] J. Kober, J. A. Bagnell, and J. Peters, “Reinforcement learning in robotics: A survey,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1238–1274, 2013.
  * [16] G. D. Konidaris and G. M. Hayes, “An architecture for behavior-based reinforcement learning,” _Adaptive Behavior_ , vol. 13, no. 1, pp. 5–32, 2005.
  * [17] K. Frans, J. Ho, X. Chen, P. Abbeel, and J. Schulman, “Meta learning shared hierarchies,” _arXiv preprint arXiv:1710.09767_ , 2017.
  * [18] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, “Curriculum learning,” in _Proceedings of the 26th annual international conference on machine learning_. ACM, 2009, pp. 41–48.
  * [19] R. Dubey, P. Agrawal, D. Pathak, T. L. Griffiths, and A. A. Efros, “Investigating human priors for playing video games,” _arXiv preprint arXiv:1802.10217_ , 2018.
  * [20] A. Nair, B. McGrew, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Overcoming exploration in reinforcement learning with demonstrations,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 6292–6299.
  * [21] M. Bojarski, D. Del Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang _et al._ , “End to end learning for self-driving cars,” _arXiv preprint arXiv:1604.07316_ , 2016.
  * [22] J. Nakanishi, J. Morimoto, G. Endo, G. Cheng, S. Schaal, and M. Kawato, “Learning from demonstration and adaptation of biped locomotion,” _Robotics and autonomous systems_ , vol. 47, no. 2-3, pp. 79–91, 2004.
  * [23] T. Hester, M. Vecerik, O. Pietquin, M. Lanctot, T. Schaul, B. Piot, D. Horgan, J. Quan, A. Sendonaris, I. Osband _et al._ , “Deep q-learning from demonstrations,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [24] R. S. Sutton, A. G. Barto _et al._ , _Introduction to reinforcement learning_. MIT press Cambridge, 1998, vol. 2, no. 4.
  * [25] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu, “Asynchronous methods for deep reinforcement learning,” in _International conference on machine learning_ , 2016, pp. 1928–1937.
  * [26] J. Schulman, P. Moritz, S. Levine, M. Jordan, and P. Abbeel, “High-dimensional continuous control using generalized advantage estimation,” _arXiv preprint arXiv:1506.02438_ , 2015.
  * [27] J. Gomes, S. M. Oliveira, and A. L. Christensen, “An approach to evolve and exploit repertoires of general robot behaviours,” _Swarm and Evolutionary Computation_ , vol. 43, pp. 265–283, 2018.
  * [28] G. Billings and M. Johnson-Roberson, “Silhonet: An rgb method for 6d object pose estimation,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3727–3734, Oct 2019.
  * [29] L. Martínez, J. R. del Solar, L. Sun, J. P. Siebert, and G. Aragon-Camarasa, “Continuous perception for deformable objects understanding,” _Robotics and Autonomous Systems_ , vol. 118, pp. 220 – 230, 2019.