††institutetext: School of Physical Sciences, Indian Association for the
Cultivation of Science, 2A and 2B Raja S.C. Mullick Road, Kolkata 700 032

# Revisiting big bang nucleosynthesis with a new particle species : effect of
co-annihilation with nucleons

Deep Ghosh<EMAIL_ADDRESS>
###### Abstract

In big bang nucleosynthesis (BBN), the light matter abundance is dictated by
the neutron-to-proton ($n/p$) ratio which is controlled by the standard weak
processes in the early universe. Here, we study the effect of an extra
particle species ($\chi$) which co-annihilates with neutron (proton), thereby
potentially changing the ($n/p$) ratio in addition to the former processes. We
find a novel interplay between the co-annihilation and the weak interaction in
deciding the ($n/p$) ratio and the yield of $\chi$. Large co-annihilation
strength ($G_{D}$) in comparison to the weak coupling ($G_{F}$), potentially
can alter the number of nucleons in the thermal bath modifying the ($n/p$)
ratio from its standard evolution. We find that the standard BBN prediction is
restored for $G_{D}/G_{F}\lesssim 10^{-2}$, while the mass of $\chi$ being
much smaller than the neutron mass. When the mass of $\chi$ is comparable to
the neutron mass, we can allow large $G_{D}/G_{F}~{}(\gtrsim 10^{2})$ values,
as the thermal abundance of $\chi$ becomes Boltzmann-suppressed. Therefore,
the ($n/p$) ratio is restored to its standard value via dominant weak
processes in later epochs. We also discuss the stability of the new particle
in an effective theory framework for co-annihilation. Further, the co-
annihilation interaction generates elastic scattering of $\chi$ and nucleons
at the next-to-leading order. This provides a way to probe the scenario in
direct detection experiments, if $\chi$ is accidentally stable over
cosmological timescale.

##  1 Introduction and Summary

Big bang nucleosynthesis (BBN) is one of the great achievements of both the
standard model (SM) of particle physics and hot big bang cosmology. The
observed primordial abundances of light matter in the universe agree with the
theoretical prediction of BBN to a good approximation. One of the important
quantities observed in this context is the primordial abundance of helium,
i.e. the mass fraction of helium, $Y_{P}=0.2449\pm 0.0040$ Aver:2015iza ,
which is sensitive to the neutron-to-proton ($n/p$) ratio at the initial epoch
of BBN. Neutrons and protons are in thermal equilibrium with the cosmic plasma
after the QCD phase transition via SM weak processes. Weak processes become
inefficient compared to the expansion of the universe at around temperature,
$T=1\,\rm MeV$. Subsequently the ($n/p$) ratio is frozen at a value around
($1/6-1/7$)Kolb:1990vq . This provides the initial condition for generating
the light nuclei in correct abundance at later epochs. It is apparent that any
alteration to this ratio due to some new physics phenomena would change the
prediction of BBN to a great extent. The addition of a new particle species
affects the ($n/p$) ratio broadly in two ways.

  1. I.

Additional particle species contributing significantly to the energy density
of the universe during BBN, changes the expansion rate of the universe which
in turn, delays or hastens the freeze-out of neutrons.

  2. II.

Chemical processes involving new particle species can alter the ($n/p$) ratio
by removing or adding extra nucleons to the thermal bath.

In scenario I, the inclusion of an extra relativistic species at the onset of
BBN is accounted as the effective number of extra neutrino species, defined as
the following.

$\displaystyle\Delta N_{eff}=\frac{8}{7}\frac{\rho_{\chi}}{\rho_{\gamma}}$ (1)

where, $\rho_{\chi}$ is the energy density of new particle species, $\chi$ and
$\rho_{\gamma}$ is the photon counterpart. To note, the above definition is
applicable untill the freeze-out of electron-positron annihilation, as the
relative heating of the photon bath has not yet taken place Weinberg:2008zzc .
Though the nuclear reactions associated with BBN happen for a range of the
photon temperature, for an order-of-magnitude estimate of the temperature and
the lower bound on the mass of $\chi$ (discussed in the next section), it is
sufficient to take the onset temperature for BBN to be, $T_{\gamma}\sim
1{~{}\rm MeV}$.

This parametrization works well for the species thermally decoupled from the
SM plasma before BBN. In ref.Kolb:1986nf ; Boehm:2013jpa  the authors have
shown that an additional particle species, strongly coupled either to photons
or to neutrinos via elastic scatterings during BBN alters neutrino-to-photon
temperature ratio ($T_{\nu}/T_{\gamma}$), thereby changing the standard weak
interaction rates. As a result the freeze-out time of neutrons gets affected,
in turn the final abundances of the light nuclei are altered.

In scenario II, the standard BBN (SBBN) prediction can potentially be altered
due to the infusion of extra nucleons either from the decay Kawasaki:2017bqm ;
Holtmann:1998gd ; Reno:1987qw ; Dimopoulos:1988ue ; Kusakabe:2008kf ;
Yeh:2024ors ; Jedamzik:2006xz  or from the annihilation Kolb:1981cx ;
Kawasaki:2015yya ; Depta:2019lbe ; Hisano:2009rc ; Ellis:2011sv  of new
particles. In most cases, the decay or annihilation products either inject
electromagnetic energy into the cosmic soup via energetic electrons, photons
or introduce additional inelastic scatterings between nucleons and nascent
antinucleons and other hadrons. The newly formed light nuclei can undergo
photo-dissociation due to electromagnetic energy injection, whereas
scatterings of hadrons can directly alter the $(n/p)$ ratio, thereby modifying
the light nuclei abundances at later epochs Nollett:2013pwa ; Kawasaki:1994sc
; Protheroe:1994dt ; Kusakabe:2014ola ; AlbornozVasquez:2012emy . Besides, the
co-annihilation Griest:1990kh  of dark sector (with more than one species)
particles can give rise to similar effects in light nuclei abundances
Jittoh:2008eq ; Jittoh:2010wh ; DAgnolo:2018wcn ; Ellis:2011sv .

In this article we have introduced a new particle species $\chi$, which co-
annihilates with nucleons. In the co-annihilation process detailed in Sec.2, a
neutron (proton) and a $\chi$ ($\bar{\chi}$) particle are removed from the
thermal bath. Therefore, the ($n/p$) ratio and the yield of $\chi$
($\bar{\chi}$) can simultaneously be affected by the co-annihilation at the
relevant epoch. The freeze-out value of the ($n/p$) ratio and the abundance of
$\chi$ are decided by the relative strength of the standard weak processes and
the newly added co-annihilation, which we have discussed in great details in
Sec.3. In addition, we also discuss the stability of $\chi$ in the context of
the effective theory for co-annihilation. In sec.4 we discuss the elastic
scattering of nucleon (electron) and $\chi$ generated radiatively (detailed in
Appendix A) from the co-annihilation interaction, leading to a possible probe
of accidentally stable $\chi$ in direct detection experiments.

We now summarize our findings regarding the modification of the SBBN scenario
due to newly added co-annihilation process. The effect of co-annihilation
depends on the initial relative abundances of the neutron (proton) and $\chi$
($\bar{\chi}$) in the cosmic soup. The number density of neutron (proton) is
decided by the observed baryon asymmetry, whereas the number density of $\chi$
is assumed to be thermal. Hence, the number density is decided by the mass of
$\chi$ ($m_{\chi}$) and the temperature of the cosmic soup. For
$m_{\chi}<\mathcal{O}(m_{n})$, $m_{n}$ being the neutron mass, $\chi$ freezes
out relativistically, therefore the ambient number density becomes too large
compared to the number density of nucleons. Consequently, the co-annihilation
processes ($n+\chi\rightarrow p+e^{-}$ and $p+\bar{\chi}\rightarrow n+e^{+}$)
keep the ($n/p$) ratio fixed at its relativistic equilibrium value, i.e.
$R_{0}=1$ even at late times altering its SBBN value. This puts a constraint
on the co-annihilation strength, i.e. $G_{D}/G_{F}\lesssim 10^{-2}$, $G_{F}$
being the weak interaction.

The constraint on the $G_{D}/G_{F}$ is significantly relaxed for
$m_{\chi}\sim\mathcal{O}(m_{n})$, unlike the previous case. This is due to an
interesting interplay between the weak interaction and the co-annihilation in
deciding the evolution of the ($n/p$) ratio and the yield of $\chi$
($\bar{\chi}$) over different epochs. Initially ($T>>1\,\rm MeV$) there is a
large number of $\chi$ ($\bar{\chi}$) and for large $G_{D}/G_{F}$ ($\gtrsim
10^{2}$) values the co-annihilation dominates over weak processes. However,
the combined co-annihilation reactions of $\chi$ and $\bar{\chi}$ particles
keep the ($n/p$) ratio at its SBBN value, $R_{0}=1$. Soon for such a large co-
annihilation strength the number density of new particle becomes Boltzmann-
suppressed. Consequently, the co-annihilation rate per neutron (proton)
becomes sub-dominant and the SBBN prediction of the ($n/p$) ratio is
eventually satisfied via weak processes in the later epoch. We notice that for
$m_{\chi}\sim\mathcal{O}(m_{n})$, we can allow large values of $G_{D}/G_{F}$
without altering SBBN predictions.

In this co-annihilating scenario, $\chi$ can decay to SM states in absence of
any symmetry protection, potentially affecting BBN predictions. In the
effective model of co-annihilation discussed in the next section, the tree-
level decay of $\chi$ is forbidden on the kinematic ground. It might decay via
some loop-induced processes, although the lifetime of $\chi$ is larger than
the duration of BBN, preserving the main outcome of the previous analysis.
However, if the particle is accidentally stable, we can probe the BBN-allowed
parameter space with $m_{\chi}\sim\mathcal{O}(m_{n})$ and $G_{D}/G_{F}\gtrsim
10^{2}$ in future direct detection experiments.

##  2 Co-annihilation with nucleons

In the SBBN, the initial number densities of neutron and proton are controlled
by the following electro-weak processes in which scattering processes freeze
out at $T\sim 1\,\rm MeV$ and the $(n/p)$ ratio becomes fixed. The frozen-out
value of the $(n/p)$ ratio changes slightly due to occasional decays of
neutron until the helium atom is produced, in which most of the neutrons are
trapped.

A. $n+\nu_{e}\rightleftharpoons p+e^{-}$, B. $n+e^{+}\rightleftharpoons
p+\bar{\nu_{e}}$, C. $n\rightleftharpoons p+e^{-}+\bar{\nu_{e}}$

We alter the standard scenario by incorporating a co-annihilation process,
i.e. $n+\chi\rightleftharpoons p+e^{-}$, where $\chi$ is considered as a
charge-neutral Dirac fermion. Such a process can be motivated by the lepto-
quark mediator models Belfatto:2021ats  in solving B-physics anomalies as well
as solving the neutron decay anomaly Strumia:2021ybk . Moreover, the process
considered here can be instrumental for detecting sub-GeV particles (dark
matter or exotic neutrinos) on the beta-decaying nuclear targets Dror:2019dib
; Long:2014zva ; Cocco:2007za . Apart from being a fundamental particle,
$\chi$ can well be a composite particle like neutron, which have been studied
in Refs.Cline:2021itd ; McKeen:2020oyr ; Carenza:2022pjd  and references
therein. Here, we want to investigate the cosmological implication of the
process being agnostic about the exact ultra-violet completion of the low
energy effective theory, taken as

$\displaystyle\mathcal{L}\supset\frac{G_{D}}{\sqrt{2}}\left(\bar{\psi}_{p}\Gamma^{a}\psi_{n}\right)\,\left(\bar{\psi}_{e}\Gamma_{a}\psi_{\chi}\right)$
(2)

where $\Gamma_{a}$ in general, represents all possible independent
combinations of Dirac matrices. As a proto-type example of our scenario, We
assume the standard $V-A$ structure for $\Gamma^{a}$, which is reminiscent of
SM weak interaction. $\psi_{i}$ is the fermion field for the $i^{th}$
particle. We want to study only the effect of the co-annihilation on the SBBN
case, therefore we take $m_{\chi}<m_{n}+m_{p}+m_{e}\sim 1.88\,\rm GeV$ to
prevent the tree-level decay i.e. $\chi\rightarrow\bar{n}+p+e^{-}$ on the
kinematic ground.111There can be loop-induced decay modes of $\chi$ allowed by
Eq.2, which have been discussed in Sec.3.2. In addition, to exclude the
scenario I discussed in Sec.1, $\chi$ should not contribute significantly to
the energy density of the universe during BBN. If $\chi$ is a decoupled but
internally thermalized relativistic species, the measurement of $N_{eff}$
Fields:2019pfx  sets an upper bound on the temperature ($T_{\chi}$) of $\chi$
using Eq.1 as,

$\displaystyle\left(T_{\chi}/T_{\gamma}\right)\lesssim 0.6\hskip
28.45274pt\text{at 68\% C.L.}$ (3)

If $\chi$ is thermalised with the photon bath during BBN, the lower bound on
$m_{\chi}$ comes to be approximately $1.68\,\rm MeV$. For more detailed
analysis on this frontier see Ref.An:2022sva , and references therein. In the
subsequent analysis, we take the number density of $\chi$ to be thermal to
start with. Thus we get a very general mass window for a stable thermalised
extra fermion species as,

$\displaystyle 1.68\,{\rm MeV}\lesssim m_{\chi}\leq 1.88\,\rm{GeV}$ (4)

Now, the co-annihilation process is active once neutrons and protons are
available in the bath after the QCD phase transition, $T\sim 150\,\rm MeV$. To
check the thermalization condition for $\chi$ via the co-annihilation with
neutrons in the expanding universe, we need $\Gamma(T)>H(T)$. Here,
$\Gamma(T)$ is the co-annihilation rate per $\chi$ particle and $H(T)$ is the
Hubble constant at temperature $T$. The interaction rate depends on the co-
annihilation cross-section and the number density of neutrons at a certain
epoch. The number density of neutrons is estimated from the observed baryon
asymmetry, i.e.

$\displaystyle\Delta_{B}=\frac{n_{B}-n_{\bar{B}}}{s}=\frac{n_{n}+n_{p}}{s}=Y_{n}+Y_{p}\simeq
10^{-10}$ (5)

Where, $Y_{i}=n_{i}/s$, $n_{i}$ is the number density of $i^{th}$ particle and
$s$ is the entropy density of the universe. For $T>>1\,\rm MeV$ the co-moving
number density of neutron can well be taken as $Y_{n}=Y_{p}=\Delta_{B}/2$. The
s-wave contribution to the thermally averaged co-annihilation cross-section is
given by,

$\displaystyle\expectationvalue{\sigma
v}\cong\frac{G^{2}_{D}}{2\pi}\left(m^{2}_{\chi}+2m_{\chi}m_{n}\right)$ (6)

in which we have assumed electron to be massless and neutron and proton to
have identical masses. Assuming usual radiation-dominated universe for the
relevant epoch, we arrive at an approximate condition for chemical equilibrium
of $\chi$ in the non-relativistic regime as the following.

$\displaystyle G_{D}\gtrsim 4\times 10^{-4}\,{\rm GeV^{-2}}\left(\frac{1\,{\rm
GeV^{2}}}{m^{2}_{\chi}+2m_{\chi}m_{n}}\right)^{1/2}\left(\frac{0.1{\rm
GeV}}{T}\right)^{1/2}\left(\frac{10}{g_{*}}\right)^{1/2}$ (7)

Thus, the thermal number density of $\chi$ is ensured around $T\sim 100~{}{\rm
MeV}(\sim m_{n}/10)$, as shown in Figs.1, 2. The initial production of $\chi$
should happen anytime before $T\sim 100~{}\rm MeV$ either from the SM thermal
bath or via the decay of the inflaton field. In addition, $\chi$ can have a
non-zero chemical potential related to the lepton number similar to electrons.
Therefore, it can take part in leptogenesisFukugita:1986hr  and subsequent
baryogenesis via the sphaleron effect Rubakov:1996vz , if its production
happens at sufficiently high temperature, i.e. $T\gtrsim 100~{}{\rm GeV}$. For
simplicity, we assume vanishing chemical potential of $\chi$ throughout the
cosmic history.

To study the neutron freeze-out in this modified scenario we need to write
Boltzmann equations for neutrons, protons and $\chi$. Due to high entropy
density per baryon the light nuclei are not synthesized immediately after the
neutron freeze-out, e.g. the synthesis of $He^{4}$ takes place around
$T=0.1\,\rm MeV$. Hence, for $T\gtrsim 0.1\,\rm MeV$ the relevant processes
are only the inter-conversions between neutrons and protons via three weak
processes (A,B,C) and newly added co-annihilation. As we are considering a
thermal number density of the new particle species with vanishing chemical
potential, there are two relevant co-annihilation processes involving $\chi$
and $\bar{\chi}$, i.e.

$\displaystyle 1.~{}n+\chi\rightarrow
p+e^{-}~{},~{}~{}2.~{}p+\bar{\chi}\rightarrow n+e^{+}$ (8)

The total co-moving number density of baryons is conserved within the relevant
epoch, which implies

$\displaystyle\frac{dY_{p}}{dx}+\frac{dY_{n}}{dx}=0\,\hskip
28.45274pt\text{for $T\gtrsim 0.1\,\rm MeV$}$ (9)

where $x=m_{n}/T$. Now, defining the ($n/p$) ratio as $R=Y_{n}/Y_{p}$ and
$Y_{\chi}$ ($Y_{\bar{\chi}}$) as the co-moving number density of $\chi$
($\bar{\chi}$), we write Boltzmann equations using Eq.9 as the following.

$\displaystyle\frac{dR}{dx}$
$\displaystyle=\frac{1+R}{Hx}\bigg{[}\lambda(p\rightarrow
n)-\lambda(n\rightarrow
p)R\bigg{]}-\frac{(1+R)s}{Hx}\bigg{[}R\expectationvalue{\sigma
v}_{1}\left(Y_{\chi}-\frac{R_{0}}{R}\,Y^{0}_{\chi}\right)-\expectationvalue{\sigma
v}_{2}\left(Y_{\bar{\chi}}-\frac{R}{R_{0}}\,Y^{0}_{\chi}\right)\bigg{]}$
$\displaystyle\frac{dY_{\chi}}{dx}$
$\displaystyle=-\frac{s\expectationvalue{\sigma
v}_{1}}{Hx}\frac{\Delta_{B}R}{1+R}\bigg{[}\left(Y_{\chi}-\frac{R_{0}}{R}\,Y^{0}_{\chi}\right)\bigg{]}$
$\displaystyle\frac{dY_{\bar{\chi}}}{dx}$
$\displaystyle=-\frac{s\expectationvalue{\sigma
v}_{2}}{Hx}\frac{\Delta_{B}}{1+R}\bigg{[}\left(Y_{\bar{\chi}}-\frac{R}{R_{0}}\,Y^{0}_{\chi}\right)\bigg{]}$
(10)

where $Y^{0}_{\chi}$ is the equilibrium co-moving number density for $\chi$,
$R_{0}=Y^{0}_{n}/Y^{0}_{p}\simeq e^{-Q/T}$ and $Q=m_{n}-m_{p}=1.293\,\rm MeV$.
$\lambda(n\rightarrow p)$ is the neutron to proton conversion rate via SM weak
processes defined assuming zero chemical potential for neutrinos as
Wagoner:1966pv ,

$\displaystyle\lambda(n\rightarrow
p)=\frac{1}{\tau}\int^{\infty}_{1}d\epsilon\bigg{[}\frac{\sqrt{\epsilon^{2}-1}~{}\epsilon~{}(q+\epsilon)^{2}}{\left(\exp\left(\frac{m_{e}\epsilon}{T}\right)+1\right)\left(\exp\left(-\frac{(q+\epsilon)m_{e}}{T_{\nu}}\right)+1\right)}+\frac{\sqrt{\epsilon^{2}-1}~{}\epsilon~{}(\epsilon-q)^{2}}{\left(\exp\left(-\frac{m_{e}\epsilon}{T}\right)+1\right)\left(\exp\left(\frac{m_{e}(\epsilon-q)}{T_{\nu}}\right)+1\right)}\bigg{]}$
(11)

where, $\tau$ is the neutron lifetime and $T_{\nu}$ is the neutrino
temperature and $q=(m_{n}-m_{p})/m_{e}=2.53$. The proton to neutron conversion
rate, $\lambda(p\rightarrow n)$ is achieved replacing $q$ by $-q$ in the
expression of $\lambda(n\rightarrow p)$. In Eq.11 $T_{\nu}$ can be written as
a function of photon temperature, $T_{\gamma}=T$ and $T$ is replaced by the
scaled temperature variable $x$. Further, $\expectationvalue{\sigma v}_{i}$ is
thermally averaged cross-section of process $i$ defined in Eq.8. The s-wave
contribution to these thermally averaged cross-section can be read off from
Eq.6 as,

$\displaystyle\expectationvalue{\sigma
v}_{1}\cong\frac{G^{2}_{D}}{2\pi}(m^{2}_{\chi}+2m_{\chi}m_{n})~{}~{};~{}~{}\expectationvalue{\sigma
v}_{2}\cong\frac{G^{2}_{D}}{2\pi}(m^{2}_{\chi}+2m_{\chi}m_{p})$ (12)

It is apparent from above Boltzmann equations that both the ($n/p$) ratio and
the number density of $\chi$ are inter-related due to added co-annihilation.
In particular, $Y_{\chi}$ ($Y_{\bar{\chi}}$) follows its equilibrium form
($Y^{0}_{\chi}$) as long as $R=R_{0}$, i.e. the ($n/p$) ratio maintains its
equilibrium value.

##  3 Results

###  3.1 Study of the evolution equations

We shall now study the evolution of $R$ and $Y_{\chi}~{}(Y_{\bar{\chi}})$ by
solving Eq.10 numerically keeping $G_{D}$ and $m_{\chi}$ as free parameters.
The evolution equations for our scenario are solved upto $x=5000$, which
corresponds to $T\simeq 0.2\,\rm MeV$. As discussed earlier, around
$T=0.1\,\rm MeV$ other nucleons start building up significantly, prompting to
include all nuclear reactions into our Boltzmann equations. The temperature
range of the current study is adequate to capture general features of co-
annihilation.

It is apparent that in a co-annihilation-type process, the reaction rate per
particle is different for each of the colliding particles. Therefore, the co-
annihilation process is very much sensitive to the initial relative abundances
of two particle species, which is decided by $m_{\chi}$. This is because the
number density of $\chi$ is thermal and the neutron number density is set by
the baryon asymmetry. To demonstrate the effect of new physics on BBN we take
two benchmark points (BPs) for $m_{\chi}$, i.e. $m_{\chi}=50\,\rm MeV$ and
$m_{\chi}=1.0\,\rm GeV$.

Figure 1: Evolution of the ($n/p$) ratio ($R$) (Left panel) and the co-moving
number density ($Y_{\chi}$ $(Y_{\bar{\chi}})$) of the additional species
(Right panel) shown for $m_{\chi}=50\,\rm MeV$, varying the scaled co-
annihilation strength, $G_{D}/G_{F}$.

BP - I : Relativistic freeze-out of $\chi$

In Fig.1 we have shown the evolution of $R$ (left panel) and $Y_{\chi}$ (right
panel) for $m_{\chi}=50\,\rm MeV$ for different values $G_{D}$, scaled by
$G_{F}=1.166\times 10^{-5}\,\rm GeV^{-2}$. In the right panel, the black
dashed line denotes the equilibrium number density ($Y^{0}_{\chi}$) of $\chi$
and other lines correspond to the yields with different $G_{D}/G_{F}$ values.
We find that the additional species, $\chi$ freezes out approximately at its
equilibrium value, almost independent of the variation of $G_{D}/G_{F}$. For
instance, $G_{D}/G_{F}$ varying from $10^{-2}$ to $10^{3}$, the freeze-out
temperature is approximately same, $T_{F}\simeq 23\,\rm MeV$ which in terms of
the scaled temperature becomes $m_{\chi}/T_{F}\simeq 2.2$. This indicates the
relativistic freeze-out of $\chi$ as $m_{\chi}/T_{F}\lesssim 3$ Kolb:1990vq .

In the left panel, the black dashed line denotes the equilibrium value of the
($n/p$) ratio ($R_{0}$) and the cyan dashed line corresponds to that of the
SBBN scenario, i.e. $G_{D}/G_{F}=0$. As long as the co-annihilation keeps
$\chi$ in chemical equilibrium, the $(n/p)$ ratio follows the equilibrium
value ($R_{0}$) because there is no extra loss or gain of neutrons via the co-
annihilation. After the freeze-out of $\chi$, the ($n/p$) ratio departs from
its standard evolution for $G_{D}/G_{F}=\\{10,10^{2},10^{3}\\}$, whereas for
$G_{D}/G_{F}=10^{-2}$ (red solid line in the left panel) it is unaltered from
the SBBN case. This feature can be understood from a comparison between the
weak interaction rate ($\Gamma_{w}$) and the co-annihilation rate
($\Gamma_{c}$) per neutron as the following.

$\displaystyle\frac{\Gamma_{c}}{\Gamma_{w}}\simeq\left(\frac{G_{D}}{G_{F}}\right)^{2}$
(13)

From Eq.13 we note that for $G_{D}/G_{F}=10^{-2}$ the weak interaction
dominates over the co-annihilation throughout the evolution history of neutron
and $\chi$. For $G_{D}/G_{F}=\\{10,10^{2},10^{3}\\}$ the ($n/p$) ratio remains
at its relativistic value due to dominant co-annihilation processes. In
particular, the processes, $\chi+n\rightarrow p+e^{-}$,
$\bar{\chi}+p\rightarrow n+e^{+}$ keep the ($n/p$) ratio around unity even
after the chemical equilibrium of the co-annihilation is broken. This is due
to a large number of $\chi$ particles present in the cosmic soup compared to
neutrons (protons). In fact, in this case $Y_{\chi}\simeq 10^{-2}$ and
$Y_{p}=Y_{n}\simeq 10^{-10}$, which imply that there are $10^{8}$ number of
$\chi$ particles available in $1$ neutron (proton) in the bath. Therefore, the
occasional co-annihilation is sufficient to disrupt drastically the neutron
freeze-out process as in SBBN. To sum up, for the relativistic freeze-out
case, the co-annihilation with larger interaction strength than the weak
interaction strength enhances the ($n/p$) ratio by one order of magnitude from
the SBBN value at temperature, $T\simeq 10\,\rm MeV$ (see the left panel of
Fig.1). This eventually jeopardizes the predictions of light nuclei abundance,
albeit the weak processes are still on.

Figure 2: Evolution of the ($n/p$) ratio ($R$) (Left panel) and the co-moving
number density ($Y_{\chi}$) of the additional species (Right panel) have been
shown for $m_{\chi}=1.0\,\rm GeV$,varying the scaled co-annihilation strength,
$G_{D}/G_{F}$. The dashed line denotes the co-moving number density of
$\bar{\chi}$ in the right panel.

BP - II : Non-relativistic freeze-out of $\chi$

Similar to the previous one, we now show the evolution of $R$ (left panel) and
$Y_{\chi}$ ($Y_{\bar{\chi}}$) (right panel) in Fig.2 for $m_{\chi}=1.0\,\rm
GeV$ for different values of $G_{D}/G_{F}$. Unlike the previous case, the
freeze-out of $\chi$ depends on the co-annihilation strength. For example, the
freeze-out temperature of $\chi$ for $G_{D}/G_{F}=10^{2}$ (red solid line) is
$T_{F}\simeq 50\,\rm MeV$, whereas for $G_{D}/G_{F}=200$ (magenta dashed line
for the corresponding yield) the freeze-out happens bit late, i.e. at
$T_{F}\simeq 20\,\rm MeV$. To note, for $m_{\chi}=1\,\rm GeV$ the freeze out
of $\chi$ happens in the non-relativistic regime, as $m_{\chi}/T_{F}>>3$ in
our example scenarios. In particular, there is a non-trivial relation between
freeze-out temperature and the co-annihilation strength, i.e. $T_{F}\propto
G^{-2}_{D}$ evident from Eq.7. This should be contrasted with the freeze-out
temperature in the standard WIMP scenario of which the dependence on the
coupling strength is rather weak, i.e. a logarithmic dependence
Rubakov:2017xzr . The difference manifests from the fact that the freeze-out
condition is determined by the baryon asymmetry in our case, whereas in the
WIMP scenario it is decided by the thermal densities with vanishing chemical
potentials. For $G_{D}/G_{F}=\\{500,10^{3}\\}$ (blue dotted and orange dot-
dashed line respectively) the co-annihilation keeps $Y_{\chi}$ in its
equilibrium form ($Y^{0}_{\chi}$, black dashed line) for the relevant epoch.

In the left panel, We notice that the SBBN scenario for the ($n/p$) ratio
($G_{D}/G_{F}=0$ (cyan dashed line)) is reproduced for largish value of
$G_{D}/G_{F}$ values, i.e. $G_{D}/G_{F}\geq 10^{2}$. This is a sheer contrast
to the relativistic case where largish value of $G_{D}/G_{F}$ destroys the
SBBN prediction by changing the ($n/p$) ratio significantly. This can be
understood by a novel interplay between the weak interaction and the co-
annihilation throughout the evolution history of neutron(proton) and $\chi$.
To understand the situation, let’s look at the relative strengths of the co-
annihilation and the weak interaction rate per neutron which is given by the
following expression, noticeably different from Eq.13.

$\displaystyle\frac{\Gamma_{c}}{\Gamma_{w}}\simeq\left(\frac{G_{D}}{G_{F}}\right)^{2}\left(\frac{m^{2}_{\chi}+2m_{\chi}m_{n}}{T^{2}}\right)\left(\frac{m_{\chi}}{T}\right)^{3/2}e^{-m_{\chi}/T}$
(14)

At around $T\simeq 90\,\rm MeV$, for $G_{D}/G_{F}=10^{2}$,
$\Gamma_{c}/\Gamma_{w}\simeq 10^{3}$, i.e. the co-annihilation rates per
neutron (proton) are large by several orders of magnitude from the standard
weak processes. Therefore, at this epoch, $\chi+n\rightarrow p+e^{-}$,
$\bar{\chi}+p\rightarrow n+e^{+}$ decide the ($n/p$) ratio. These processes
together do not the alter the the ($n/p$) ratio. As a result the SBBN scenario
is restored initially, i.e. $R_{0}=1$. Very soon, the weak interaction becomes
dominant as the number density of $\chi$ is Boltzmann-suppressed. The dominant
co-annihilations remove $\chi$ particles substantially from the thermal bath,
contrary to the relativistic scenario. Therefore, the available $\chi$
($\bar{\chi}$) particles per neutron (proton) become scarce to facilitate co-
annihilation any further in the later epoch. In particular, at $T\sim 40\,\rm
MeV$ for $G_{D}/G_{F}=10^{2}$, $\Gamma_{c}/\Gamma_{w}\simeq 10^{-2}$.

To note, throughout the evolution history, the standard evolution of the
($n/p$) ratio is maintained by two different kinds of reaction. At the early
epoch, the co-annihilation processes keep neutrons and protons in thermal
equilibrium, subsequently dominant weak processes control the evolution of the
same. Hence, in the non-relativistic scenario, largish values of $G_{D}/G_{F}$
are allowed from the BBN constraints. In passing, we also note that in the
non-relativistic regime, the freeze-out number densities of $\chi$ and
$\bar{\chi}$ for $G_{D}/G_{F}=\\{100,200\\}$ are slightly different. This is
due to the fact that the interaction rate of $\chi$ with neutrons is smaller
than the interaction rate of $\bar{\chi}$ with protons. In the non-
relativistic regime, the number density of neutron is smaller than that of
protons due to the mass difference. Now, for $\expectationvalue{\sigma
v}_{1}\simeq\expectationvalue{\sigma v}_{2}$ (see Eq.12), the surviving number
of $\chi$ particles is larger than $\bar{\chi}$ particles at the freeze-out.
Therefore, a small asymmetry is manifested in number densities of $\chi$ and
$\bar{\chi}$ particles.

We can now summarize our discussion regarding the interplay between the weak
interaction and the co-annihilation for both the cases, i.e. the relativistic
and the non-relativistic freeze-out of $\chi$ using an instructive diagram
shown in the left panel of Fig.3. The freeze-out value of the ($n/p$) ratio is
denoted by $R_{F}$, which is calculated at $m_{n}/T_{F}=5000$ varying
$G_{D}/G_{F}$ continuously over several orders of magnitude for different
$m_{\chi}$ values. In SBBN scenario, $R_{F}(\rm BBN)\approx 1/7$
Bernstein:1988ad ; Mukhanov:2003xs , denoted by the black dashed line in the
left panel including the effect of the neutron decay. As suggested by the
previous two example scenarios, the mass of $\chi$ controls two different
features in determining $R_{F}$.

Figure 3: Left panel : The freeze-out value of the ($n/p$) ratio, $R_{F}$
shown as a function of the scaled co-annihilation strength, $G_{D}/G_{F}$ for
different values of $m_{\chi}$. Right panel : Allowed parameter space of
$m_{\chi}$ and $G_{D}/G_{F}$ shown by the shaded area.

  1. I.

For $m_{\chi}<<\mathcal{O}(m_{n})$, there is a large number of $\chi$ compared
to neutrons, which makes weak processes inefficient for most of the time. Thus
$R_{F}$ deviates from its SBBN value for $G_{D}/G_{F}\gtrsim 10^{-2}$. This is
illustrated for two masses, i.e. $m_{\chi}=100\,\rm MeV$ (magenta dot-dashed
line) and $m_{\chi}=50\,\rm MeV$ (orange solid line) in the left panel of
Fig.3. We note, for $G_{D}/G_{F}\approx 10^{2}$, $R_{F}\approx 1$, which
completely ruins the SBBN predictions.

  2. II.

For $m_{\chi}\sim\mathcal{O}(m_{n})$ or larger, $R_{F}$ remains at its SBBN
value for two regions, i.e. for $G_{D}/G_{F}\lesssim 10^{-1}$ and for
$G_{D}/G_{F}\gtrsim 10^{2}$. The SBBN prediction is altered only for the
intermediate region, i.e. $10^{-1}\lesssim G_{D}/G_{F}\lesssim 10^{2}$. The
situation is depicted through two illustrative masses, i.e. $m_{\chi}=1\,\rm
GeV$ (red dashed line) and $m_{\chi}=1.5\,\rm GeV$ (blue dotted line). For
$G_{D}/G_{F}\lesssim 10^{-1}$, the weak interaction dominates over the co-
annihilation throughout the evolution history, therefore the SBBN scenario is
trivially satisfied. For $G_{D}/G_{F}\gtrsim 10^{2}$, the relative abundance
of $\chi$ to neutrons (protons) ($Y_{\chi}/Y_{n/p}$) becomes too small to
modify $R_{F}(\rm BBN)$ via the interplay between the weak and the co-
annihilation processes discussed earlier. In the intermediate region, the co-
annihilation can not keep $Y_{\chi}$ in its equilibrium as indicated in Eq.7.
Therefore, $Y_{\chi}/Y_{n}$ becomes large enough to ruin BBN predictions.

To sum up ( as depicted by the shaded area in the right panel of Fig.3), when
$m_{\chi}<<\mathcal{O}(m_{n})$, the BBN prediction allows only small values of
$G_{D}/G_{F}$. The parameter space for $G_{D}/G_{F}$ is enhanced for
$m_{\chi}\sim\mathcal{O}(m_{n})$ or larger, as both small and large values of
$G_{D}/G_{F}$ can reproduce the BBN predictions. We also show the two mass
bounds of $\chi$ (as in Eq.4) by the dot-dashed line in the right panel
figure.

###  3.2 Probable Loop-induced decays of $\chi$

The result shown in the Fig.3 is modified if $\chi$ decays during the BBN
dynamics. Typically, a particle species decaying to electron-positrons,
photons or neutrinos modify the hubble expansion rate thus modifying the
primordial abundances. In addition, it can change the neutrino-to-photon
temperature ratio via energy injection to either of two sectors after the
freeze-out of electron-positron pair annihilation. For detailed analysis
regarding bounds on the lifetime and the abundance of a decaying particle
species during BBN see Refs.Cyburt:2002uv ; Forestell:2018txr ; Alves:2023jlo
; Hambye:2021moy  and references therein. The co-annihilation effects become
prominent when the decay of $\chi$ is sub-dominat during first three minutes
of the cosmic history. As stated in Sec.2, the tree-level decay of $\chi$ is
avoided on the kinematic ground. However, there can be loop-induced decays as
shown in Fig.4. Within the framework of the effective theory, new interaction
terms are generated at quantum level due to these loops, which should be
regularized adding appropriate counter terms. Therefore, the decay width of
$\chi$ would depend on the residual finite terms of the new interactions
generated radiatively after renormalization at a certain order in the
perturbation theory.

Figure 4: Loop-induced decays of $\chi$

Now, we make an order of magnitude estimate of the lifetime of $\chi$
considering only the dominant one-loop decay channel, $\chi(p_{1})\rightarrow
e^{-}(p_{2})+e^{+}(p_{3})+\nu_{e}(p_{4})$ to check the validity of the
parameter space shown in the right panel of Fig.3. The contribution from the
quadratically divergent one-loop diagram shown in Fig.4 is given by,

$M_{loop}=(i~{}\frac{G_{D}G_{F}}{2})~{}\frac{m^{2}_{\chi}}{16\pi^{2}},$ (15)

where we have put $m_{\chi}$ as the cut-off scale. Now, the decay width of the
three body decay is calculated assuming vanishing masses of electron and
neutrino as well as having same energies, $E_{i}=m_{\chi}/3$.

$\displaystyle\Gamma_{\chi}\approx\frac{m_{\chi}}{512\pi^{3}}~{}|J|^{2}~{}|M_{loop}|^{2},$
(16)

where the contribution from the fermion currents is encoded in $|J|^{2}\simeq
0.6~{}m^{4}_{\chi}$. Then, the lifetime ($\tau_{\chi}=\Gamma^{-1}_{\chi}$) of
$\chi$ becomes in terms of $m_{\chi}$ and $(G_{D}/G_{F})$,

$\displaystyle\tau_{\chi}\approx
10^{5}sec\left(\frac{1~{}GeV}{m_{\chi}}\right)^{9}\left(\frac{G_{F}}{G_{D}}\right)^{2}.$
(17)

Figure 5: Constraints on $m_{\chi}$ and $G_{D}/G_{F}$ considering BBN results
as well as the dominant decay mode of $\chi$. The region below the magenta
line indicates the parameter space where the lifetime ($\tau_{\chi}$) of
$\chi$ is greater than the duration of BBN ($\tau_{BBN}\sim 3$ min). In the
left panel regions marked by the black dashed line are viable parameter space
considering the probable decay of $\chi$. The right panel is the zoomed-in
version of the left one, showing larger values for $m_{\chi}$ and
$G_{D}/G_{F}$.

Consequently, our new physics parameters are further constrained as shown in
Fig.5, modifying the earlier results as described in Fig.3 (right panel). In
particular, the pink shaded area surrounded by the magenta solid line shows
the region in which the lifetime of $\chi$ is greater than the total duration
of BBN, i.e. $\tau_{BBN}=3$ minutes. It is apparent that the for
$G_{F}/G_{D}\lesssim 10^{-2}$ the allowed parameter space remains unaffected
as shown in Fig.3 (right panel). However, for $G_{D}/G_{F}\gtrsim 10^{2}$, we
find the modified mass range of $\chi$ to be $m_{\chi}\in[0.3,0.7]~{}\rm GeV$.

$\chi$ can be cosmologically stable in the regime, $G_{D}/G_{F}\lesssim
10^{-2}$ for $m_{\chi}<<\mathcal{O}(m_{n})$. For example,
$G_{D}/G_{F}=10^{-2}$ and $m_{\chi}=100~{}\rm MeV$, $\tau_{\chi}\approx
10^{20}$ sec, which is greater than the current model-independent bound on the
lifetime of dark matter (DM), i.e. $5\times 10^{18}$ sec or $160$ Gyr
Audren:2014bca . However, for such small masses $\chi$ undergoes relativistic
freeze-out as explained in 3.1, therefore its abundance over-saturates the
density of DM in the current universe. For being a valid DM candidate, its
density has to be reduced by some other number-changing processes subsequent
to the co-annihilation with nucleons. One of the well-studied scenarios in
this regard is the cannibalism Carlson:1992fn ; Hochberg:2014dra ;
Pappadopulo:2016pkp ; Ho:2022tbw  of $\chi$ decoupled from the SM thermal bath
Ghosh:2022asg . The cannibalism can be achieved incorporating self-interaction
of $\chi$ in the effective theory.

For masses $m_{\chi}\sim\mathcal{O}(m_{n})$, and $G_{D}/G_{F}\gtrsim 10^{2}$,
$\chi$ decays after BBN, potentially leaving imprints on the cosmic micro-wave
background radiation Slatyer:2016qyl ; Acharya:2019uba  depending on the
lifetime of $\chi$.

In passing, we note that the accidental stability of $\chi$ may be ensured if
there is an exact cancellation between diagrams coming from other possible
interactions of $\chi$, e.g. similar loops with mesons. The another way to
avoid the decay is to think of an extended dark sector Carenza:2022pjd  where
$\chi$ either forms a stable bound state or becomes confined within a DM
nuclei or an atom Cline:2021itd , i.e. reminiscent of neutron surviving inside
the helium nuclei. In a nutshell, we need additional structures in our theory
to make $\chi$ absolutely stable. All these require a rather involved study of
their own, therefore are left for elsewhere.

Figure 6: DM-nucleon (left panel) and DM-electron (right panel) elastic
scattering shown by the magenta solid line in BBN-allowed parameter space
depicted in the right panel of Fig.5. The exclusion limit on DM-nucleon
scattering given by CRESST-III CRESST:2017cdd  is shown by the red dashed
line. The blue shaded region indicates the neutrino floor where the DM-nucleon
scattering becomes indistinguishable with coherent neutrino-nucleus
scattering. For DM-electron scattering, similar exclusion limit is shown using
XENON1T results XENON:2019gfn . The corresponding neutrino floor for liquid
xenon is also shown as in Ref.Carew:2023qrj .

##  4 Elastic scattering of nucleon (electron) and $\chi$

The elastic scattering between the nucleon (electron) and $\chi$ is generated
at the one-loop level as shown in Fig.7, given the co-annihilation interaction
in Eq.2. For both $\chi-n$ and $\chi-e^{-}$ scattering, the scattering
amplitudes take the following generic form.

$\displaystyle
iM=\frac{iG^{2}_{D}}{32\pi^{2}}\left(AP_{\rho}P_{\sigma}+B\frac{\Lambda^{2}}{4}\eta_{\rho\sigma}\right)\bigg{[}\bar{u}(p_{3})\,\Gamma^{\mu}\,\gamma^{\rho}\Gamma^{\nu}u(p_{1})\bigg{]}\,\bigg{[}\bar{u}(p_{4})\,\Gamma_{\mu}\,\gamma^{\sigma}\Gamma_{\nu}u(p_{2})\bigg{]},$
(18)

where $A$, $B$ are the loop factors and $\Lambda$ is the relevant cut-off
scale. $\Gamma^{\mu}$’s are of $V-A$ structure as used before and
$P^{\mu}=(p_{1}+p_{2})^{\mu}$. To note, loop factors would be different for
two scatterings. For details see Appendix A. We have taken scatterings with
neutrons only although similar scattering can be considered with protons as
well. Now, if $\chi$ becomes accidentally stable due to some hidden sector
dynamics as indicated earlier, these elastic scatterings become instrumental
in probing the co-annihilation scenario in terrestrial direct detection
experiments. In particular, for $m_{\chi}>\mathcal{O}(GeV)$, the DM-nucleon
scatterings are important Goodman:1984dc , whereas the DM-electron scatterings
put stringent constraints for sub-GeV DM masses Essig:2011nj . In Fig.6, we
have shown these scattering cross-sections, allowed from BBN prediction and
compared with the present experimental bounds.

At first, the spin-independent part of DM-nucleon scattering is extracted from
the generic expression shown in Eq.18, as the spin-dependent part is always
suppressed in presence of the former one. Therefore, the relevant cross-
section (detailed in Appendix A) is given by

$\displaystyle\sigma_{\chi
n}\approx\frac{G^{4}_{D}}{\pi^{5}}\bigg{(}m_{\chi}m_{n}(m_{\chi}+m_{n})\bigg{)}^{2}$
(19)

In our scenario, $\chi$ might account for low mass DM ($m_{\chi}\lesssim
3~{}\rm GeV$) for which the bound on the spin-independent DM-nucleon cross-
section is rather weak for XENON1T-like experiments XENON:2019gfn . This type
of scenario can be probed at CRESST-III experiments CRESST:2019jnq  (and via
Migdal effect Bell:2019egg ; Dolan:2017xbu ) of which the current bound on the
DM-nucleon scattering cross-section is $\sigma_{\chi n}\lesssim 10^{-38}~{}\rm
cm^{2}$ for $m_{\chi}\lesssim 1\,\rm GeV$.

In the left panel of Fig.6, we have shown the detection prospect of BBN-
allowed parameter space for $\chi$ particles. There are two regions of
interest as in previous section, namely $m_{\chi}<<\mathcal{O}(m_{n})$ and
$m_{\chi}\sim\mathcal{O}(m_{n})$. For $m_{\chi}<<\mathcal{O}(m_{n})$, the BBN
prediction allows only small $G_{D}/G_{F}(\lesssim 10^{-2})$ values, thereby
$\sigma_{\chi n}$ is well below the neutrino floor (shaded region below the
blue solid line). For example, with $m_{\chi}=0.05~{}\rm GeV$ and
$G_{D}/G_{F}=10^{-2}$, the cross-section becomes $\sigma_{\chi n}=5\times
10^{-61}~{}\rm cm^{2}$. For $m_{\chi}\sim\mathcal{O}(m_{n})$, the BBN
prediction allows for both small and large values of $G_{D}/G_{F}$. Therefore,
the detection possibility of $\chi$ becomes more optimistic for
$m_{\chi}\in[0.3,0.7]~{}\rm GeV$ and $G_{D}/G_{F}\gtrsim 10^{2}$. For a
benchmark point, $m_{\chi}=0.6~{}\rm GeV$ and $G_{D}/G_{F}=250$, we get
$\sigma_{\chi n}\approx 9.8\times 10^{-41}~{}\rm cm^{2}$, which is evidently
above the neutrino floor. The magenta solid line in Fig.6 depicts the coveted
sensitivity of the future experiments to probe the co-annihilation scenario,
being consistent with the cosmological observations. The current exclusion
limit on the DM-nucleon cross-section coming from the preliminary results of
CRESST-III phase-I experiments CRESST:2017cdd  have been indicated by the red
dashed line and the shaded region above the line.

Now, the DM-electron scattering cross-section can be calculated (detailed in
Appendix A ) in terms of DM-nucleon scatterings as both type of scatterings
are generated radiatively from the same effective co-annihilation operator.
The DM-electron scattering cross-section is suppressed compared to the DM-
nucleon scattering due to small electron mass, i.e.

$\displaystyle\sigma_{\chi
e}\approx\frac{m^{2}_{e}}{m^{2}_{\chi}m^{2}_{n}}(m_{\chi}+m_{n})^{2}\sigma_{\chi
n}.$ (20)

In the right panel of Fig.6, we have shown the BBN-allowed parameter space for
DM-electron scattering, contrasted with results from ionization signal of
XENON1T experiments XENON:2019gfn . Unlike the DM-nucleon scattering, the
allowed contour in $\sigma_{\chi e}-m_{\chi}$ plane is within neutrino floor,
making the distinction between the signal and background events difficult. The
dominant background comes from coherent neutrino-nucleus scattering, shown as
the neutrino floor. We have adopted the neutrino floor line as in Fig.1 (left
panel) of Ref.Carew:2023qrj . However, the neutrino floor (sometimes called as
neutrino ‘fog’OHare:2021utq ) may not be hard, as with long exposure
experiments the energy profile of signal events might be discernible from its
background Wyenberg:2018eyv ; Carew:2023qrj .

In a nutshell, if $\chi$ survives accidentally till date, it can be probed in
future direct detection experiments with improved sensitivities and long
exposure time.

## Acknowledgment

We thank Satyanarayan Mukhopadhyay for various helpful conversations since the
outset of this work. We also thank Rohan Pramanick, Sourav Gope and Utpal
Chattopadhyay for computational help and Sougata Ganguly, Avirup Ghosh and
Tarak Nath Maity for interesting discussions. This work is supported by the
Institute Fellowship provided by the Indian Association for the Cultivation of
Science (IACS), Kolkata.

##  Appendix A Calculation of loop-induced scattering of neutron (electron)
and $\chi$

To start with, we had only operator that enables the co-annihilation of $\chi$
and neutron producing a proton and an electron in the final state. This
operator itself induces an effective operator at one-loop level, responsible
for the elastic scattering of neutron (electron) and $\chi$ as shown in Fig.7,
i.e. relevant for direct detection experiments.

Figure 7: Feynman graphs for elastic scattering of neutron (electron) and
$\chi$ generated at one-loop level from the co-annihilation operator described
in Eq.2.

The amplitude for elastic scattering, $n(p_{1})+\chi(p_{2})\rightarrow
n(p_{3})+\chi(p_{4})$ calculated from Eq.2 is given by,

$\displaystyle
iM=\left(-i\frac{G_{D}}{\sqrt{2}}\right)^{2}\int\frac{d^{4}k}{(2\pi)^{4}}\bigg{[}\bar{u}(3)\,\Gamma^{\mu}\,\frac{i(\not{k}+m_{p})}{k^{2}-m^{2}_{p}}\,\Gamma^{\nu}u(1)\bigg{]}\,\bigg{[}\bar{u}(4)\,\Gamma_{\mu}\,\frac{i(\not{P}-\not{k}+m_{e})}{(P-k)^{2}-m^{2}_{e}}\,\Gamma_{\nu}u(2)\bigg{]}$
(21)

where $u(i)=u(p_{i})$ and $P^{2}=(p_{1}+p_{2})^{2}=s$, while $p_{i}$ is the
on-shell four momentum. For $\Gamma^{a}=\gamma^{a}(I-\gamma_{5})$, the form
used for all calculations, we get $\Gamma^{a}\Gamma^{b}=0$. Then the amplitude
takes rather simpler form, i.e.

$\displaystyle
iM=\left(\frac{G_{D}}{\sqrt{2}}\right)^{2}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{k_{\rho}\left(P-k\right)_{\sigma}}{(k^{2}-m^{2}_{p})((P-k)^{2}-m^{2}_{e})}\bigg{[}\bar{u}(3)\,\Gamma^{\mu}\,\gamma^{\rho}\Gamma^{\nu}u(2)\bigg{]}\,\bigg{[}\bar{u}(4)\,\Gamma_{\mu}\,\gamma^{\sigma}\Gamma_{\nu}u(1)\bigg{]}$
(22)

The above integral is divergent in general, which is expected in an effective
theory with coupling of negative mass dimension. Therefore, we put a cut-off
scale, $\Lambda$ and calculate the amplitude as

$\displaystyle
iM=\frac{iG^{2}_{D}}{32\pi^{2}}\left(AP_{\rho}P_{\sigma}+B\frac{\Lambda^{2}}{4}\eta_{\rho\sigma}\right)\bigg{[}\bar{u}(3)\,\Gamma^{\mu}\,\gamma^{\rho}\Gamma^{\nu}u(1)\bigg{]}\,\bigg{[}\bar{u}(4)\,\Gamma_{\mu}\,\gamma^{\sigma}\Gamma_{\nu}u(2)\bigg{]}$
(23)

where $A$ and $B$ given by Eq.24 are complex numbers in general, to be
calculated numerically.

$\displaystyle A$
$\displaystyle=\int^{1}_{0}dx\,x(1-x)F_{1}(x),~{}~{}B=\int^{1}_{0}dx\,F_{2}(x)$
$\displaystyle F_{1}(x)$
$\displaystyle=\log\left(\frac{\Lambda^{2}+\Delta}{\Delta}\right)-\frac{\Lambda^{2}}{\Lambda^{2}+\Delta}$
$\displaystyle F_{2}(x)$
$\displaystyle=1+\frac{\Delta}{\Lambda^{2}+\Delta}-\frac{2\Delta}{\Lambda^{2}}\log\left(\frac{\Lambda^{2}+\Delta}{\Delta}\right)$
$\displaystyle\Delta$ $\displaystyle=m^{2}_{e}(1-x)+m^{2}_{p}\,x-x(1-x)P^{2}$
(24)

Now, this is the most general structure for elastic scattering at hand from
our effective theory. In direct detection experiments, the stringent bounds
comes from the spin-independent operators. Therefore, we extract the spin-
independent part ($M_{SI}$) from the above expression as the following using
the relation,
$\Gamma^{\mu}\gamma^{\rho}\Gamma^{\nu}=4\gamma^{\mu}\eta^{\rho\nu}(I-\gamma_{5})-2\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}(I-\gamma_{5})$.

$\displaystyle iM$ $\displaystyle=$
$\displaystyle\frac{iG^{2}_{D}}{32\pi^{2}}\left(AP_{\rho}P_{\sigma}+B\frac{\Lambda^{2}}{4}\eta_{\rho\sigma}\right)\bigg{[}\bar{u}(3)\bigg{(}4\gamma^{\mu}\eta^{\rho\nu}(I-\gamma_{5})-2\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}(I-\gamma_{5})\bigg{)}u(1)\bigg{]}$
$\displaystyle\hskip
113.81102pt\bigg{[}\bar{u}(4)\bigg{(}4\gamma_{\mu}\delta^{\sigma}_{\nu}(I-\gamma_{5})-2\gamma_{\mu}\gamma_{\nu}\gamma^{\sigma}(I-\gamma_{5})\bigg{)}u(2)\bigg{]}$
$\displaystyle iM_{SI}$ $\displaystyle=$
$\displaystyle\frac{iG^{2}_{D}}{2\pi^{2}}\left(AP^{2}+B\Lambda^{2}\right)\bigg{[}\bar{u}(3)\gamma^{\mu}u(1)\bigg{]}\bigg{[}\bar{u}(4)\gamma_{\mu}u(2)\bigg{]}$
(25)

In direct detection experiments, this elastic scattering is considered in the
non-relativistic regime, in which the momentum transfer is small compared to
particle masses. Therefore, we can approximate Dirac fermion as Lin:2019uvt ,
$u(k)\approx\sqrt{2m}(\xi_{s}~{}\xi_{s})^{T}$, where $\xi_{s}$ is a spinor
with $\sum_{s}\xi_{s}\xi^{\dagger}_{s}=I$. Consequently, the vector operator
in 25 can be approximated as the scalar one, i.e.

$\displaystyle\bar{u}(p)\gamma^{\mu}u(p^{\prime})\rightarrow
2m~{}\xi^{\dagger}_{s^{\prime}}\xi_{s}$ (26)

Now, setting the cut-off scale to be $\Lambda=(m_{n}+m_{\chi})$, which is also
the center-of-mass energy ($\sqrt{P^{2}}$), we can write the spin-independent
effective operator as the following Lin:2019uvt ; Lisanti:2016jxe .

$\displaystyle
M_{SI}=\frac{G^{2}_{D}}{\pi^{2}}(m_{n}+m_{\chi})^{2}\left(4m_{n}m_{\chi}~{}\xi^{\dagger}_{s^{\prime}}(n)\xi_{s}(n)\xi^{\dagger}_{t^{\prime}}(\chi)\xi_{t}(\chi)\right),$
(27)

where $A,B$ is taken as $\mathcal{O}(1)$ numbers. Now, the scattering cross-
section in the non-relativistic assumption is given by,

$\displaystyle\sigma_{\chi n}$
$\displaystyle=\frac{1}{4m_{\chi}m_{n}v}\int\frac{d^{3}p_{3}}{(2\pi)^{3}2m_{n}}\frac{d^{3}p_{2}}{(2\pi)^{3}2m_{\chi}}(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4})\overline{|M_{SI}|^{2}}$
$\displaystyle=\frac{\overline{|M_{SI}|^{2}}}{32\pi
m^{2}_{\chi}m^{2}_{n}v^{2}}\int^{2\mu_{\chi
n}v}_{0}p_{3}~{}dp_{3}\int~{}d\cos\theta~{}\delta\left(\cos\theta-\frac{p_{3}}{2\mu_{\chi
n}v}\right)$ (28)

where the spin-averaged probability of the elastic scattering,
$\overline{|M_{SI}|^{2}}=\frac{G^{4}_{D}}{\pi^{4}}(4m_{\chi}m_{n}(m_{\chi}+m_{n})^{2})^{2}$
and $v$ is the initial velocity of $\chi$ in the rest frame of nucleon. The
momentum of the recoiled nucleon in the non-relativistic approximation is
given by,$|p_{3}|=2\mu_{\chi n}v\cos\theta$, hence comes the upper limit of
$p_{3}$ shown in Eq.28. $\mu_{\chi n}(=m_{\chi}m_{n}/(m_{\chi}+m_{n}))$ is the
reduced mass of the system . Finally, the spin-independent $\chi-n$ elastic
cross-section ($\sigma_{\chi n}$) turns out to be

$\displaystyle\sigma_{\chi
n}\approx\frac{G^{4}_{D}}{\pi^{5}}\bigg{(}m_{\chi}m_{n}(m_{\chi}+m_{n})\bigg{)}^{2}.$
(29)

Now, we can also calculate $\chi-e$ scattering at the one-loop level using
same methodology, which produces similar result as in 22, replacing $m_{e}$ by
$m_{n}$ in one of the propagators. The loop factors are also modified as
$\Delta$ becomes $\tilde{\Delta}=m^{2}_{n}(1-x)+m^{2}_{p}x-x(1-x)P^{2}\approx
m^{2}_{n}-x(1-x)P^{2}$. Following similar steps as earlier ( 25-28), we get
the spin-independent $\chi-e$ elastic cross-section ($\sigma_{\chi e}$) to be,

$\displaystyle\sigma_{\chi
e}\approx\frac{G^{4}_{D}}{\pi^{5}}\frac{m^{2}_{\chi}m^{2}_{e}}{(m_{\chi}+m_{e})^{2}}(m_{\chi}+m_{n})^{4}.$
(30)

To note, we have set same cut-off scale as before, which brings the factor of
$(m_{\chi}+m_{n})^{4}$. The ratio of two scattering cross-sections for
$m_{\chi}>>m_{e}$ becomes

$\displaystyle\frac{\sigma_{\chi e}}{\sigma_{\chi
n}}\approx\frac{m^{2}_{e}}{m^{2}_{\chi}m^{2}_{n}}(m_{\chi}+m_{n})^{2}.$ (31)

## References

  * (1) E. Aver, K. A. Olive and E. D. Skillman, “The effects of He I $\lambda$10830 on helium abundance determinations,” JCAP 07, 011 (2015) doi:10.1088/1475-7516/2015/07/011 [arXiv:1503.08146 [astro-ph.CO]].
  * (2) E. W. Kolb and M. S. Turner, “The Early Universe,” Front. Phys. 69, 1-547 (1990) doi:10.1201/9780429492860
  * (3) S. Weinberg, “Cosmology,” Oxford, UK: Oxford Univ. Pr. (2008) 593 p
  * (4) E. W. Kolb, M. S. Turner and T. P. Walker, “The Effect of Interacting Particles on Primordial Nucleosynthesis,” Phys. Rev. D 34, 2197 (1986) doi:10.1103/PhysRevD.34.2197
  * (5) C. Boehm, M. J. Dolan and C. McCabe, “A Lower Bound on the Mass of Cold Thermal Dark Matter from Planck,” JCAP 08, 041 (2013) doi:10.1088/1475-7516/2013/08/041 [arXiv:1303.6270 [hep-ph]].
  * (6) M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu, “Revisiting Big-Bang Nucleosynthesis Constraints on Long-Lived Decaying Particles,” Phys. Rev. D 97, no.2, 023502 (2018) doi:10.1103/PhysRevD.97.023502 [arXiv:1709.01211 [hep-ph]].
  * (7) E. Holtmann, M. Kawasaki, K. Kohri and T. Moroi, “Radiative decay of a longlived particle and big bang nucleosynthesis,” Phys. Rev. D 60, 023506 (1999) doi:10.1103/PhysRevD.60.023506 [arXiv:hep-ph/9805405 [hep-ph]].
  * (8) M. H. Reno and D. Seckel, “Primordial Nucleosynthesis: The Effects of Injecting Hadrons,” Phys. Rev. D 37, 3441 (1988) doi:10.1103/PhysRevD.37.3441
  * (9) S. Dimopoulos, R. Esmailzadeh, L. J. Hall and G. D. Starkman, “Limits on Late Decaying Particles From Nucleosynthesis,” Nucl. Phys. B 311, 699-718 (1989) doi:10.1016/0550-3213(89)90173-9
  * (10) M. Kusakabe, T. Kajino, T. Yoshida, T. Shima, Y. Nagai and T. Kii, “New Constraints on Radiative Decay of Long-Lived X Particles in Big Bang Nucleosynthesis with New Rates of Photodisintegration Reactions of 4He,” Phys. Rev. D 79, 123513 (2009) doi:10.1103/PhysRevD.79.123513 [arXiv:0806.4040 [astro-ph]].
  * (11) T. H. Yeh, K. A. Olive and B. D. Fields, “Limits on Non-Relativistic Matter During Big-Bang Nucleosynthesis,” [arXiv:2401.08795 [astro-ph.CO]].
  * (12) K. Jedamzik, “Big bang nucleosynthesis constraints on hadronically and electromagnetically decaying relic neutral particles,” Phys. Rev. D 74, 103509 (2006) doi:10.1103/PhysRevD.74.103509 [arXiv:hep-ph/0604251 [hep-ph]].
  * (13) E. W. Kolb and R. J. Scherrer, “Massive Neutrinos and Primordial Nucleosynthesis,” Phys. Rev. D 25, 1481 (1982) doi:10.1103/PhysRevD.25.1481
  * (14) M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu, “Revisiting Big-Bang Nucleosynthesis Constraints on Dark-Matter Annihilation,” Phys. Lett. B 751, 246-250 (2015) doi:10.1016/j.physletb.2015.10.048 [arXiv:1509.03665 [hep-ph]].
  * (15) P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg and S. Wild, “BBN constraints on the annihilation of MeV-scale dark matter,” JCAP 04, 029 (2019) doi:10.1088/1475-7516/2019/04/029 [arXiv:1901.06944 [hep-ph]].
  * (16) J. Hisano, M. Kawasaki, K. Kohri, T. Moroi and K. Nakayama, “Cosmic Rays from Dark Matter Annihilation and Big-Bang Nucleosynthesis,” Phys. Rev. D 79, 083522 (2009) doi:10.1103/PhysRevD.79.083522 [arXiv:0901.3582 [hep-ph]].
  * (17) M. Kusakabe, M. K. Cheoun and K. S. Kim, “General limit on the relation between abundances of D and 7Li in big bang nucleosynthesis with nucleon injections,” Phys. Rev. D 90, no.4, 045009 (2014) doi:10.1103/PhysRevD.90.045009 [arXiv:1404.3090 [astro-ph.CO]].
  * (18) D. Albornoz Vasquez, A. Belikov, A. Coc, J. Silk and E. Vangioni, “Neutron injection during primordial nucleosynthesis alleviates the primordial 7Li problem,” Phys. Rev. D 86, 063501 (2012) doi:10.1103/PhysRevD.86.063501 [arXiv:1208.0443 [astro-ph.CO]].
  * (19) K. M. Nollett and G. Steigman, Phys. Rev. D 89, no.8, 083508 (2014) doi:10.1103/PhysRevD.89.083508 [arXiv:1312.5725 [astro-ph.CO]].
  * (20) R. J. Protheroe, T. Stanev and V. S. Berezinsky, Phys. Rev. D 51, 4134-4144 (1995) doi:10.1103/PhysRevD.51.4134 [arXiv:astro-ph/9409004 [astro-ph]].
  * (21) M. Kawasaki and T. Moroi, Astrophys. J. 452, 506 (1995) doi:10.1086/176324 [arXiv:astro-ph/9412055 [astro-ph]].
  * (22) K. Griest and D. Seckel, “Three exceptions in the calculation of relic abundances,” Phys. Rev. D 43, 3191-3203 (1991) doi:10.1103/PhysRevD.43.3191
  * (23) T. Jittoh, K. Kohri, M. Koike, J. Sato, T. Shimomura and M. Yamanaka, “Big-bang nucleosynthesis and the relic abundance of dark matter in a stau-neutralino coannihilation scenario,” Phys. Rev. D 78, 055007 (2008) doi:10.1103/PhysRevD.78.055007 [arXiv:0805.3389 [hep-ph]].
  * (24) T. Jittoh, K. Kohri, M. Koike, J. Sato, T. Shimomura and M. Yamanaka, “Stau relic density at the Big-Bang nucleosynthesis era consistent with the abundance of the light element nuclei in the coannihilation scenario,” Phys. Rev. D 82, 115030 (2010) doi:10.1103/PhysRevD.82.115030 [arXiv:1001.1217 [hep-ph]].
  * (25) R. T. D’Agnolo, C. Mondino, J. T. Ruderman and P. J. Wang, “Exponentially Light Dark Matter from Coannihilation,” JHEP 08, 079 (2018) doi:10.1007/JHEP08(2018)079 [arXiv:1803.02901 [hep-ph]].
  * (26) J. Ellis, B. D. Fields, F. Luo, K. A. Olive and V. C. Spanos, “An Enhanced Cosmological Li6 Abundance as a Potential Signature of Residual Dark Matter Annihilations,” Phys. Rev. D 84, 123502 (2011) doi:10.1103/PhysRevD.84.123502 [arXiv:1109.0549 [astro-ph.CO]].
  * (27) B. Belfatto, D. Buttazzo, C. Gross, P. Panci, A. Strumia, N. Vignaroli, L. Vittorio and R. Watanabe, “Dark Matter abundance via thermal decays and leptoquark mediators,” [arXiv:2111.14808 [hep-ph]].
  * (28) A. Strumia, “Dark Matter interpretation of the neutron decay anomaly,” JHEP 02, 067 (2022) doi:10.1007/JHEP02(2022)067 [arXiv:2112.09111 [hep-ph]].
  * (29) J. A. Dror, G. Elor and R. Mcgehee, “Absorption of Fermionic Dark Matter by Nuclear Targets,” JHEP 02, 134 (2020) doi:10.1007/JHEP02(2020)134 [arXiv:1908.10861 [hep-ph]].
  * (30) A. J. Long, C. Lunardini and E. Sabancilar, “Detecting non-relativistic cosmic neutrinos by capture on tritium: phenomenology and physics potential,” JCAP 08, 038 (2014) doi:10.1088/1475-7516/2014/08/038 [arXiv:1405.7654 [hep-ph]].
  * (31) A. G. Cocco, G. Mangano and M. Messina, “Probing low energy neutrino backgrounds with neutrino capture on beta decaying nuclei,” JCAP 06, 015 (2007) doi:10.1088/1475-7516/2007/06/015 [arXiv:hep-ph/0703075 [hep-ph]].
  * (32) J. M. Cline, “Dark atoms and composite dark matter,” doi:10.21468/SciPostPhysLectNotes.52 [arXiv:2108.10314 [hep-ph]].
  * (33) D. McKeen, M. Pospelov and N. Raj, “Cosmological and astrophysical probes of dark baryons,” Phys. Rev. D 103, no.11, 115002 (2021) doi:10.1103/PhysRevD.103.115002 [arXiv:2012.09865 [hep-ph]].
  * (34) P. Carenza, R. Pasechnik, G. Salinas and Z. W. Wang, “Glueball Dark Matter Revisited,” Phys. Rev. Lett. 129, no.26, 26 (2022) doi:10.1103/PhysRevLett.129.261302 [arXiv:2207.13716 [hep-ph]].
  * (35) B. D. Fields, K. A. Olive, T. H. Yeh and C. Young, “Big-Bang Nucleosynthesis after Planck,” JCAP 03, 010 (2020) [erratum: JCAP 11, E02 (2020)] doi:10.1088/1475-7516/2020/03/010 [arXiv:1912.01132 [astro-ph.CO]].
  * (36) R. An, V. Gluscevic, E. Calabrese and J. C. Hill, “What does cosmology tell us about the mass of thermal-relic dark matter?,” [arXiv:2202.03515 [astro-ph.CO]].
  * (37) M. Fukugita and T. Yanagida, “Baryogenesis Without Grand Unification,” Phys. Lett. B 174, 45-47 (1986) doi:10.1016/0370-2693(86)91126-3
  * (38) V. A. Rubakov and M. E. Shaposhnikov, “Electroweak baryon number nonconservation in the early universe and in high-energy collisions,” Usp. Fiz. Nauk 166, 493-537 (1996) doi:10.1070/PU1996v039n05ABEH000145 [arXiv:hep-ph/9603208 [hep-ph]].
  * (39) R. V. Wagoner, W. A. Fowler and F. Hoyle, “On the Synthesis of elements at very high temperatures,” Astrophys. J. 148, 3-49 (1967) doi:10.1086/149126
  * (40) V. A. Rubakov and D. S. Gorbunov, “Introduction to the Theory of the Early Universe: Hot big bang theory,” doi:10.1142/10447
  * (41) J. Bernstein, L. S. Brown and G. Feinberg, “COSMOLOGICAL HELIUM PRODUCTION SIMPLIFIED,” Rev. Mod. Phys. 61, 25 (1989) doi:10.1103/RevModPhys.61.25
  * (42) V. F. Mukhanov, “Nucleosynthesis without a computer,” Int. J. Theor. Phys. 43, 669-693 (2004) doi:10.1023/B:IJTP.0000048169.69609.77 [arXiv:astro-ph/0303073 [astro-ph]].
  * (43) R. H. Cyburt, J. R. Ellis, B. D. Fields and K. A. Olive, “Updated nucleosynthesis constraints on unstable relic particles,” Phys. Rev. D 67, 103521 (2003) doi:10.1103/PhysRevD.67.103521 [arXiv:astro-ph/0211258 [astro-ph]].
  * (44) L. Forestell, D. E. Morrissey and G. White, “Limits from BBN on Light Electromagnetic Decays,” JHEP 01, 074 (2019) doi:10.1007/JHEP01(2019)074 [arXiv:1809.01179 [hep-ph]].
  * (45) J. R. Alves, L. Angel, L. Guedes, R. M. P. Neves, F. S. Queiroz, D. R. da Silva, R. Silva and Y. Villamizar, “Updated Big Bang Nucleosynthesis Bounds on Long-lived Particles from Dark Sectors,” [arXiv:2311.07688 [hep-ph]].
  * (46) T. Hambye, M. Hufnagel and M. Lucca, “Cosmological constraints on the decay of heavy relics into neutrinos,” JCAP 05, no.05, 033 (2022) doi:10.1088/1475-7516/2022/05/033 [arXiv:2112.09137 [hep-ph]].
  * (47) B. Audren, J. Lesgourgues, G. Mangano, P. D. Serpico and T. Tram, “Strongest model-independent bound on the lifetime of Dark Matter,” JCAP 12, 028 (2014) doi:10.1088/1475-7516/2014/12/028 [arXiv:1407.2418 [astro-ph.CO]].
  * (48) E. D. Carlson, M. E. Machacek and L. J. Hall, “Self-interacting dark matter,” Astrophys. J. 398, 43-52 (1992) doi:10.1086/171833
  * (49) A. Ghosh, S. Gope and S. Mukhopadhyay, “Cannibal dark matter decoupled from the standard model: Cosmological constraints,” Phys. Rev. D 106, no.10, 103515 (2022) doi:10.1103/PhysRevD.106.103515 [arXiv:2206.11046 [hep-ph]].
  * (50) Y. Hochberg, E. Kuflik, T. Volansky and J. G. Wacker, “Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 113, 171301 (2014) doi:10.1103/PhysRevLett.113.171301 [arXiv:1402.5143 [hep-ph]].
  * (51) D. Pappadopulo, J. T. Ruderman and G. Trevisan, Phys. Rev. D 94, no.3, 035005 (2016) doi:10.1103/PhysRevD.94.035005 [arXiv:1602.04219 [hep-ph]].
  * (52) S. Y. Ho, “An asymmetric SIMP dark matter model,” JHEP 10, 182 (2022) doi:10.1007/JHEP10(2022)182 [arXiv:2207.13373 [hep-ph]].
  * (53) T. R. Slatyer and C. L. Wu, “General Constraints on Dark Matter Decay from the Cosmic Microwave Background,” Phys. Rev. D 95, no.2, 023010 (2017) doi:10.1103/PhysRevD.95.023010 [arXiv:1610.06933 [astro-ph.CO]].
  * (54) S. K. Acharya and R. Khatri, “CMB anisotropy and BBN constraints on pre-recombination decay of dark matter to visible particles,” JCAP 12, 046 (2019) doi:10.1088/1475-7516/2019/12/046 [arXiv:1910.06272 [astro-ph.CO]].
  * (55) M. W. Goodman and E. Witten, “Detectability of Certain Dark Matter Candidates,” Phys. Rev. D 31, 3059 (1985) doi:10.1103/PhysRevD.31.3059
  * (56) R. Essig, J. Mardon and T. Volansky, “Direct Detection of Sub-GeV Dark Matter,” Phys. Rev. D 85, 076007 (2012) doi:10.1103/PhysRevD.85.076007 [arXiv:1108.5383 [hep-ph]].
  * (57) A. H. Abdelhameed et al. [CRESST], “First results from the CRESST-III low-mass dark matter program,” Phys. Rev. D 100, no.10, 102002 (2019) doi:10.1103/PhysRevD.100.102002 [arXiv:1904.00498 [astro-ph.CO]].
  * (58) F. Petricca et al. [CRESST], “First results on low-mass dark matter from the CRESST-III experiment,” J. Phys. Conf. Ser. 1342, no.1, 012076 (2020) doi:10.1088/1742-6596/1342/1/012076 [arXiv:1711.07692 [astro-ph.CO]].
  * (59) N. F. Bell, J. B. Dent, J. L. Newstead, S. Sabharwal and T. J. Weiler, “Migdal effect and photon bremsstrahlung in effective field theories of dark matter direct detection and coherent elastic neutrino-nucleus scattering,” Phys. Rev. D 101, no.1, 015012 (2020) doi:10.1103/PhysRevD.101.015012 [arXiv:1905.00046 [hep-ph]].
  * (60) M. J. Dolan, F. Kahlhoefer and C. McCabe, “Directly detecting sub-GeV dark matter with electrons from nuclear scattering,” Phys. Rev. Lett. 121, no.10, 101801 (2018) doi:10.1103/PhysRevLett.121.101801 [arXiv:1711.09906 [hep-ph]].
  * (61) E. Aprile et al. [XENON], “Light Dark Matter Search with Ionization Signals in XENON1T,” Phys. Rev. Lett. 123, no.25, 251801 (2019) doi:10.1103/PhysRevLett.123.251801 [arXiv:1907.11485 [hep-ex]].
  * (62) B. Carew, A. R. Caddell, T. N. Maity and C. A. J. O’Hare, “Neutrino fog for dark matter-electron scattering experiments,” Phys. Rev. D 109, no.8, 083016 (2024) doi:10.1103/PhysRevD.109.083016 [arXiv:2312.04303 [hep-ph]].
  * (63) C. A. J. O’Hare, “New Definition of the Neutrino Floor for Direct Dark Matter Searches,” Phys. Rev. Lett. 127, no.25, 251802 (2021) doi:10.1103/PhysRevLett.127.251802 [arXiv:2109.03116 [hep-ph]].
  * (64) J. Wyenberg and I. M. Shoemaker, “Mapping the neutrino floor for direct detection experiments based on dark matter-electron scattering,” Phys. Rev. D 97, no.11, 115026 (2018) doi:10.1103/PhysRevD.97.115026 [arXiv:1803.08146 [hep-ph]].
  * (65) T. Lin, “Dark matter models and direct detection,” PoS 333, 009 (2019) doi:10.22323/1.333.0009 [arXiv:1904.07915 [hep-ph]].
  * (66) M. Lisanti, “Lectures on Dark Matter Physics,” doi:10.1142/9789813149441_0007 [arXiv:1603.03797 [hep-ph]].


# Characteristic approach
to the soliton resolution

Piotr Bizoń  Institute of Theoretical Physics, Jagiellonian University, Kraków
<EMAIL_ADDRESS>,  Bradley Cownden  Institute of Theoretical Physics,
Jagiellonian University, Kraków<EMAIL_ADDRESS>and  Maciej
Maliborski  University of Vienna, Faculty of Mathematics, Oskar-Morgenstern-
Platz 1, 1090 Vienna, Austria, and University of Vienna, Gravitational
Physics, Boltzmanngasse 5, 1090 Vienna, Austria
<EMAIL_ADDRESS>
###### Abstract.

As a toy model for understanding the soliton resolution phenomenon we consider
a characteristic initial boundary value problem for the 4$d$ equivariant Yang-
Mills equation outside a ball. Our main objective is to illustrate the
advantages of employing outgoing null (or asymptotically null) foliations in
analyzing the relaxation processes due to the dispersal of energy by
radiation. In particular, within this approach it is evident that the endstate
of evolution must be non-radiative (meaning vanishing flux of energy at future
null infinity). In our toy model such non-radiative configurations are given
by a static solution (called the half-kink) plus an alternating chain of $N$
decoupled kinks and antikinks. We show numerically that the configurations
$N=0$ (static half-kink) and $N=1$ (superposition of the static half-kink and
the antikink which recedes to infinity) appear as generic attractors and we
determine a codimension-one borderline between their basins of attraction. The
rates of convergence to these attractors are analyzed in detail.

##  1\. Introduction

According to the soliton resolution conjecture, global-in-time generic
solutions of nonlinear dispersive wave equations resolve for
$t\rightarrow\infty$ into a superposition of decoupled nonlinear bound states
(solitons) and radiation [1]. There are numerous physical manifestations of
this phenomenon, ranging from the formation of solitons in optical fibers
(modelled by a nonlinear Schrödinger equation) to the formation of stationary
black holes in binary black hole mergers (modelled by the Einstein equation).

The past decade has seen a significant progress in mathematical understanding
of the soliton resolution, especially for radial solutions of the energy
critical nonlinear wave equations (see [2] for a survey and references
therein). Notably, the soliton resolution was recently proved for the
equivariant wave maps $\mathbb{R}^{2+1}\rightarrow\mathbb{S}^{2}$ [3, 4] and
the equivariant Yang-Mills (YM) equation in $4+1$ dimensions [4] (both for the
global and blowup solutions). These remarkable results are abstract in the
sense that they enumerate all possible asymptotic scenarios but do not settle
which scenarios are actually realized.

As a toy model for more quantitative description of the soliton resolution
phenomenon, in this paper we consider the equivariant YM equation in $4+1$
dimensions

$W_{tt}=W_{rr}+\frac{1}{r}W_{r}+\frac{2}{r^{2}}\,W(1-W^{2}),$ (1)

where $W(t,r)$ is the YM potential. As the spatial domain we take the exterior
of the unit ball, i.e. $r\geq 1$, and impose the Dirichlet condition on the
boundary $W(t,r=1)=0$. The associated conserved energy is

$E[W]=\frac{1}{2}\int_{1}^{\infty}\left(W_{t}^{2}+W_{r}^{2}+\frac{(1-W^{2})^{2}}{r^{2}}\right)rdr.$
(2)

Finiteness of energy requires that $|W(t,\infty)|=1$. Since the singular point
$r=0$ is outside the domain, it is easy to see that solutions starting from
smooth, finite-energy initial data ($W(0,r),W_{t}(0,r)$), which are compatible
with the boundary condition, remain smooth for all times. Our goal is to
describe their asymptotic behavior for $t\rightarrow\infty$.

The key role in our analysis will be played by the half-kink

$Q(r)=\frac{r^{2}-1}{r^{2}+1}\,,$ (3)

which is the unique (modulo sign) static solution of equation (1) satisfying
the boundary condition $Q(1)=0$. The half-kink is a global minimizer of energy
(see the Bogomolnyi inequality (20) below) and thereby a natural candidate for
an attractor. Indeed, we will see that on any compact interval $[1,R)$ every
smooth solution $W(t,r)$ converges to $Q(r)$ or $-Q(r)$ as
$t\rightarrow\infty$. However, for sufficiently large energies a nontrivial
coherent structure can simultaneously develop in the asymptotic region
$(R,\infty)$. This behavior is intimately related to the energy criticality of
the model and is absent in supercritical dimensions; cf. the soliton
resolution for the equivariant wave maps exterior to a ball in $3+1$
dimensions [5, 6]. The paper is organized as follows. In section 2 we first
recall from [4] the formulation of the soliton resolution conjecture for
equation (1) in the whole space. Then we present an analogous conjecture in
our model and support one special case by a heuristic argument based on the
method of collective coordinates. In section 3 we reformulate the initial-
boundary problem in terms of null foliations of constant retarded time and
compactified spatial domain. Using this formulation, in section 4 we consider
the late-time behavior (the quasinormal ringdown and the polynomial tail) for
the linearized problem. Finally, in section 5 we present numerical evidence
supporting the soliton resolution conjecture.

##  2\. Soliton resolution

Let us first recall what is known about equation (1) posed on the whole space
$r\geq 0$. In this case, equation (1) is invariant under scaling, i.e. if
$W(t,r)$ is a solution, so is $W_{\lambda}(t,r)=W(t/\lambda,r/\lambda)$ for
any positive number $\lambda$. The conserved energy

$E_{0}[W]=\frac{1}{2}\int_{0}^{\infty}\left(W_{t}^{2}+W_{r}^{2}+\frac{(1-W^{2})^{2}}{r^{2}}\right)rdr$
(4)

is scale invariant, i.e. $E_{0}[W_{\lambda}]=E_{0}[W]$, which is a distinctive
feature of the critical dimension $d=4$ allowing for the existence of
nontrivial static solutions in the presence of scaling symmetry. These static
solutions, hereafter called kinks (also referred to in the literature as
instantons, solitons, or bubbles), form a one-parameter family

$Q_{\lambda}(r)=\frac{r^{2}-\lambda^{2}}{r^{2}+\lambda^{2}}$ (5)

with energy $E_{0}[Q_{\lambda}]=\frac{4}{3}$ which is the minimum energy for
solutions interpolating between different vacuum states $W=\pm 1$ at the
origin and at infinity. Obviously, the antikink $-Q_{\lambda}(r)$ is also the
solution with the same energy. Jendrej and Lawrie proved (see Theorem 1 in
[4]) that any finite-energy solution of equation (1) posed on the whole space
tends (modulo sign) either to the vacuum $W=1$ or to an alternating chain of
$N$ rescaled kinks and antikinks111Strictly speaking [4] deals with 2$d$
equivariant wave maps which split into equivariance classes indexed by a
positive integer $k$. The case $k=2$ is essentially equivalent to the 4$d$
equivariant YM. The only qualitative difference is that for wave maps there
are infinitely many topological sectors, while for YM there are only two
sectors (modulo sign). For this reason, in the case of wave maps the chain of
kinks and antikinks in Theorem 1 in [4] need not be alternating.

$1+\sum_{j=1}^{N}(-1)^{N+j}\left(Q_{\lambda_{j}(t)}(r)-1\right)\,.$ (6)

Here $\lambda_{j}(t)$ are continuous positive functions such that for each
$j=1,...,N$

$\frac{\lambda_{j}(t)}{\lambda_{j+1}(t)}\rightarrow
0\quad\text{as}\quad\begin{cases}t\rightarrow\infty\quad\text{(for global-in-
time solutions)}\\\ t\rightarrow T\quad\text{(for blowup at finite time
$T$)},\end{cases}$

where by convention $\lambda_{N+1}(t)=t$ (in the global case) or
$\lambda_{N+1}(t)=T-t$ (in the blowup case), corresponding to the non-existing
self-similar expansion and collapse. We return now to our toy model and make
the soliton resolution conjecture:

######  Conjecture 1.

Any smooth, finite-energy solution $W(t,r)$ of equation (1) subject to the
boundary condition $W(t,1)=0$ tends for $t\rightarrow\infty$ (modulo sign)
either to the half-kink or to the rescaled half-kink plus an alternating chain
of $N$ rescaled kinks and antikinks:

$\begin{cases}\quad Q(r)&\text{if $N=0$,}\\\
-Q_{\mu(t)}(r)+\sum\limits_{j=1}^{N}(-1)^{j+1}Q_{\lambda_{j}(t)}(r)+1&\text{if
$N$ is odd,}\\\ \quad
Q_{\mu(t)}(r)+\sum\limits_{j=1}^{N}(-1)^{j}Q_{\lambda_{j}(t)}(r)&\text{if
$N\geq 2$ is even.}\end{cases}$ (7)

Here $\lambda_{j}(t)$ are continuous positive functions such that for each
$j=1,...,N$

$\lambda_{j}(t)\rightarrow\infty\quad\text{and}\quad\frac{\lambda_{j}(t)}{\lambda_{j+1}(t)}\rightarrow
0\quad\text{as}\quad t\rightarrow\infty,$ (8)

where by convention $\lambda_{N+1}(t)=t$. The function $\mu(t)$ is determined
by the functions $\lambda_{j}(t)$ through the boundary condition $W(t,1)=0$
which implies that $\mu(t)\rightarrow 1$ as $t\rightarrow\infty$.

In the rest of the paper we confirm this conjecture for $N=0$ and $N=1$ and
determine the rate of convergence to the attractors. In addition, we find a
borderline between the basins of attraction using bisection along an
interpolating one-parameter family of initial data. Before presenting the
results of numerical simulations, we wish to put forward a heuristic argument
based on the method of collective coordinates that helps to understand some
aspects of asymptotic dynamics. According to (7), in the case $N=1$ the
attractor has the following form:

$W(t,r)=1-Q_{\mu(t)}(r)+Q_{\lambda(t)}(r),\qquad\mu^{2}(t)=\frac{\lambda^{2}(t)-1}{\lambda^{2}(t)+3},$
(9)

where the formula for $\mu(t)$ follows from the boundary condition $W(t,1)=0$.
Inserting this ansatz into the lagrangian

$L=\frac{1}{2}\int_{1}^{\infty}\left(W_{t}^{2}-W_{r}^{2}-\frac{(1-W^{2})^{2}}{r^{2}}\right)rdr$
(10)

and integrating over $r$, in the limit of large $\lambda$ we get the effective
lagrangian (we retain only the first two leading terms)

$L_{\text{eff}}=\left(\frac{4}{3}-\frac{32}{\lambda^{4}}\right)\,{\dot{\lambda}}^{2}-\left(2-\frac{16}{\lambda^{2}}\right)\,,$
(11)

hence

$\left(\frac{4}{3}-\frac{32}{\lambda^{4}}\right)\,{\dot{\lambda}}^{2}+\left(2-\frac{16}{\lambda^{2}}\right)=E_{\text{eff}}=\text{const}.$
(12)

Thus, the $\lambda$-particle starting at some large $\lambda(0)$ with velocity
$\dot{\lambda}(0)>0$ escapes to infinity if $E_{\text{eff}}\geq 2$, while if
$E_{\text{eff}}<~{}2$ it reaches a turning point in a finite time. The ODE
(12) provides a qualitative picture of the attractive interaction between the
anti-half-kink and the expanding outer kink. This approximation ceases to work
when the outer kink starts shrinking because the PDE solution is no longer
close to the ansatz (9) (see Fig. 4 below). At the quantitative level the
predictions of the effective model should be taken with caution because the
ansatz (9) neglects radiation. In particular, according to Conjecture 1 the
expansion rate $\dot{\lambda}(t)$ in (7) must go to zero as
$t\rightarrow\infty$, whereas the ODE yields $\dot{\lambda}(\infty)>0$ if
$E_{\text{eff}}>2$.

##  3\. Characteristic formulation

We now reformulate our problem as the characteristic initial boundary value
problem. To this end we define new coordinates

$u=t-r,\qquad x=r^{-\frac{1}{2}}.$

Then, on the interval $0<x\leq 1$, the YM potential $\mathsf{w}(u,x)=W(t,r)$
satisfies

$\displaystyle-4x\,\mathsf{w}_{xu}+4\mathsf{w}_{u}=x^{4}\,\mathsf{w}_{xx}+x^{3}\,\mathsf{w}_{x}+8x^{2}\mathsf{w}(1-\mathsf{w}^{2}),$
(13) $\displaystyle\mathsf{w}(u,1)=0,\qquad\mathsf{w}(0,x)=g(x),$ (14)

where the function $g(x)$ is assumed to be smooth and satisfying the finite
energy condition $g(0)=1$. For such data, the solution $\mathsf{w}(u,x)$
remains smooth for all future times $u>0$. Moreover, the results by Chruściel
and collaborators [7, 8] imply222We are grateful to Piotr Chruściel and Roger
Tagné Wafo for checking that the hypotheses of theorems on propagation of
polyhomogeneity in [7, 8] hold for our equation. This, together with the
absence of $\log{x}$ terms in the formal polyhomogeneous series, shows that
the solution is smooth in $x$.  that the following asymptotic expansion holds
near $x=0$:

$\mathsf{w}(u,x)=1+\sum_{n=1}c_{n}(u)x^{n}\,.$ (15)

Inserting this expansion into equation (13) and equating the coefficients of
the same powers of $x$, we obtain an infinite system of ordinary differential
equations for the coefficients $c_{n}(u)$. This system can be solved
recursively one-by-one starting from the radiation coefficient $c_{1}(u)$
which is free. For $n=2$ we get $\dot{c}_{2}(u)=0$, hence the coefficient
$c_{2}$ is constant (so called Newman-Penrose constant). For large $n$ the
nonzero coefficients $c_{n}(u)$ grow polynomially for $u\rightarrow\infty$
which is a reflection of the well-known fact that the expansion (15) is not
uniform; see e.g. [9]. Multiplying equation (13) by $x^{-3}\mathsf{w}_{u}$ we
get the local conservation law

$\partial_{u}\left(\frac{x}{2}\mathsf{w}_{x}^{2}+\frac{2}{x}(1-\mathsf{w}^{2})^{2}\right)=\partial_{x}\left(\frac{2}{x^{2}}\,\mathsf{w}_{u}^{2}+x\,\mathsf{w}_{u}\,\mathsf{w}_{x}\right).$
(16)

Integrating this over $x$ and using (15), we obtain the energy loss formula

$\frac{d\mathcal{E}}{du}=-\dot{c}_{1}^{2}(u)\,,$ (17)

where

$\mathcal{E}[\mathsf{w}]:=\int_{0}^{1}\left(\frac{1}{4}\,\mathsf{w}_{x}^{2}+\frac{1}{x^{2}}(1-\mathsf{w}^{2})^{2}\right)\,x\,dx$
(18)

is the Bondi-type energy (hereafter just called energy). Note that
$\mathcal{E}[\mathsf{w}]$ is equal to the potential part of the total
conserved energy $E[W]$. In terms of $x$ the half-kink reads

$q(x):=Q\left(\frac{1}{x^{2}}\right)=\frac{1-x^{4}}{1+x^{4}}.$ (19)

It is the global minimizer of $\mathcal{E}[\mathsf{w}]$ as follows from the
Bogomolnyi inequality

$\displaystyle\mathcal{E}[\mathsf{w}]=\int_{0}^{1}\left[\frac{1}{2}\mathsf{w}_{x}+\frac{1}{x}(1-\mathsf{w}^{2})\right]^{2}\,x\,dx-\int_{0}^{1}\partial_{x}(\mathsf{w}-\frac{1}{3}\mathsf{w}^{3})\,dx\geq\frac{2}{3},$
(20)

which is saturated on $\mathsf{w}=q(x)$, i.e. $\mathcal{E}[q]=\frac{2}{3}$.
Since $\mathcal{E}[\mathsf{w}]$ is non-increasing and bounded below, there
exists a limit

$\mathcal{E}_{\infty}=\lim_{u\rightarrow\infty}\mathcal{E}(\mathsf{w}(u))\geq\mathcal{E}(q)=\frac{2}{3}\,.$

According to Conjecture 1 the endstate (which clearly must be non-radiative)
has the form (7), hence the final energy is a sum of the energies of the half-
kink and $N$ kinks/antikinks

$\mathcal{E}_{\infty}=\frac{2}{3}+\frac{4}{3}N.$ (21)

In section 5 we describe the relaxation to the $N=0$ and $N=1$ attractors. In
our numerical simulations we have not observed $N\geq 2$ attractors which
suggests that, if they exist, they are nongeneric.

##  4\. Linearized dynamics near the half-kink

Let $\mathsf{w}=q(x)+xf(u,x)$. Substituting this into equation (13) we obtain

$f_{ux}+\frac{1}{4}\partial_{x}(x^{3}f_{x})-U(x)f-6q(x)x^{2}f^{2}-2x^{3}f^{3}=0,\qquad
f(u,1)=0,$ (22)

where

$U(x)=\frac{x(15-66x^{4}+5x^{8})}{4(1+x^{4})^{2}}=\frac{15}{4}x+\mathcal{O}(x^{5}).$
(23)

Dropping the nonlinear terms and the $\mathcal{O}(x^{5})$ term in the
potential, we get the linear equation (corresponding to linearization around
$\mathsf{w}=1$ rather than $q$)

$f_{ux}+\frac{1}{4}\partial_{x}(x^{3}f_{x})-\frac{15}{4}xf=0.$ (24)

This equation has an explicit solution

$f_{0}(u,x)=u^{-\frac{5}{2}}(ux^{2}+2)^{-\frac{5}{2}}.$ (25)

General solutions of equation (24) for initial data with the vanishing NP
constant behave similarly to $f_{0}(u,x)$ for $u\rightarrow\infty$, i.e. they
decay as $u^{-5}$ in the interior ($x>0$) and $u^{-5/2}$ at future null
infinity ($x=0$). This can be shown directly, or by defining
$F(t,r)=x^{5}f(u,x)$ and rewriting equation (24) in terms of the original
coordinates $(t,r)$. Then, $F(t,r)$ satisfies the free radial wave equation in
$6+1$ dimensions $F_{tt}-F_{rr}-\frac{5}{r}F_{r}=0$, for which the late-time
pointwise decay $F(t,r)\sim(t-r)^{-5/2}(t+r)^{-5/2}$ is well known [10]. The
rate of decay of linear perturbations about $q(x)$ is the same because the
term $\mathcal{O}(x^{5})$ in the potential (23) is asymptotically negligible.
The numerical verification of this claim is depicted in Fig. 1 where we also
show that solutions with nonzero NP constant exhibit slower decay.

Figure 1. Linear evolution about the half-kink $q(x)$ for sample initial data
with zero and nonzero NP constant, respectively: $f(0,x)=\cos^{2}\\!\left(\pi
x/2\right)$ (top row) and $f(0,x)=\cos^{2}\\!\left(\pi x/2\right)+(1-x)x$
(bottom row). On the left, we plot $\log|f(u,x_{i})|$ at different grid points
(solid lines labelled by the grid point number $i$, where $x_{256}=0$; see
section 5 for the details) together with the theoretical decay rates
(dashed/dotted lines). On the right, we plot the local power index
$p(u,x_{j}):=u\,\partial_{u}f(u,x_{i})/f(u,x_{i})$ evaluated at different grid
points.

Having determined the late-time linear tail, now we turn to the computation of
quasinormal modes. They govern the relaxation to the half-kink for
intermediate times before the tail is uncovered. It is convenient to rewrite
the linear part of equation (22) in terms of $y=x^{2}$. Substituting

$f(u,x)=e^{su}\,v(y),$ (26)

we get the eigenvalue problem

$2sv^{\prime}+(y^{2}v^{\prime})^{\prime}-V(y)v=0,\qquad
V(y)=\frac{15}{4}-\frac{48+5y^{2}}{(1+y^{2})^{2}}\,y^{2}$ (27)

with $v(1)=0$. Following Leaver [11] we seek solutions of (27) in terms of the
power series

$v(y)=\sum_{n\geq 1}a_{n}(1-y)^{n},\qquad a_{1}=1.$ (28)

Since the nearest singularity from $y=1$ is located at $y=0$, this power
series is absolutely convergent for $y\in(0,1]$. The eigenvalues $s_{n}$,
called quasinormal frequencies, are selected by the condition that the power
series is absolutely convergent at $y=0$; as follows from (26), the
corresponding solutions $f_{n}(u,y)$, called the quasinormal modes, are purely
outgoing333An alternative way of selecting the outgoing solution by a certain
Gevrey-class regularity condition has been recently proposed by Gajic and
Warnick [13]; see also [14]. For a very interesting discussion of
hyperboloidal approach to quasinormal modes and Leaver’s method we refer the
reader to [15, 16].. Inserting (28) into equation (27) we get a 7-term
recurrence relation. Among its six linearly independent solutions, four
solutions decay as $a_{n}\sim 2^{-\frac{n}{2}}$ for $n\rightarrow\infty$,
hence they do not affect the convergence properties of the series at $y=0$.
Using the method of successive approximations [12] one can show that the
remaining two solutions have the following asymptotic expansions

$a^{(+)}_{n}\sim
n^{-\frac{3}{4}}e^{\sqrt{8sn}}\sum_{k=0}c^{(+)}_{k}n^{-\frac{k}{2}},\qquad
a^{(-)}_{n}\sim
n^{-\frac{3}{4}}e^{-\sqrt{8sn}}\sum_{k=0}c^{(-)}_{k}n^{-\frac{k}{2}},$ (29)

where the coefficients $c^{(\pm)}_{k}$ can be determined successively by
plugging the expansions (29) into the recurrence relation. We conclude that
for $n\rightarrow\infty$

$a_{n}=C_{+}(s)a^{(+)}_{n}+C_{-}(s)a^{(-)}_{n}+\mathcal{O}(2^{-\frac{n}{2}}).$
(30)

For $|\arg(s)|<\pi$, the series $\sum|a^{(+)}_{n}|$ diverges while the series
$\sum|a^{(-)}_{n}|$ converges, therefore the quasinormal frequencies are given
by the roots of the coefficient $C_{+}(s)$. There are several alternative ways
to find these roots. The most frequently used is Leaver’s method of continued
fractions [11]. It is stable and accurate but tedious in the case at hand
because the recurrence relation must first be reduced to three terms by
Gaussian elimination. Employing this method we found exactly one pair of
complex conjugate frequencies $s\approx-0.364322\pm 0.476858i$. To verify this
result, we have reproduced it by two different methods: a brute force
evaluation of a dominant solution by forward recurrence and an algebraic
method introduced in [17]. We skip the details of these straightforward but
dull computations. We confirmed the above perturbative analysis by the direct
numerical integration of the linearized equation; see Fig. 2.

Figure 2. The solution depicted in blue on the top left plot in Fig. 1 is
shown here for early times when the relaxation to the half-kink is governed by
the quasinormal mode. The fit (dashed red line) of a superposition of an
exponentially damped oscillation and the polynomial tail to the numerical data
(black line) gives the quasinormal frequency $s=-0.364271\pm 0.476856\,i$, in
very good agreement with the perturbative calculation. The pure quasinormal
ringdown is depicted for reference by the dashed blue line.

######  Remark.

It is instructive to compare the above computation of the quasinormal modes
for the half-kink with an analogous computation for the vacuum solution
$\mathsf{w}=1$ of equation (13) with the boundary condition
$\mathsf{w}(u,1)=1$ (as mentioned above, this is equivalent to the free wave
equation in $6+1$ dimensions). For the ansatz
$\mathsf{w}=1+\sqrt{y}\,e^{su}v(y)$ (where $y=x^{2}$), we obtain the same
eigenvalue equation as (27) but with the potential $V=15/4$. Repeating the
above analysis, we get a three-term recurrence relation having two linearly
independent solutions $a^{(\pm)}_{n}$ with the asymptotic behavior (29), hence
as before the quasinormal frequencies are given by the roots of the
coefficient $C_{+}(s)$ of the dominant solution. We remark that in this case
the analysis based on the recurrence relation is purely academic because the
eigenvalue equation can be solved exactly and the outgoing solution is given
by $v_{\text{out}}=y^{-\frac{1}{2}}e^{s/y}K_{2}(s/y)$, where $K_{2}(z)$ is the
modified Bessel function of the second kind. Thus, the quantization condition
for the quasinormal modes is $K_{2}(s)=0$, which has exactly one pair of
complex conjugate zeros on the principal branch $s=-1.281373\pm 0.4294849i$
[18]. We verified that the roots of $C_{+}(s)$ are the same, which provides a
reassuring benchmark test for Leaver’s method.

##  5\. Numerical results

In this section we corroborate Conjecture 1 with direct numerical simulations
of the initial boundary value problem (13)-(14). As in section 4, we write
$\mathsf{w}=q(x)+xf(u,x)$ and then solve equation (22) numerically using the
method of lines. To this end, we first discretize equation (22) in space using
the pseudospectral approach. For numerical convenience, we rescale the spatial
domain to the interval $[-1,1]$ using $z=2x-1$ and work with function values
$\\{f_{j}(u)\equiv f(u,z_{j})\\}$ evaluated at $K$ Chebyshev points of the
second kind $\\{z_{j}=\cos\left(\frac{(j-1)\pi}{K-1}\right)\\}$, $1\leq j\leq
K$. Spatial derivatives $\partial_{z}$ and $\partial_{z}^{2}$ are replaced by
the corresponding spectral differentiation matrices $D^{(1)}_{K}$ and
$D^{(2)}_{K}$ [21] and then both the derivatives and nonlinear terms are
evaluated using the grid function $\\{f_{i}\\}$. The resulting semi-discrete
system takes the following schematic form

$\partial_{u}f_{1}=0,\qquad\partial_{u}(D^{(1)}_{K}f)_{j}=F_{j}(D^{(2)}_{K}f,D^{(1)}_{K}f,f,z)\,,\quad
j=2,\ldots,K,$ (31)

where the boundary condition $f_{1}=0$ replaces the equation at the grid point
$z_{1}=1$. We bring this system to an explicit form by solving the linear
system for $u$ derivatives of $\\{f_{j}\\}$. This requires inverting the
operator $D^{(1)}_{N}$ with the first row replaced by a condition
$\partial_{u}f_{1}=0$. Note that this linear operator is invertible and non-
degenerate. The solution uses the LU decomposition of the resulting matrix.
For efficient time integration we use an implicit scheme. We employ the BDF
method (variable-order, variable-coefficient, in fixed-leading-coefficient
form) which, for the best performance, we limit to the second order (higher-
order methods struggled to find the optimal step size/order, probably due to
the stiffness of the equation). We used the IDA code [22], as available in
Wolfram Mathematica [23], in which we set the error tolerances to very
conservative values (typically $10^{-13}$) so that the spatial resolution
determines the errors in the numerical solution. Tests of the final algorithm
show the spectral (exponential) convergence with increasing $K$. Using the
above method we have simulated the evolution of various initial data. Here we
illustrate the results for a sample one-parameter family

$\mathsf{w}(0,x)=1+bx^{4}-(1+b)x^{6}\,,$ (32)

where $b$ is a free parameter. For this data the NP constant is equal to zero.
To see how a nonzero NP constant affects the dynamics we look in parallel at
the evolution of initial data (32) with the additional term $x(1-x)$. The
energy of initial data (32) attains the minimum value $0.6672$ at
$b\approx-2.1022$. For $-7.7295\lesssim b\lesssim 2.5933$ we have
$\mathcal{E}<2$, hence according to (21) the half-kink $q$ is the only
possible attractor (the $N=0$ case in our terminology). In agreement with
this, we observe rapid convergence to $q$ through a short ringdown and then
the late-time tail; see Fig. 3. Notice that the linear decay rates determined
in the previous section, namely $u^{-5}$ in the interior ($x>~{}0$) and
$u^{-5/2}$ along the future null infinity ($x=0$), are not propagated by the
nonlinear flow444Interestingly, if the NP constant is nonzero, then the linear
and nonlinear tails decay at the same rate: $u^{-3}$ for $x>0$ and $u^{-1/2}$
for $x=0$; see the bottom rows in Fig. 1 and 3. for which the decay rates are
slower by one power of $u$. Most important for us is the decay rate for the
radiation coefficient $c_{1}(u)\sim u^{-3/2}$, however establishing this fact
rigorously is a task that goes beyond the scope of this paper. We remark that
similar nonlinear tails (but only in the interior) have been studied in the
literature for semilinear wave equations in high even spatial dimensions (in
particular, for the quadratic wave equation in $6+1$ dimensions which is
relevant in our context); see [19] and [20].

Figure 3. Pointwise decay to the half-kink for the initial data (32) with
$b=1$ (top row) and for the same data with an extra term $x(1-~{}x)$ which
generates a nonzero NP constant (bottom row). Compare with the analogous plots
of the linear decay shown in Fig. 1.

Next, we consider initial data with energy greater than $2$. For moderate
positive values of $b$ the solution again quickly converges to the half-kink,
however for larger values of $b$ we observe formation of the superposition of
the anti-half-kink and the expanding kink (the $N=1$ attractor in our
terminology). We find that the transition between these two scenarios occurs
at $b_{0}\approx 12.458288341217909$. For marginally subcritical solutions
(i.e. for $b=b_{0}-\varepsilon$ with small positive $\varepsilon$) a
superposition of the anti-half-kink and the expanding kink appears for
intermediate times but at a later time the expansion stops, the kink starts
shrinking and is quickly annihilated. In this process the energy of the kink
is rapidly radiated away and the solution settles down to the half-kink. This
behavior is shown in Figs. 4 and 5 where we plot the snapshots of marginally
subcritical and supercritical solutions and the corresponding energies,
respectively.

Figure 4. Snapshots from the evolution for the initial data (32) with
marginally subcritical (blue) and supercritcal (red) values of $b$ near
$b_{0}$. The half-kink and anti-half-kink are plotted with dashed lines. For
the subcritical evolution the turning point is at $u_{R}\approx 1000$.  Figure
5. The energies of solutions from Fig. 4. In the inset we zoom into the
intermediate phase of subcritical evolution when the energy of the outer kink
is radiated away.

For a more quantitative description of the expanding phase of subcritical
solutions, let $x_{0}(u)$ be the zero of the solution $\mathsf{w}(u,x)$ and
$u_{R}$ be the return time when the expansion stops. We find that
$u_{R}\sim\varepsilon^{-1}$ and $x_{0}(u_{R})\sim\varepsilon^{1/4}$; see Fig.
6. Translating these scaling relations to the original variables and comparing
with the ansatz (9) we get $t_{R}\sim\varepsilon^{-1}$ and
$\lambda(t_{R})\sim\varepsilon^{-1/2}$. This is in agreement with the ODE
approximation (12) which predicts that for marginally subthreshold effective
energies we have $\lambda(t_{R})\sim(2-E_{\text{eff}})^{-1/2}$.

Figure 6. Left panel: for marginally subcritical solutions we plot the zero of
$\mathsf{w}(u,x)$, denoted by $x_{0}(u)$, for several values of
$\log(b_{0}-b)$ (depicted by different colors). For almost critical data the
fit gives $x_{0}\sim u^{-0.2494}$. This is in accord with the ODE
approximation (12) which gives $\lambda(t)\sim t^{1/2}$ for motion on the
separatrix with $E_{\text{eff}}=2$. Right panel: $x_{0}(u)$ at the return time
$u_{R}$ as the function of $b_{0}-b$. The fit gives
$x_{0}(u_{R})\sim(b_{0}-b)^{0.2495}$ (red line).

Increasing $b$ we find that above $b_{1}\approx 47.90418049175238$ the
solution again settles down on the half-kink after a few rapid nonlinear
oscillations. A similar transition occurs below
$b_{-1}\approx-53.94479194728988$. As $|b|$ grows further the endstate keeps
flipping back and forth between the $N=0$ and $N=1$ attractors. We conjecture
that there are infinitely many critical values $b_{n}$ ($n\in\mathbb{Z}$) at
which the curve of initial data (32) intersects the $N=0$ and $N=1$ basins of
attraction. In the case $N=1$ of Conjecture 1, it remains to verify that the
speed of expansion of the outer kink goes asymptotically to zero, i.e.
$\frac{\lambda(t)}{t}\rightarrow 0$ for $t\rightarrow\infty$. This is shown in
Fig. 7. Unfortunately, we are not in position to say more about the dynamics
of $\lambda(t)$. The computation of the precise asymptotic behavior of
$\lambda(t)$, which takes into account the loss of energy by radiation, is a
challenging open problem that we leave to future work555See [24, 25] for the
derivation of the modulation equation for $\lambda(t)$ for the blowup
solutions of the $4+1$ YM equation in the whole space..

Figure 7. The speed of expansion of the outer kink for supercritical evolution
with $b=20$. To compute the function $\lambda(t)$, we first find the zero
$x_{0}(u)$ of the solution $\mathsf{w}(u,x)$, then translate the result to the
variables $(t,r)$, and finally compare it with the ansatz (9).

_Acknowledgement._ PB wishes to thank Piotr Chruściel, Peter Hintz and Jacek
Jendrej for useful comments. Special thanks are due to Arthur Wasserman for
very helpful suggestions. The work of PB and BC was supported by the National
Science Centre grant no. 2017/26/A/ST2/00530. MM acknowledges the support of
the Austrian Science Fund (FWF), Project P 29517-N27 and the START-Project
Y963.

## References

  * [1] T. Tao, _Why are solitons stable?,_ Bull. Amer. Math. Soc. 46, 1-33 (2009)
  * [2] C.E. Kenig, _Asymptotic simplification for solutions of the energy critical nonlinear wave equation,_ J. Math. Phys. 62, 011502 (2021)
  * [3] T. Duyckaerts, C. Kenig, Y. Martel, F. Merle, _Soliton resolution for critical co-rotational wave maps and radial cubic wave equation,_ arXiv:2103.01293
  * [4] J. Jendrej, A. Lawrie, _Soliton resolution for equivariant wave maps_ , arXiv:2106.10738
  * [5] P. Bizoń, T. Chmaj, M. Maliborski, _Equivariant wave maps exterior to a ball_ ,
Nonlinearity 25, 1299-1309(2012)

  * [6] C.E. Kenig, A. Lawrie, W. Schlag, _Relaxation of wave maps exterior to a ball to harmonic maps for all data,_ Geom. Funct. Anal. 24, 610-647 (2014)
  * [7] P.T. Chruściel, S. Łȩski, _Polyhomogeneous solutions of nonlinear wave equations without corner conditions,_ J. Hyperbolic Diff. Equations 3, 81-141 (2006)
  * [8] P.T. Chruściel, R. Tagné Wafo, _Solutions of quasi-linear wave equations polyhomogeneous at null infinity in high dimensions,_ J. Hyperbolic Diff. Equ. 8, 269-346 (2011)
  * [9] P. Bizoń, H. Friedrich, _A remark about wave equations on the extreme Reissner-Nordström black hole exterior_ , Class. Quantum Grav. 30, 065001 (2013)
  * [10] L. Hörmander, _Lectures on Nonlinear Hyperbolic Differential Equations,_ Springer 1997.
  * [11] E. Leaver, _Quasinormal modes of Reissner-Nordström black holes,_ Phys. Rev D41, 2986–2997 (1990)
  * [12] R. Wong, H. Li, _Asymptotic expansions for second-order linear difference equations,_
J. Comput. Appl. Math. 41, 65-94 (1992)

  * [13] D. Gajic, C. Warnick, _A model problem for quasinormal ringdown on asymptotically flat or extremal black holes,_ J. Math. Phys. 61, 12501 (2020)
  * [14] J. Galkowski, M. Zworski, _Outgoing solutions via Gevrey-2 properties,_ arXiv:2004.07868
  * [15] M. Ansorg, R.P. Macedo, _Spectral decomposition of black-hole perturbations on hyperboloidal slices,_ Phys. Rev. D 93, 124016 (2016)
  * [16] R.P. Macedo, J.J. Jaramillo, M. Ansorg, _Hyperboloidal slicing approach to quasi-normal mode expansions: the Reissner-Nordström case,_ Phys. Rev. D 98, 124005 (2018)
  * [17] P. Bizoń, M. Maliborski, _Dynamics at the threshold for blowup for supercritical wave equations outside a ball,_ Nonlinearity 33, 3195-3205 (2020)
  * [18] R. Parnes, _Complex zeros of the modified Bessel function $K_{n}(Z)$,_ Mathematics of Computation 26, 949-953 (1972)
  * [19] R. Agemi, K. Kubota, H. Takamura, _On certain integral equations related to nonlinear wave equations,_ Hokkaido Math. J. 23, 241-276 (1994)
  * [20] P. Hintz, A. Vasy, _Semilinear wave equations on asymptotically de Sitter, Kerr-de Sitter and Minkowski spacetimes,_ Analysis & PDE 8, 1807-1890 (2015)
  * [21] L.N. Trefethen, _Spectral Methods in MATLAB_ , SIAM, Philadelphia, 2000
  * [22] A. Hindmarsh, A. Taylor. _User Documentation for IDA: A Differential-Algebraic Equation Solver for Sequential and Parallel Computers_ Lawrence Livermore National Laboratory report, UCRL-MA-136910, 1999
  * [23] Wolfram Research, Inc., _Mathematica 12.3.1_ , Champaign, Illinios, 2021
  * [24] P. Bizoń, Yu. Ovchinnikov, M.I. Sigal, _Collapse of an instanton_ , Nonlinearity 17, 1179-1191 (2004)
  * [25] P. Raphaël, I. Rodnianski, _Stable blow up dynamics for the critical co-rotational wave maps and equivariant Yang-Mills problems,_ Publ. math. IHES 115, 1-122 (2012)