# Implications of NANOGrav results and UV freeze-in in a fast-expanding
Universe

Basabendu Barman<EMAIL_ADDRESS>Department of Physics, IIT Guwahati,
Guwahati-781039, India  Amit Dutta Banik<EMAIL_ADDRESS>Key
Laboratory of Quark and Lepton Physics (MoE) and Institute of Particle
Physics, Central China Normal University, Wuhan 430079, China  Avik Paul
<EMAIL_ADDRESS>Astroparticle Physics and Cosmology Division, Saha
Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Kolkata 700064, India

###### Abstract

Recent pulsar timing data reported by the NANOGrav collaboration indicates the
existence of a stochastic gravitational wave (GW) background at a frequency
$f\sim 10^{-8}\leavevmode\nobreak\ \rm Hz$. We show that a dark sector
consisting of a Standard Model (SM) gauge singlet fermion $\chi$ and a singlet
scalar $\phi$, both charged under a $Z_{4}$ symmetry, is capable of generating
such a low frequency GW via strong first order phase transition (SFOPT)
through the modification of the standard cosmological history, where we assume
faster-than-usual expansion at pre-BBN times driven by a new cosmological
species $\varphi$ whose energy density red-shifts with the scale factor as
$\rho_{\varphi}\propto a^{-\left(4+n\right)}$. Depending on the choice of the
fast expansion parameters, reheat temperature and effective scale of the
theory, it is also possible to address correct dark matter (DM) relic
abundance via freeze-in. We show that a successful first order phase
transition explaining NANOGrav results together with PLANCK observed DM
abundance put bound on the fast expansion parameters requiring $n\lesssim 4$
to explain both.

##  I Introduction

Recently, the North American Nanohertz Observatory for Gravitational Waves
(NANOGrav) has reported a possible detection of a Stochastic gravitational
wave background (SGWB) with their 12.5 year data set of the pulsar timing
array Arzoumanian  _et al._ (2020) that may be interpreted as a gravitational
wave (GW) signal with amplitude of $\sim\mathcal{O}\left(10^{-15}\right)$ at a
reference frequency of $f_{\text{yr}}=1\text{yr}^{-1}$. The NANOGrav
collaboration has found no evidence of monopolar or dipolar correlations which
may arise from reference clock or solar-system ephemeris systematics, they
also claim that there is no statistically significant evidence that this
process has quadrupolar spatial correlations which can be a smoking gun for GW
background. However, several attempts were proposed, so as to interpret the
NANOGrav observations which include supermassive black hole merger events
Arzoumanian  _et al._ (2020); Sesana  _et al._ (2004), cosmic strings Ellis
and Lewicki (2020); Blasi  _et al._ (2020); Chigusa  _et al._ (2020), magneto-
hydrodynamic (MHD) turbulence at a first-order cosmological QCD phase
transition Neronov  _et al._ (2020) or phase transition from a hidden sector
Nakai  _et al._ (2020). As shown in Addazi  _et al._ (2020) the NANOGrav data
can be well explained as a GW signal from first order phase transition (FOPT)
around the warm dark matter (DM) physics scale that lies in the keV-MeV range,
well below the electroweak scale. In Bhattacharya  _et al._ (2020a) it was
shown that NANOGrav observation can be explained by the first order GWs in the
nonstandard thermal history with an early matter dominated (EMD) era.

The expansion rate of the universe is parametrized by the Hubble parameter
which is directly related to the total energy density of the universe through
standard Friedman equations Kolb and Turner (1990). In the standard case, it
is assumed that the universe was radiation dominated before Big Bang
Nucleosynthesis (BBN). However, there are no fundamental reasons to assume
that the universe was radiation-dominated prior to BBN at $t\sim
1\leavevmode\nobreak\ \rm sec$. The absence of any direct evidence that the
universe was indeed only radiation dominated in the pre-BBN era compels us to
explore a modified cosmological scenario where the expansion rate of the
universe naturally alters from what it is in case of the standard scenario.
The influence of the presence of another species that drives the expansion of
the universe faster than the standard cosmological history on DM yield has
been studied in detail both in the context of Weakly Interacting Massive
Particle (WIMP) D’Eramo  _et al._ (2017) and Feebly Interacting Massive
Particle (FIMP) D’Eramo  _et al._ (2018). In case of DM freeze-in it has been
shown that in a fast expanding universe the DM production via decay and/or
annihilation of the bath particles is dramatically suppressed that demands a
larger coupling between the dark and the visible sector to satisfy the observe
relic abundance. This sheds some optimistic prospects on detecting DM produced
via freeze-in which otherwise is extremely difficult to detect.

Motivated from all these in this work we propose a simple framework where we
assume the presence of another species $\varphi$ whose energy density red-
shifts with the scale factor as $\rho_{\varphi}\propto a^{-\left(4+n\right)}$.
For $n>0$ the $\varphi$ energy density dominates over radiation at early
enough times. However, the equality between $\rho_{\varphi}$ and
$\rho_{\text{rad}}$ happens at a temperature $T_{R}\gtrsim T_{\text{BBN}}$.
This, as one can understand, alters the Hubble parameter leading to a faster-
than-usual-expansion of the universe in the pre-BBN era. In such a modified
cosmological background we consider a simple dark sector made up of a Standard
Model (SM) gauge singlet fermion $\chi$ and a singlet scalar $\phi$ both
charged under a $Z_{4}$ symmetry. We assume $\chi$ is lighter than $\phi$ and
hence it serves as a potential DM candidate. The lack of any renormalizable
interaction between the dark sector and the SM motivates the DM production via
Ultra-violate (UV) freeze-in Elahi  _et al._ (2015); Chen and Kang (2018);
Bernal  _et al._ (2019). The requirement of satisfying PLANCK observed Ade
_et al._ (2014); Aghanim  _et al._ (2018) relic abundance constraints the
reheat temperature, effective scale of the theory and also the fast expansion
parameters $n$ and $T_{R}$. The singlet scalar, on the other hand, triggers a
strong first order electroweak phase transition (SFOPT) that generates
detectable GW signal. We show, in a fast expanding universe such GW
signals111GW signal in a modified cosmological scenario have also been studied
in details in Bernal and Hajkarim (2019); Bhattacharya  _et al._ (2020b);
Bernal  _et al._ (2020) generated from electroweak phase transition can
explain the NANOGrav 12.5 yrs data within 2$\sigma$ confidence. The
requirement of producing right DM abundance together with GW signal satisfying
NANOGrav results put bound on the fast expansion parameters, typically
requiring $n\lesssim 4$ and
$T_{R}\sim\mathcal{O}\left(100\right)\leavevmode\nobreak\ \rm MeV$.

The paper is organized as follows: in Sec. II we describe the modification in
the Hubble parameter due to fast expansion; the effect of modified Hubble
parameter on DM yield in the context of a toy model is detailed in Sec. III;
in Sec. IV we discuss the generation of detectable GW signal from strong first
order phase transition that is capable of explaining the NANOGrav 12.5 yrs
data; we then briefly touch upon a possible UV completion for this toy model
in Sec. V and finally we conclude in Sec. VI.

##  II A Fast expanding universe

As mentioned in Sec. I, we assume the universe before BBN has two different
species: radiation and some other species $\varphi$ with energy densities
$\rho_{\rm rad}$ and $\rho_{\varphi}$ respectively. In presence of a new
species ($\varphi$) along with the radiation field, the total energy budget of
the universe is $\rho=\rho_{\text{rad}}+\rho_{\varphi}$. One can always
express the energy density of the radiation component as function of
temperature $T$ Kolb and Turner (1990)

$\displaystyle\rho_{\text{rad}}(T)=\frac{\pi^{2}}{30}g_{*}(T)T^{4}$ (1)

with $g_{*}(T)$ being the effective number of relativistic degrees of freedom
at temperature $T$. In the absence of entropy production per comoving volume
i.e. $sa^{3}=$ const., radiation energy density redshifts as
$\rho_{\text{rad}}(t)\propto a(t)^{-4}$. In case of a rapid expansion of the
universe the energy density of $\varphi$ field is expected to be redshifted
faster than the radiation. Accordingly, one can assume $\rho_{\varphi}\propto
a(t)^{-(4+n)}$. Thus $n>0$ implies that $\varphi$ energy density dominates
over radiation at early enough times.

Now, the entropy density (per comoving volume) of the universe is expressed as

$s(T)=\frac{2\pi^{2}}{45}g_{*s}(T)T^{3}$

where $g_{*s}\left(T\right)$ is the effective relativistic degrees of freedom.
A general form of $\rho_{\varphi}$ can then be constructed using the entropy
conservation $g_{\star}\left(T\right)^{1/3}aT=\text{constant}$ in a co-moving
frame as:

$\displaystyle\rho_{\varphi}(T)=\rho_{\varphi}(T_{R})\left(\frac{g_{*s}(T)}{g_{*s}(T_{R})}\right)^{(4+n)/3}\left(\frac{T}{T_{R}}\right)^{(4+n)}.$
(2)

The temperature $T_{R}$ is an unknown variable and can be considered as the
point of equality where the two fluids have equal energy densities:
$\rho_{\varphi}(T_{R})=\rho_{\text{rad}}(T_{R})$. Using this criteria, the
total energy density at any temperature $T$ reads D’Eramo  _et al._ (2017,
2018)

$\displaystyle\rho(T)$ $\displaystyle=\rho_{rad}(T)+\rho_{\varphi}(T)$ (3)
$\displaystyle=\rho_{rad}(T)\left[1+\frac{g_{*}(T_{R})}{g_{*}(T)}\left(\frac{g_{*s}(T)}{g_{*s}(T_{R})}\right)^{(4+n)/3}\left(\frac{T}{T_{R}}\right)^{n}\right]$
(4)

Now, from standard Friedman equation one can write the Hubble parameter in
terms of the energy density as

$\displaystyle\mathcal{H}=\frac{\sqrt{\rho}}{\sqrt{3}M_{\text{Pl}}},$ (5)

with $M_{\text{Pl}}=\left(8\pi G\right)^{-1/2}=2.4\times 10^{18}$ GeV being
the reduced Planck mass. At temperature higher than $T_{R}$ with the condition
$g_{*}(T)=\bar{g}_{*}$ which can be considered to be some constant, the Hubble
rate can approximately be recasted as D’Eramo  _et al._ (2017)

$\displaystyle\begin{split}&\mathcal{H}(T)=\mathcal{H}_{R}\left(T\right)\Biggl{[}1+\Biggl{(}\frac{T}{T_{R}}\Biggr{)}^{n}\Biggr{]}^{1/2}\\\
&\approx\frac{\pi\bar{g}_{*}^{1/2}}{3\sqrt{10}}\frac{T^{2}}{M_{\text{Pl}}}\left(\frac{T}{T_{R}}\right)^{n/2},\leavevmode\nobreak\
\leavevmode\nobreak\ (T\gg T_{R}),\end{split}$ (6)

where $H_{R}\left(T\right)$ is the Hubble parameter in the standard radiation
dominated universe. In case of SM $\bar{g}_{*}\equiv g_{*}\text{(SM)}=106.75$.
It is important to note from Eq. (6) that the expansion rate is larger than
what it is supposed to be in the standard cosmological background for
$T>T_{R}$ and $n>0$. It is noteworthy that $T_{R}$ can not be too small such
that it changes the standard BBN history. For a certain value of $n$, BBN
constraints provide a lower limit on $T_{R}$ D’Eramo  _et al._ (2017):

$\displaystyle T_{R}\gtrsim\left(15.4\right)^{1/n}\leavevmode\nobreak\
\text{MeV}.$ (7)

As one can see, for larger $n$ this constraint becomes more and more lose.
Before moving on to the next section we would like to comment that in order to
get values larger than $n=2$, a negative scalar potential is to be considered.
A specific structure of $n>2$ potential can be found in Ratra and Peebles
(1988); D’Eramo  _et al._ (2017) which is asymptotically free.

##  III Dark matter yield via Freeze-in

In this section we will discuss the modification of standard Boltzmann
equation Kolb and Turner (1990) governing the DM number density due to the
effect of fast expansion encoded in the Hubble parameter $\mathcal{H}$ (Eq.
(6)). The key for freeze-in Hall  _et al._ (2010); Bernal  _et al._ (2017) DM
production is to assume that DM was not present in the early universe. The DM
is then produced via annihilation and/or decay of the particles in the thermal
bath. Due to the extremely feeble coupling of the DM with the visible sector
particles the DM never really enters into thermal equilibrium. Now, the
Boltzmann equation (BEQ) for DM production via annihilation for $2\to 2$
processes of the form $a,b\to 1,2$ is then given by Edsjo and Gondolo (1997):

$\displaystyle\dot{n_{1}}+3\mathcal{H}n_{1}$ $\displaystyle=\int
d\Pi_{1}d\Pi_{2}d\Pi_{a}d\Pi_{b}\left(2\pi\right)^{4}\delta^{4}\left(p_{1}+p_{2}-p_{a}-p_{b}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \left|\mathcal{M}\right|_{a,b\to 1,2}^{2}f_{a}f_{b},$ (8)

where $d\Pi_{j}=\frac{d^{3}p_{j}}{2E_{j}\left(2\pi\right)^{3}}$ are Lorentz
invariant phase space elements, and $f_{i}$ is the phase space density of
particle $i$ with corresponding number density being:

$\displaystyle n_{i}=\frac{g_{i}}{\left(2\pi\right)^{3}}\int d^{3}pf_{i},$ (9)

with $g_{i}$ is the internal degrees of freedom (DOF). Here we have made two
crucial assumptions: the initial abundance for particle $1$ is negligible such
that $f_{1}\approx 0$ and we also neglect Pauli-blocking/stimulated emission
effects, i.e. approximating $1\pm f_{i}\approx 1$. The BEQ in Eq. (8) can be
rewritten as an integral with respect to the CM energy as Hall  _et al._
(2010); Elahi  _et al._ (2015):

$\displaystyle\dot{n_{1}}+3\mathcal{H}n_{1}$
$\displaystyle\approx\frac{T}{512\pi^{6}}\int_{s=4m_{\chi}^{2}}^{\infty}ds\leavevmode\nobreak\
d\Omega P_{ab}P_{12}\left|\mathcal{M}\right|_{ab\to 12}^{2}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{1}{\sqrt{s}}K_{1}\left(\frac{\sqrt{s}}{T}\right),$ (10)

where
$P_{ij}=\frac{1}{2\sqrt{s}}\sqrt{s-(m_{i}+m_{j})^{2}}\sqrt{s-(m_{i}-m_{j})^{2}}\to\frac{\sqrt{s}}{2}$
in the limit $m_{i,j}\to 0$. The BEQ in terms of the yield $Y_{\chi}=n_{1}/s$
can be written in the differential form as:

$\displaystyle-s(T).\mathcal{H}(T).T.\frac{dY_{\chi}^{\text{ann}}}{dT}$
$\displaystyle=\frac{T}{512\pi^{6}}\int_{s=4m_{\chi}^{2}}^{\infty}dsd\Omega
P_{ab}P_{12}\left|\mathcal{M}\right|_{ab\to 12}^{2}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\frac{1}{\sqrt{s}}K_{1}\left(\frac{\sqrt{s}}{T}\right).$ (11)

The total yield due then turns out to be:

$\displaystyle Y_{\chi}\left(T\right)$
$\displaystyle=\frac{1}{512\pi^{6}}\int_{T}^{T_{\text{max}}}\frac{dT^{{}^{\prime}}}{s(T^{{}^{\prime}}).\mathcal{H}(T^{{}^{\prime}})}\int_{s=4m_{\chi}^{2}}^{\infty}dsd\Omega\left(\frac{\sqrt{s}}{2}\right)^{2}$
$\displaystyle\left|\mathcal{M}\right|_{ab\to
12}^{2}\frac{1}{\sqrt{s}}K_{1}\left(\frac{\sqrt{s}}{T^{{}^{\prime}}}\right)$
(12)

where we have put $m_{i,j}=0$ just for convenience. The upper limit of the $T$
integration is considered to be $T_{\text{max}}=T_{\text{RH}}$ corresponding
to the reheat temperature of the universe. Eq. (12) should be divided into two
eras: one before EWSB with $T>160\leavevmode\nobreak\ \text{GeV}$ and the
other post-EWSB for $T\lesssim 160\leavevmode\nobreak\ \rm GeV$. For
$T_{\text{RH}}\gg T_{\text{EW}}$ the DM production is dominant before EWSB.
The lower limit on the $s$ integral can be approximated to zero if the DM mass
is negligible compared to reheat temperature of the universe. For a lower
reheat temperature, however, one can not ignore the masses of the particle
spectrum in the theory and for $T_{\text{RH}}\sim m_{i}$ the IR nature of the
yield becomes prominent Barman  _et al._ (2020a). Since in a fast expanding
universe the Hubble parameter gets modified according to Eq. (6), the
resulting DM yield also gets affected following Eq. (12). In the next section
we shall quantify this through analytical and numerical computation. Finally,
the relic abundance of the DM at present temperature can be obtained via:

$\displaystyle\Omega_{\text{DM}}h^{2}=\left(2.75\times
10^{8}\right)\frac{m_{\text{DM}}}{\text{GeV}}Y_{\text{DM}}^{\text{total}}(T_{0}),$
(13)

where we find $Y_{\text{DM}}^{\text{total}}(T_{0})$ by solving Eq. (12)
numerically. The important point here to note is that due to modified Hubble
rate $\mathcal{H}$ the DM yield and consequently the DM abundance is now a
function of the fast expansion parameters $n,T_{R}$ on top of DM mass, reheat
temperature and the effective scale.

###  III.1 A toy model

We consider a minimal dark sector comprising of a Dirac fermion singlet $\chi$
and a scalar singlet $\phi$. There exists a $Z_{4}$ symmetry under which only
the dark sector particles are charged as tabulated in Tab. 1. We assume the
fermion $\chi$ to be the lightest particle charged under the new symmetry and
hence a potential DM candidate. With this symmetry assigned one can then write
the effective Lagrangian for the dark sector as:

Symmetry | $\chi$ | $\phi$
---|---|---
$Z_{4}$ | i | -1
Table 1: Transformation of the new particles under $Z_{4}$ symmetry.

$\displaystyle\begin{split}&\mathcal{L}\supset\overline{\chi}\left(i\not{\partial}-m_{\chi}\right)\chi-
y_{\chi}\overline{\chi^{c}}\chi\phi+\frac{1}{\Lambda}\overline{\chi}\chi\left|H\right|^{2}\\\
&-V\left(H,\phi\right)\end{split}$ (14)

with the renormalizable scalar potential given by

$\displaystyle\begin{split}&V\left(H,\phi\right)\supset-\mu_{H}^{2}\left|H\right|^{2}+\mu_{\phi}^{2}\hskip
2.84526pt\phi^{2}+\lambda_{H}\left|H\right|^{4}+\lambda_{\phi}\hskip
2.84526pt\phi^{4}\\\ &\hskip
42.67912pt+\lambda_{H\phi}\left|H\right|^{2}\phi^{2}.\end{split}$ (15)

The $Z_{4}$ symmetry prevents us from writing any odd terms in the scalar
potential. Note that the portal coupling $\lambda_{H\phi}\sim\mathcal{O}(1)$
keeps the singlet scalar $\phi$ in equilibrium with the thermal bath. As the
new singlet scalar $\phi$ does not acquire any non-zero VEV, $\chi$ has only a
pure Dirac mass term $m_{\chi}$. Because of the $Z_{4}$ symmetry it is also
not possible to write Majorana mass terms for $\chi$ that otherwise would have
resulted in a pseudo-Dirac splitting De Simone  _et al._ (2010). The second
term in Eq. (14), due to its renormalizable nature, contributes to the DM
abundance via Infra-Red (IR) freeze-in due to $\phi\to\chi\chi$ decay which is
kinematically accessible as long as $m_{\chi}<m_{\phi}/2$. The third term
shall contribute to UV freeze-in Elahi  _et al._ (2015) for DM production via
$HH\to\chi\chi$ contact interaction before EWSB while after EWSB one can have
both UV and IR freeze-in contribution due to non-zero VEV of the Higgs that
results in renormalizable DM-SM interaction. For UV freeze-in the abundance is
sensitive to the maximum temperature of the thermal bath, which we assume to
be the reheat temperature $T_{\text{RH}}$222In principle, the maximum
temperature during reheating can be larger than $T_{\text{RH}}$ Chung  _et
al._ (1999).. This is in sharp contrast to the IR freeze-in scenario where the
two sectors communicate via renormalizable operators, and the DM abundance is
set by the IR physics i.e., the yield becomes maximum at low temperature,
typically at $T\sim m_{\text{DM}}$ Hall  _et al._ (2010); Bernal  _et al._
(2017). Now, the reheat temperature $T_{\text{RH}}$ of the universe is
somewhat loosely bounded. Typically, the lower bound on $T_{\text{RH}}$ comes
from the measurement of light element abundance during BBN, which requires
$T_{\text{RH}}\gtrsim 4.7\leavevmode\nobreak\ \rm MeV$. The upper bound, on
the other hand, may emerge from (i) cosmological gravitino problem Moroi  _et
al._ (1993); Kawasaki and Moroi (1995) in the context of supersymmetric
framework, that demands $T_{\text{RH}}\lesssim 10^{10}\leavevmode\nobreak\ \rm
GeV$ to prohibit thermal gravitino over production and (ii) simple
inflationary scenarios that require at most $T_{\text{RH}}\sim
10^{16}\leavevmode\nobreak\ \rm GeV$ Kofman  _et al._ (1997); Linde (1990) for
a successful inflation. In view of this, the maximum reheat temperature of the
universe can be regarded as a free parameter.

###  III.2 Relic abundance of the dark matter

A model-independent study of IR freeze-in in context of fast expansion has
been done in D’Eramo  _et al._ (2018). Here we would like to see the effect of
such fast expansion in the context of pure UV freeze-in before electroweak
symmetry breaking, together with IR freeze-in, once the electroweak symmetry
is broken. The requirement of obtaining right DM relic abundance via UV
freeze-in shall also put some bound on the reheat temperature, along with the
cut-off scale $\Lambda$ of the theory in the fast expanding scenario. Up to
this end we note that the contribution to DM abundance via IR freeze-in
through $\phi\to\chi\chi$ is negligible once we choose the Yukawa coupling
$y_{\chi}\lesssim 10^{-12}$ 333For example, for DM mass of 20 GeV and
$m_{\phi}=70\leavevmode\nobreak\ \rm GeV$ choosing $y_{\chi}=10^{-12}$ gives
rise to DM relic abundance of $\sim 10^{-3}$ dominantly via $\phi\to\chi\chi$
decay with a lifetime of $\sim 0.4\leavevmode\nobreak\ \rm sec$.. Such a
choice of coupling might look like somewhat “fine-tuned” but the smallness of
coupling for IR freeze-in is in itself a requirement to keep the DM out of
thermal equilibrium. The small Yukawa coupling thus prevents the dark sector
to reach thermal equilibrium within itself without hampering the rest of our
analysis. This makes UV freeze-in as the dominant process that can produce the
total observed DM relic via decay and /or annihilation which we will discuss
below.

In order to calculate the DM relic abundance by solving Eq. (12) we first
point out the possible annihilation and decay channels leading to DM pair
production above and below the electroweak symmetry breaking temperature
$T\simeq 160\leavevmode\nobreak\ \rm GeV$. Before EWSB i.e.,
$T>160\leavevmode\nobreak\ \rm GeV$, for the $SU(2)_{L}$ scalar doublet, the
propagating DOFs are the Goldstone bosons (GB):

$\displaystyle H=\begin{pmatrix}\varphi^{+}\\\ \varphi^{0}\end{pmatrix},$ (16)

which are then converted to the longitudinal DOF for the massive gauge bosons
after EWSB in unitary gauge:

$H=\begin{pmatrix}0\\\ \frac{h+v_{h}}{\sqrt{2}}\end{pmatrix},$

where $h$ is the CP-even SM-like Higgs field. After electroweak symmetry
breaking (EWSB) i.e., $T\lesssim 160\leavevmode\nobreak\ \rm GeV$ there can be
several $2\to 2$ annihilation channels (mediated by SM Higgs) that can produce
the DM via freeze-in, while the $1\to 2$ decay channels will only contribute
if $m_{h}>2m_{\chi}$ (see Fig. 1). Another SM Higgs mediated $2\to 2$ process
is shown in bottom left of Fig. 1 with the non-standard scalars $\phi$ in the
initial state. The amplitude for this process goes roughly as
$\lambda_{H\phi}^{2}v_{h}^{4}/\Lambda^{2}$. For
$\lambda_{H\phi}\lesssim\mathcal{O}\left(1\right)$ and
$\Lambda\gtrsim\mathcal{O}\left(10^{8}\right)\leavevmode\nobreak\ \rm GeV$,
however, this process has negligible contribution to the DM yield via IR
freeze-in. Using Eq. (12) it is possible to derive an approximate analytical
expression for DM yield via UV freeze-in (before EWSB). As all the SM
particles are also massless before EWSB, the squared amplitude for the 4-point
processes simply become

$\displaystyle\left|\mathcal{M}\right|^{2}\simeq\frac{s}{\Lambda^{2}}.$ (17)

Figure 1: DM yield via freeze-in before EWSB (top left)and after EWSB (top
right and bottom). Here SM stands for all the SM gauge bosons, quarks and
leptons. Here $H$ stands for the SM-like Higgs field before EWSB and $h$ is
the CP-even Higgs field after EWSB.

Figure 2: (a): Variation of DM yield $Y_{\chi}$ with $x=m_{\chi}/T$ for
different values of $n=\\{1,2,3,4\\}$ shown in green, blue, black and orange
respectively for $m_{\chi}=50\leavevmode\nobreak\ \rm GeV$,
$T_{\text{RH}}=10^{8}\leavevmode\nobreak\ \rm GeV$ and
$T_{R}=40\leavevmode\nobreak\ \rm MeV$. Each curve satisfies observed relic
abundance for different values of the cut-off scale
$\Lambda\approx\\{10^{15},10^{13},10^{11},10^{10}\\}\leavevmode\nobreak\ \rm
GeV$ respectively. (b): Same as (a) for $m_{\chi}=500\leavevmode\nobreak\ \rm
GeV$. The dashed black straight line in both plots corresponds to
$x_{\text{EW}}=m_{\chi}/T_{\text{EW}}$.

Figure 3: (a): Contours satisfying PLANCK observed relic abundance in the bi-
dimensional plane of $\Lambda$-$T_{\text{RH}}$ for different choices of
$n=\\{1,2,3\\}$ shown in green, blue and orange respectively where the DM mass
is fixed at 50 GeV and $T_{R}=40\leavevmode\nobreak\ \rm MeV$. (b): Same as
(a) for $m_{\chi}=500\leavevmode\nobreak\ \rm GeV$.

Considering the DM and $\phi$ mass to be negligible compared to
$\\{\Lambda,T_{\text{RH}}\\}$, for $T\gg T_{R}$, we can analytically obtain
the expression for DM yield before EWSB

$\displaystyle\begin{split}&Y_{\chi}=\left\\{\begin{array}[]{ll}\xi\left(\frac{T_{R}}{T_{\text{EW}}}\right)^{n/2}\frac{T_{\text{EW}}}{\Lambda^{2}\left(n-2\right)}\Biggl{[}1-\frac{T_{\text{RH}}}{T_{\text{EW}}}\Biggl{(}\frac{T_{\text{EW}}}{T_{\text{RH}}}\Biggr{)}^{n/2}\Biggr{]},\leavevmode\nobreak\
\leavevmode\nobreak\ (n\neq 2)\\\ \frac{\xi
T_{R}}{2\Lambda^{2}}\log\left[\frac{T_{\text{RH}}}{T_{\text{EW}}}\right]\leavevmode\nobreak\
\leavevmode\nobreak\ (n=2)\\\ \end{array}\right.\end{split}$ (18)

where $\xi=\frac{45M_{\text{pl}}}{1696\sqrt{106}\pi^{8}}$ is a constant of
mass dimension unity. After EWSB there will be IR dominated processes as well
which are proportional to $\left(v/\Lambda\right)^{2}$. Now, since all the SM
particles become massive after EWSB (together with the DM), hence one has to
find the DM yield numerically in this regime. The only possible decay channel
contributing to DM abundance is $h\to\chi\chi$ which is available if
$m_{\chi}<m_{h}/2$. We can obtain an approximate analytical expression for the
yield due to decay D’Eramo  _et al._ (2018) as

$\displaystyle\begin{split}&Y_{\chi}^{\text{decay}}=\int_{T_{\text{min}}}^{T_{\text{max}}}dT\frac{m_{h}^{2}\Gamma_{h\to\chi\chi}}{2\pi^{2}}\frac{K_{1}\left(m_{h}/T\right)}{s\left(T\right)\mathcal{H}\left(T\right)}\\\
&\approx\frac{\Gamma_{h\to\chi\chi}T_{R}^{n/2}}{2\pi^{2}m_{h}^{\frac{n}{2}+2}}\frac{M_{\text{pl}}}{A}\int_{0}^{\infty}dxx^{3+n/2}K_{1}\left(x\right)\\\
&=\frac{2M_{\text{pl}}\Gamma_{h\to\chi\chi}}{A\pi^{2}m_{h}^{2}}\Biggl{(}\frac{T_{R}}{m_{h}}\Biggr{)}^{n/2}\Biggl{[}2^{n/2}\Gamma\Bigl{(}\frac{6+n}{4}\Bigr{)}\Gamma\Bigl{(}\frac{10+n}{4}\Bigr{)}\Biggr{]},\end{split}$
(19)

where $A$ contains all the dimensionless constants and we assumed
$T_{\text{max}}\to\infty$ (which, in reality, should be $T_{\text{EW}}$) such
that $x_{\text{max}}\to 0$ with $T_{\text{min}}=0$. The decay width goes
roughly as $\Gamma_{h\to\chi\chi}\propto\left(v/\Lambda\right)^{2}$. In the
case of $2\to 2$ annihilation processes, in the limit where all the SM
particles are massless, the analytical form of the yield after EWSB turns out
to be

$\displaystyle\begin{split}&Y_{\chi}^{2\to
2}\sim\frac{m_{\chi}M_{\text{pl}}}{B\Lambda^{2}}\Biggl{(}\frac{T_{R}}{m_{\chi}}\Biggr{)}^{n/2}\int_{x_{\text{min}}}^{x_{\text{max}}}dxx^{n/2-2},\end{split}$
(20)

with $x_{\text{max}}=m_{\chi}/T_{0}$ and
$x_{\text{min}}=m_{\chi}/T_{\text{EW}}$, where $B$ contains all the
dimensionless constants. Taking all the masses of the particles into account,
however, such a simple analytical expression would be difficult to obtain and
one has to resort to numerical computation. The asymptotic nature of the yield
after EWSB can be determined from Eq. (19) and Eq. (20) which tells that the
yield is suppressed for larger $n$ as $x\to\infty$. Under the framework of
faster expansion one can not vary $T_{\text{RH}}$, $T_{R}$ and $n$ arbitrarily
as pointed out in D’Eramo  _et al._ (2018). The dominance of $\rho_{\varphi}$
at early universe sets an upper bound on the reheat temperature since beyond
$\rho_{\varphi}\sim M_{\text{pl}}^{4}$ this framework fails. This, together
with Eq. (6), sets an upper bound on the reheat temperature D’Eramo  _et al._
(2018)

$\displaystyle\begin{split}&H\sim\frac{T^{2}}{M_{\text{pl}}}\Bigl{(}T_{\text{RH}}/T_{R}\Bigr{)}^{n/2}\lesssim
M_{\text{pl}}\\\ &\implies T_{\text{RH}}\lesssim
M_{\text{pl}}\left(M_{\text{pl}}^{-1}T_{R}\right)^{n/n+4}.\end{split}$ (21)

Note that the above bound gets stronger for larger $n$. Additionally, for the
effective theory to remain valid $\Lambda>m_{\chi},T_{\text{RH}}$ is
mandatory. For processes appearing before and after EWSB we compute the DM
yield by performing numerical integrations using Mathematica. In Fig. 2 we
have illustrated how the DM yield varies with $x=m_{\chi}/T$ for a fixed
$T_{R}=40\leavevmode\nobreak\ \rm MeV$ and reheat temperature
$T_{\text{RH}}=10^{8}\leavevmode\nobreak\ \rm GeV$ for different choices of
$n$-values for different DM masses $m_{\chi}=\\{50,500\\}\leavevmode\nobreak\
\rm GeV$. For all these plots we have kept
$T_{\text{RH}}=10^{8}\leavevmode\nobreak\ \rm GeV$ which implies $n\lesssim 4$
following Eq. (21). First of all we see all these curves converge as
$x\to\infty$ since all of them satisfy the DM relic abundance for a fixed DM
mass of 50 GeV (on left panel) and 500 GeV (on right panel). Also note that
the DM yield builds up very early at a smaller $x\leavevmode\nobreak\ (T\sim
T_{\text{RH}})$ which is the typical feature of UV freeze-in when
$T_{\text{RH}}\gg m$ Elahi  _et al._ (2015); Barman  _et al._ (2020a). For
larger $n$, as one can notice from Eq. (12), the DM yield gets suppressed
since $\mathcal{H}$ becomes larger according to Eq. (6) for a fixed
$T_{\text{RH}}$. This shows that the freeze-in production of DM in a fast
expanding universe is dramatically suppressed which has also been reported
earlier D’Eramo  _et al._ (2018). Thus, a larger $n$ requires a smaller
$\Lambda$ to achieve right DM abundance for a fixed $T_{\text{RH}}$ as
$Y_{\chi}\propto\frac{1}{\Lambda^{2}}$. The orange curves $(n=4)$ in Fig. 2
correspond to $\Lambda\sim 10^{9}\leavevmode\nobreak\ \rm GeV$ which is about
five orders less than that for the red curves (standard case) corresponding to
which $\Lambda\sim 10^{15}\leavevmode\nobreak\ \rm GeV$. This exactly what is
reflected in Fig. 3 where we show DM abundance satisfying contours in
$\Lambda$-$T_{\text{RH}}$ plane for DM mass of
$m_{\chi}=\\{50,500\\}\leavevmode\nobreak\ \rm GeV$ in the left and in the
right panel respectively considering all processes appearing before and after
EWSB. For a 50 GeV DM $h\to\chi\chi$ decay contributes significantly after
EWSB via IR freeze-in Barman  _et al._ (2020b). As a result we see the
horizontal part in the red and green curves for $T_{\text{RH}}\lesssim
10\leavevmode\nobreak\ \rm TeV$. This arises due to the IR nature of the
freeze-in where the DM abundance depends no more on the reheat temperature,
rather on the mass of the decaying mother particle (which is Higgs in this
case). For larger $n$ although the curves look absolutely flat (e.g., the blue
one) but as one can see from the inset figure their variation over $\Lambda$
is almost negligible compared to the standard case or $n=1$ case. For a
heavier DM mass, typically $m_{\chi}\gg T_{\text{EW}}$, most of the
contribution to the DM yield comes from processes before EWSB. As a
consequence we observe typical UV nature in the contours on the right panel of
Fig. 3 where the cut-off scale rises linearly with the reheat temperature. For
small enough $T_{\text{RH}}\ll 1\leavevmode\nobreak\ \rm TeV$ the IR nature of
the freeze-in might show up even for heavier DM but we choose to confine
ourselves to $T_{\text{RH}}\gtrsim 1\leavevmode\nobreak\ \rm TeV$.

We conclude this section by inferring that one has to increase the effective
dimensionful coupling $\propto\Lambda^{-1}$ with respect to the standard
cosmological background in order to get the right DM abundance in a fast
expanding background. This implies that for non-renormalizable interactions
with dimension larger than five, it is possible to bring down
$\Lambda\leavevmode\nobreak\ \sim\text{TeV}$ for $n>0$ that provides larger
detection prospect for UV freeze-in. For dimension five interactions, since
for larger $n$ the DM yield gets more and more suppressed, hence we stick only
up to $n=4$. For such a choice of $n$, as we have illustrated, it is possible
to obtain observed DM relic abundance for different choices of $T_{R}$ and
$T_{\text{RH}}$ obeying Eq. (7) and Eq. (21) by suitably modifying the cut-off
scale $\Lambda$. As we shall see in the next section, for such a choice of
$n,T_{R}$ it is possible to generate a strong first order phase transition
that is capable of explaining the NANOGrav 12.5 yrs data.

##  IV Stochastic gravitational wave signal in a fast expanding universe

In general, an extended scalar sector may give rise to a strong first order
phase transition (SFOPT)resulting in stochastic gravitational wave (GW) with
detectable frequencies (see, for example Kakizaki  _et al._ (2015)). In our
case the presence of an extra scalar singlet $\phi$ does that job. In the
context of the present model we would like to explore how the frequency and
amplitude of such a stochastic GW signal can get modified due to the fast
expansion of the universe. Our particular interest is to explore if a simple
dark sector can explain the NANOGrav data when the cosmological background is
modified. As we shall see, the stochastic GW signal that can explain the
reported NANOGrav data Arzoumanian  _et al._ (2020) can only be generated in a
fast expanding universe for some specific choices of $T_{R},n$.

####  IV.0.1 Finite Temperature Effective Potential

To construct the finite temperature effective potential, we add two
temperature corrected terms with the tree-level potential. The effective
potential at a finite temperature $T$ then can be written as Wainwright (2012)

$\displaystyle V_{\text{eff}}=V_{\text{tree-
level}}+V_{1-loop}^{T=0}+V_{1-loop}^{T\neq 0}\,\,,$ (22)

where $V_{\text{tree-level}}$, $V_{1-loop}^{T=0}$ and $V_{1-loop}^{T\neq 0}$
are the zero temperature tree level potential, zero temperature Coleman-
Weinberg (CW) one-loop potential and the finite temperature one-loop potential
respectively. As mentioned earlier, the scalar sector under consideration
contains a real scalar singlet on top of the SM-like Higgs doublet which leads
to Eq. (15). The SM Higgs $H$ can be represented as

$\displaystyle H=\left(\begin{array}[]{c}G^{+}\\\
\dfrac{1}{\sqrt{2}}\left(v+h+iG^{0}\right)\\\ \end{array}\right)\,\,,$ (25)

where $G^{+}$ and $G_{0}$ are the charged and the neutral Goldstone bosons
after the spontaneous symmetry breaking. To explore the electroweak phase
transition we replace the scalar fields $H$ and $\phi$ with their VEVs $v$,
$v_{\phi}$ respectively. Now Eq. (15) can be expressed as

$\displaystyle\begin{split}&V_{\rm tree-
level}^{T=0}=-\dfrac{1}{2}\mu_{H}^{2}v^{2}+\dfrac{1}{2}\mu_{\phi}^{2}v_{\phi}^{2}+\dfrac{1}{4}\lambda_{H}v^{4}\\\
&+\dfrac{1}{4}\lambda_{\phi}v_{\phi}^{4}+\dfrac{1}{4}\lambda_{H\phi}v^{2}v_{\phi}^{2}.\end{split}$
(26)

Note that, although the scalar $\phi$ do not acquire any vacuum expectation
value (VEV) at $T=0$ (to keep the $Z_{2}$ symmetry intact), but its VEV can be
generated at a finite temperature. Here the classical VEVs $v,v_{\phi}$ change
with temperature i.e., behave like a field and at the zero temperature it
tends to the classical fixed values. The zero-temperature CW one-loop
effective potential is given by Wainwright (2012); Coleman and Weinberg (1973)

$\displaystyle
V_{1-loop}^{T=0}=\pm\dfrac{1}{64\pi^{2}}\sum_{i}n_{i}m_{i}^{4}\left[\log\dfrac{m_{i}^{2}}{Q^{2}}-C_{i}\right]\,\,,$
(27)

where the ‘$+$’ and ‘$-$’ symbols appear due to bosons and fermions
respectively. The quantity $Q$ denotes the renormalization scale which we take
$Q=v=246.22$ GeV in our work. The summation $i$ is over all the particle
species associated in the scalar potential with
$i\in(h,\phi,G^{\pm},G^{0},W,Z,t$). In Eq. (27), the quantities $n_{i}$,
$m_{i}$ and $C_{i}$ represents the number of degrees of freedom, the field-
dependent masses and renormalization-scheme-dependent numerical constant of
the $i$th particle species respectively. The associated degrees of freedom of
the particles are $(n_{W^{\pm}})_{L}=4$, $(n_{W^{\pm}})_{T}=2$,
$(n_{Z})_{L}=2$, $(n_{Z})_{T}=1$, $n_{t}=12$, $n_{G^{\pm}}=2$ and
$n_{h,\phi,G^{0}}=1$. The values of the renormalization-scheme-dependent
constant $C_{i}$ are $(C_{W,Z})_{T}=1/2$ for the transverse component of
$W,\hskip 2.84526ptZ$ boson, $(C_{W,Z})_{L}=3/2$ for the longitudinal
component of $W,\hskip 2.84526ptZ$ boson and for the other particles,
$C_{h,\phi,G^{+},G^{-},G^{0},t}=3/2$. We perform our calculations by
considering the Landau gauge where the Goldstone bosons are massless at $T=0$
and also the theory is free from ghost contributions Basler  _et al._ (2017).
With this, the finite temperature one-loop effective potential can be
expressed as Wainwright (2012)

$\displaystyle V_{1-loop}^{T\neq
0}=\dfrac{T^{4}}{2\pi^{2}}\sum_{i}n_{i}J_{\pm}\left[\dfrac{m_{i}^{2}}{T^{2}}\right]\,\,,$
(28)

where

$\displaystyle
J_{\pm}\left(\dfrac{m_{i}^{2}}{T^{2}}\right)=\pm\int_{0}^{\infty}dy\hskip
2.84526pty^{2}\log\Biggl{(}1\mp
e^{-\sqrt{y^{2}+\dfrac{m_{i}^{2}}{T^{2}}}}\Biggr{)}\,\,.$ (29)

We include the thermal correction to the boson masses by applying daisy
resummation method Arnold and Espinosa (1993) as $\mu^{\prime
2}_{H}(T)=\mu_{H}^{2}+c_{1}T^{2}$ and $\mu^{\prime
2}_{\phi}(T)=\mu_{\phi}^{2}+c_{2}T^{2}$ where

$\displaystyle
c_{1}=\dfrac{6\lambda_{H}+2\lambda_{H\phi}}{12}+\dfrac{3g^{2}+{g^{\prime}}^{2}}{16}+\dfrac{y_{t}^{2}}{4}\,\,$
(30)

and

$\displaystyle c_{2}=\dfrac{6\lambda_{\phi}+2\lambda_{H\phi}}{12}\,\,,$ (31)

with the couplings $g$, $g^{\prime}$ and $y_{t}$ as the $SU(2)_{L}$ gauge
coupling, $U(1)_{Y}$ gauge coupling and top Yukawa coupling of the SM
respectively. We also include the thermal corrections to the $W$ and $Z$ boson
masses with the modified thermal masses Basler  _et al._ (2017)

$\displaystyle
m_{W}^{2}\left(T\right)=\dfrac{1}{4}g^{2}v^{2}+2g^{2}T^{2}\,\,,$ (32)

and

$\displaystyle\begin{split}&m_{Z}^{2}\left(T\right)=\dfrac{1}{8}\left(g^{2}+{g^{\prime}}^{2}\right)v^{2}+\left(g^{2}+{g^{\prime}}^{2}\right)T^{2}\\\
&+\dfrac{1}{8}\sqrt{\left[\left(g^{2}-{g^{\prime}}^{2}\right)^{2}\left(64T^{2}+16T^{2}v^{2}\right)+\left(g^{2}+{g^{\prime}}^{2}\right)^{2}v^{4}\right]}\,\,.\end{split}$
(33)

With this set-up we then proceed to analyze the GW production due to strong
first order phase transition (SFOPT).

####  IV.0.2 Gravitational wave production from SFOPT

Here we discuss the possible production mechanism of GW from the strong first-
order electroweak phase transition originating mainly via bubble collisions
Kosowsky  _et al._ (1992a); Paul  _et al._ (2019); Kosowsky and Turner (1993);
Huber and Konstandin (2008); Kosowsky  _et al._ (1992b); Kamionkowski  _et
al._ (1994); Caprini  _et al._ (2008), sound waves induced by the bubbles
running through the cosmic plasma Hindmarsh  _et al._ (2014); Giblin and
Mertens (2013, 2014); Hindmarsh  _et al._ (2015) and turbulence induced by the
bubble expansions in the cosmic plasma Caprini and Durrer (2006); Kahniashvili
_et al._ (2008a, b, 2010); Caprini  _et al._ (2009). The total GW intensity
$\Omega_{\text{GW}}{\rm h}^{2}$ for a particular frequency $f$ can be
estimated by considering the contributions of the above three processes and
can be expressed as Kosowsky  _et al._ (1992a); Paul  _et al._ (2019);
Kosowsky and Turner (1993); Huber and Konstandin (2008); Kosowsky  _et al._
(1992b); Kamionkowski  _et al._ (1994); Caprini  _et al._ (2008); Hindmarsh
_et al._ (2014); Giblin and Mertens (2013, 2014); Hindmarsh  _et al._ (2015);
Caprini and Durrer (2006); Kahniashvili  _et al._ (2008a, b, 2010); Caprini
_et al._ (2009); Barman  _et al._ (2020c); Paul  _et al._ (2020)

$\displaystyle\Omega_{\text{GW}}{\rm h}^{2}=\Omega_{\text{col}}{\rm
h}^{2}+\Omega_{\text{SW}}{\rm h}^{2}+\Omega_{\text{turb}}{\rm h}^{2}\,\,.$
(34)

The expression for the GW intensity from the component of bubbles collision
$\Omega_{\text{col}}{\rm h}^{2}$ is given by

$\displaystyle\begin{split}&\Omega_{\text{col}}{\rm h}^{2}=1.67\times
10^{-5}\left(\dfrac{\beta}{\mathcal{H}}\right)^{-2}\dfrac{0.11v_{w}^{3}}{0.42+v_{w}^{2}}\\\
&\left(\dfrac{\kappa\alpha}{1+\alpha}\right)^{2}\left(\dfrac{g_{*}}{100}\right)^{-\frac{1}{3}}\dfrac{3.8\left(\dfrac{f}{f_{col}}\right)^{2.8}}{1+2.8\left(\dfrac{f}{f_{\text{col}}}\right)^{3.8}}\,\,\,,\end{split}$
(35)

where the parameter $\beta$ has the form

$\displaystyle\beta=\left[\mathcal{H}T\dfrac{d}{dT}\left(\dfrac{S_{3}}{T}\right)\right]\bigg{|}_{T_{n}}\,\,,$
(36)

where $S_{3}\left(T\right)$ is the Euclidean action of the critical bubble. In
general, the nucleation of the bubbles occur at a temperature $T_{n}$ where
the condition $S_{3}\left(T_{n}\right)/T_{n}\approx 140$ Wainwright (2012) is
satisfied. In Eq. (36) $\mathcal{H}$ denotes the Hubble parameter at the
nucleation temperature $T_{n}$ which is given by Eq. (6) for the fast
expanding universe. The bubble wall velocity $v_{w}$ is computed using the
most general expression for the wall velocity Steinhardt (1982); Kamionkowski
_et al._ (1994); Dev  _et al._ (2019); Paul  _et al._ (2019)

$\displaystyle
v_{w}=\dfrac{1/\sqrt{3}+\sqrt{\alpha^{2}+2\alpha/3}}{1+\alpha}\,\,.$ (37)

The quantity $\kappa$ in Eq. (LABEL:eq:66) is given by

$\displaystyle\kappa=1-\dfrac{\alpha_{\infty}}{\alpha}\,\,,$ (38)

with Shajiee and Tofighi (2019); Caprini  _et al._ (2016)

$\displaystyle\begin{split}&\alpha_{\infty}=\dfrac{30}{24\pi^{2}g_{*}}\left(\dfrac{v_{n}}{T_{n}}\right)^{2}\Biggl{[}6\left(\dfrac{m_{W}}{v}\right)^{2}+3\left(\dfrac{m_{Z}}{v}\right)^{2}\\\
&+6\left(\dfrac{m_{t}}{v}\right)^{2}\Biggr{]}\,\,.\end{split}$ (39)

where $v_{n}$ is the Higgs VEV at the nucleation temperature $T_{n}$, $m_{W}$,
$m_{Z}$ and $m_{t}$ are the masses of the gauge bosons $W$, $Z$ and the top
quark $t$ respectively. The phase transition parameter $\alpha$ can be
expressed as

$\displaystyle\alpha=\left[\dfrac{\rho_{\text{vac}}}{\rho^{*}_{\text{rad}}}\right]\bigg{|}_{T_{n}}\,\,.$
(40)

where $\rho_{*}^{\text{rad}}$ is the background energy density of the plasma
and $\rho_{\text{vac}}$ is the energy density difference between false and
true vacuum during the electroweak phase transition. The quantities
$\rho_{\text{vac}}$ and $\rho_{*}^{\text{rad}}$ have the following form

$\displaystyle\rho_{\text{vac}}=\left[\left(V_{\text{eff}}^{\text{high}}-T\dfrac{dV_{\text{eff}}^{\text{high}}}{dT}\right)-\left(V_{\text{eff}}^{\text{low}}-T\dfrac{dV_{\text{eff}}^{\text{low}}}{dT}\right)\right]\,\,,$
(41)

and

$\displaystyle\rho^{*}_{\text{rad}}=\dfrac{g_{*}\pi^{2}T_{n}^{4}}{30}\,\,.$
(42)

In Eq. (LABEL:eq:66) the peak frequency $f_{\text{col}}$ from the bubble
collisions is given by

$\displaystyle\begin{split}&f_{\text{col}}=16.5\times 10^{-6}\hskip
2.84526pt\text{Hz}\left(\dfrac{0.62}{v_{w}^{2}-0.1v_{w}+1.8}\right)\left(\dfrac{\beta}{\mathcal{H}}\right)\\\
&\left(\dfrac{T_{n}}{100\hskip
2.84526pt\text{GeV}}\right)\left(\dfrac{g_{*}}{100}\right)^{\frac{1}{6}}\,\,.\end{split}$
(43)

Next, the expression for GW intensity from the component of sound wave (SW) is
given by

$\displaystyle\begin{split}&\Omega_{\text{SW}}{\rm h}^{2}=2.65\times
10^{-6}\left(\dfrac{\beta}{\mathcal{H}}\right)^{-1}v_{w}\left(\dfrac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\dfrac{g_{*}}{100}\right)^{-\frac{1}{3}}\\\
&\left(\dfrac{f}{f_{\text{SW}}}\right)^{3}\left[\dfrac{7}{4+3\left(\dfrac{f}{f_{\text{SW}}}\right)^{2}}\right]^{\frac{7}{2}}\,\,,\end{split}$
(44)

where the parameter $\kappa_{v}$

$\displaystyle\kappa_{v}=\dfrac{\alpha_{\infty}}{\alpha}\left[\dfrac{\alpha_{\infty}}{0.73+0.083\sqrt{\alpha_{\infty}}+\alpha_{\infty}}\right]\,\,.$
(45)

The peak frequency $f_{\text{SW}}$ from the sound wave contribution reads

$\displaystyle f_{\text{SW}}=1.9\times 10^{-5}\hskip
2.84526pt\text{Hz}\left(\dfrac{1}{v_{w}}\right)\left(\dfrac{\beta}{\mathcal{H}}\right)\left(\dfrac{T_{n}}{100\hskip
2.84526pt\text{GeV}}\right)\left(\dfrac{g_{*}}{100}\right)^{\frac{1}{6}}\,\,.$
(46)

For the contribution of SW to the total GW intensity depends on the Hubble
time scale. If it survives more than a Hubble time then the expression in Eq.
(LABEL:eq:75) will be valid, otherwise we need to include a factor called
suppression factor to the SW component of the GW intensity. Following Caprini
_et al._ (2016); Ellis  _et al._ (2018, 2019) we compute the suppression
factor $\dfrac{\mathcal{H}R_{*}}{\bar{U}_{f}}$ (where $\bar{U}_{f}$ is the
root-mean-square (RMS) fluid velocity and $R_{*}$ is the mean bubble
separation) to check whether the SW components lasts more than a Hubble time.

Figure 4: Plot of calculated GW intensities as a function of frequency for
different choices of of $\\{n,T_{R}\\}$. The sensitivity curves of NANOGrav11,
NANOGrav12.5 (purple coloured rectangular box), future GW detectors such as
ALIA, BBO, DECIGO, LISA, aLIGO and aLIGO+ are also shown in the same plane.

Lastly, the expression for GW intensity from the turbulence in the plasma
$\Omega_{\text{turb}}{\rm h}^{2}$ is given by

$\displaystyle\begin{split}&\Omega_{\text{turb}}{\rm h}^{2}=3.35\times
10^{-4}\left(\dfrac{\beta}{\mathcal{H}}\right)^{-1}v_{w}\left(\dfrac{\epsilon\kappa_{v}\alpha}{1+\alpha}\right)^{\frac{3}{2}}\left(\dfrac{g_{*}}{100}\right)^{-\frac{1}{3}}\\\
&\dfrac{\left(\dfrac{f}{f_{\text{turb}}}\right)^{3}\left(1+\dfrac{f}{f_{\text{turb}}}\right)^{-\frac{11}{3}}}{\left(1+\dfrac{8\pi
f}{h_{*}}\right)}\,\,,\end{split}$ (47)

with

$\displaystyle h_{*}=16.5\times 10^{-6}\hskip
2.84526pt\text{Hz}\left(\dfrac{T_{n}}{100\hskip
2.84526pt\text{GeV}}\right)\left(\dfrac{g_{*}}{100}\right)^{\frac{1}{6}}\,\,.$
(48)

Here $\epsilon=0.1$ and the peak frequency $f_{\text{turb}}$ from the
turbulence mechanism has the form

$\displaystyle f_{\text{turb}}=2.7\times 10^{-5}\hskip
2.84526pt\text{Hz}\left(\dfrac{1}{v_{w}}\right)\left(\dfrac{\beta}{\mathcal{H}}\right)\left(\dfrac{T_{n}}{100\hskip
2.84526pt\text{GeV}}\right)\left(\dfrac{g_{*}}{100}\right)^{\frac{1}{6}}\,\,,$
(49)

BP | $m_{h}$ | $m_{\phi}$ | $\lambda_{H}$ | $\lambda_{\phi}$ | $\lambda_{H\phi}$
---|---|---|---|---|---
| in GeV | in GeV |  |  |
1 | 125.09 | 70 | 0.129 | 0.1 | 0.1
Table 2: Chosen benchmark point to demonstrate strong first-order phase transition. BP | $v_{c}$ | $T_{c}$ | $\dfrac{v_{c}}{T_{c}}$ | $v_{n}$ | $T_{n}$ | $\alpha$
---|---|---|---|---|---|---
| in GeV | in GeV |  | in GeV | in GeV |
1 | 140.50 | 134 | 1.04 | 151.20 | 134.58 | 0.0028
Table 3: Values of thermal parameters for the chosen BP.

To investigate the observational signatures of such stochastic GW, our
estimated model-dependent GW intensities are compared with the future space-
based and ground-based detectors such as ALIA, BBO, DECIGO, aLIGO, aLIGO+ and
pulsar timing arrays such as NANOGrav11 and NANOGrav12.5. We are typically
interested in the GW spectra from FOPT in a fast expanding universe which lie
within the NANOGrav sensitivity. For calculating the thermal parameters
related to the phase transition we use the publicly available CosmoTransition
package Wainwright (2012) and compute the total GW intensities using Eqs.
(34)-(48). We choose a benchmark point (BP) with
$m_{\phi}=70\leavevmode\nobreak\ \rm GeV$ which is safe from Higgs invisible
decay constraint. Note that all the scalar couplings giving rise to strong
first order PT are $\lambda_{i}\sim 0.1$. As we have checked, a larger scalar
mass does not alter our results.

The thermal parameters $\\{v_{c}$, $T_{c}$, $v_{n}$, $T_{n}$, $\alpha\\}$ for
the chosen BP are tabulated in Tab. 3. The strength of the FOPT depends on the
order parameter $v_{c}/T_{c}$ and the transition is said to be strong enough
if it satisfies the condition $v_{c}/T_{c}\geqslant 1$, which in our case
turns out to be $v_{c}/T_{c}=1.04$. In Tab. 4 we present the results for GW
intensity and corresponding peak frequency for different values of $n$ and
$T_{R}$. From Tab. 4 and Fig. 4 we observe that the GW intensity increases and
peak shifts to the lower frequency region for smaller values of $T_{R}$ at a
fixed $n$. The GW intensities as a function of frequency are plotted in Fig. 4
for different choices of $n$ and $T_{R}$. We compare our model-dependent GW
intensities with the sensitivity plots of PTAs such as NANOGrav11,
NANOGrav12.5 and future space-based, ground-based gravitational wave detectors
such as ALIA, BBO, DECIGO, LISA, aLIGO and aLIGO+ 111In this work we use the
power-law-integrated sensitivity curves following Thrane and Romano (2013);
Dev  _et al._ (2019). One can also use the other methods to represent the
sensitivity curves Schmitz (2020); Moore  _et al._ (2015); Alanne  _et al._
(2020). From the Fig. 4 one observes that for $n=3$,
$T_{R}=40\leavevmode\nobreak\ \rm MeV$ and $n=4$,
$T_{R}=200\leavevmode\nobreak\ \rm MeV$ the GW intensity lies within the
sensitivity curves of NANOGrav-12.5. There are other choices of
$\\{n,T_{R}\\}$ that also give rise to detectable GW intensities which fall
within the reach of detectors like ALIA, BBO, DECIGO, LISA and NANOGrav11.

Figure 5: The blue and black solid lines correspond to the predictions of
$\gamma_{\text{CP}}$ and $A_{\text{CP}}$ (with $A_{\text{CP}}$ the
characteristic GW strain amplitude at $f=\text{yr}^{-1}$ and
$\gamma_{\text{CP}}$ the spectral index of the pulsar timing-residual cross-
power spectral density) obtained for different fixed values of $\\{n,T_{R}\\}$
which gives rise to $\Omega_{\text{GW}}h^{2}$ required to explain the NANOGrav
result. The dashed and solid red curves correspond to to the $1\sigma$ and
$2\sigma$ $\gamma_{\text{CP}}$-$A_{\text{CP}}$ contours obtained by NANOGrav
12.5 yrs dataset.

The results of PTA GW searches are generally reported in terms of a pulsar
timing-residual cross-power spectral density which has a frequency dependence
of the form $f^{-\gamma_{\text{CP}}}$. Corresponding power spectrum of the
characteristic GW strain is usually approximated as a power-law of the form

$\displaystyle
h_{c}\left(f\right)=A_{\text{CP}}\Biggl{(}\frac{f}{f_{\text{yr}}}\Biggr{)}^{\alpha_{\text{CP}}}$
(50)

with $\alpha_{\text{CP}}=\left(3-\gamma_{\text{CP}}\right)/2$ which is -2/3
for for a population of inspiraling supermassive black holes (SMBHB) in
circular orbits whose evolution is dominated by GW emission. Most importantly,
The reference frequency is given by $f=\text{yr}^{-1}$. The power spectrum is
related to the GW energy density via Vagnozzi (2020)

$\displaystyle\Omega_{\text{GW}}=\frac{2\pi^{2}}{3\mathcal{H}^{2}}f^{2}h_{c}^{2}\left(f\right).$
(51)

$n$ | $T_{R}$ | $\dfrac{\beta}{\mathcal{H}}$ | $f_{{\rm{peak}}}$ | $\Omega_{\text{GW}}{\rm h}^{2}_{{\rm{peak}}}$
---|---|---|---|---
| in MeV |  | in Hz |
- | - | 5264.70 |  3.45$\times 10^{-2}$ |  1.95$\times 10^{-19}$
1 | 16 | 57.53 |  3.76$\times 10^{-4}$ |  2.72$\times 10^{-14}$
2 | 10 | 0.39 |  2.57$\times 10^{-6}$ |  3.49$\times 10^{-11}$
3 | 40 | 0.027 |  1.78$\times 10^{-7}$ |  7.30$\times 10^{-9}$
4 | 200 | 0.012 |  7.68$\times 10^{-8}$ |  3.91$\times 10^{-8}$
4 | 100 | 0.0029 |  1.92$\times 10^{-8}$ |  6.26$\times 10^{-7}$
Table 4: The calculated values of peak frequencies and corresponding GW
intensity for the standard cosmology scenario (first row) and non-standard
scenarios with different choices of $T_{R}$.

The GW search results due to PTA are reported in terms of joint
$A_{\text{CP}}-\gamma_{\text{CP}}$ posterior distributions or as the posterior
distribution at a fiducial value of $\gamma_{\text{CP}}$ Arzoumanian  _et al._
(2020); Vagnozzi (2020). NANOGrav collaboration has obtained joint constraints
on $A_{\text{CP}}$ and $\gamma_{\text{CP}}$ by fitting the power-law in Eq.
(50) Arzoumanian  _et al._ (2020). Since Eq. (51) is a function of $n,T_{R}$
(through $\mathcal{H}$), hence it is possible to project limits on the fast
expansion parameters by exploiting Eq. (50). Now, as we have shown in Fig. 4,
corresponding to $n=3,4$ the GW intensities overlap with the NANOGrav
sensitivity curves. We thus cast $n-T_{R}$ constraints from Fig. 4 onto the
$A_{\text{CP}}-\gamma_{\text{CP}}$ plane with corresponding 1$\sigma$ and
2$\sigma$ contours from NANOGrav result Arzoumanian  _et al._ (2020). Here we
observe that for a fixed $n=4$ a larger $T_{R}$ shifts the straight line
contour (orange) to the edge of the 2$\sigma$ contour (red thick), while
smaller $n$ (black) is closer to the 1$\sigma$ contour (red dashed). From Fig.
4 and 5 we see that NANOGrav results not only support a faster expansion in
the pre-BBN era, also favour $n\lesssim 4$ for certain choices of $T_{R}$.
Based on the analysis in Sec. III and Sec. IV we can infer that while right
relic abundance for the DM can be obtained for any $n>0$ by properly tuning
the cut-off scale, reheat temperature and $T_{R}$ but explaining NANOGrav
results from a SFOPT together with the right DM abundance necessarily demands
$n\lesssim 4$ and
$T_{R}\lesssim\mathcal{O}\left(100\right)\leavevmode\nobreak\ \rm MeV$.

##  V Possible UV completion

Here we would like to elucidate a possible UV complete framework of the
effective Lagrangian we discussed in Eq. (14). There might be several
possibilities to obtain an effective interaction of the form in Eq. (14). For
an instance, the authors in Dutta Banik  _et al._ (2018) have discussed one
such alternative. However, in their case the $Z_{4}$ is broken by the non-zero
VEV of the singlet scalar $\phi$ leading to a remnant $Z_{2}$ symmetry under
which the extra fermions are odd (DM). On the other hand, in Konar  _et al._
(2020) a $Z_{2}$ symmetry has been imposed under which all the fields in the
dark sector are odd. Here we would like to furnish the simplest possible model
remembering the fact that in the present scenario the singlet scalar does not
acquire any tree-level VEV. This is in contrast to the set-up discussed in
Fairbairn and Hogan (2013) where the singlet fermionic WIMP-like DM is
stabilized by a $Z_{2}$ symmetry and the singlet scalar acquires a VEV
providing a strongly first order electroweak phase transition.

Fields | $SU(2)_{L}$ | $U(1)_{Y}$ | $Z_{4}$
---|---|---|---
$\chi$ | 1 | 0 | $i$
$\psi^{T}:\left(\psi^{0},\psi^{-}\right)$ | 1 | 1/2 | $-i$
$\phi$ | 1 | 0 | $-1$
Table 5: Charges of different fields appearing in a possible UV completion of
the Lagrangian in Eq. (14).

We assume a dark sector that is made up of a vector like lepton (VLL) singlet
$\chi$ and a $SU(2)_{L}$ VLL doublet $\psi$ with the charges as mentioned in
Tab. 5. As mentioned earlier, we consider the new particles transform under a
$Z_{4}$ symmetry with their charges as in Tab. 5. With the given
$\text{SM}\times Z_{4}$ symmetry the renormalizable Lagrangian can be written
as

$\displaystyle\begin{split}&\mathcal{L}\supset
m_{\chi}\overline{\chi}\chi+M_{\psi}\overline{\psi}\psi+Y\overline{\psi}\widetilde{H}\chi+y_{\chi}\overline{\chi^{c}}\chi\phi+\text{H.c.}\end{split}$
(52)

At a scale $\mu\ll M_{\psi}$ the heavy VLLs can be integrated out, resulting
in operators of the form
$\left(H^{\dagger}H\right)\left(\overline{\chi}\chi\right)/M_{\psi}$. If the
reheat temperature of the universe is smaller than $M_{\psi}$ then the
additional VLLs are not present in the thermal bath and the resulting DM-SM
interaction mimics a dimension five operator with cut-off scale $\Lambda\equiv
M_{\psi}$.

##  VI Conclusion

In this paper we have studied a simple framework where it is possible to
address both dark matter (DM) relic abundance via freeze-in and the recent
NANOGrav result from a strong first order phase transition (SFOPT) in a fast
expanding universe. In order to accomplish them together, we consider a dark
sector consisting of a SM gauge singlet fermion $\chi$ and a scalar singlet
$\phi$ both charged under a $Z_{4}$ symmetry. We consider $\chi$ to be the
lightest stable fermion in the dark sector that can serve as a viable DM
candidate. Being singlet under SM gauge symmetry the dark fermion does not
have any interactions in the renormalizable level with the SM fields. The
scalar, on the other hand, remains in equilibrium with the thermal bath via
non-negligible portal coupling and is responsible in generating SFOPT.

We focus on the DM genesis via freeze-in, typically UV freeze-in, as the dark
sector and the visible sector are conjoined through non-renormalizable
interaction of dimension five. The DM yield is computed by taking into account
all processes leading to DM production both before and after the electroweak
symmetry breaking (EWSB). While the processes before EWSB (dominantly UV) are
only $2\to 2$ scattering channels producing DM pair in the final state, after
EWSB a host of annihilation and decay channels appear which contribute to both
UV and IR freeze-in of the DM. From the requirement of obtaining the observed
DM abundance we put constraints on the fast expansion parameters $T_{R},n$ for
fixed DM masses. This also put constraints on the other two relevant free
parameters of the model, namely the effective scale $\Lambda$ and the reheat
temperature $T_{\text{RH}}$ which we consider to $T_{\text{RH}}\gtrsim
1\leavevmode\nobreak\ \rm TeV$. We see, in accordance with previously reported
results, in a faster-than-standard expanding universe, freeze-in production is
hugely suppressed for any $n>0$. This, in turn, implies that to produce enough
DM via freeze-in to match observations, one has to bring down the effective
scale of interaction (equivalently, the DM-SM coupling) substantially with
respect to the standard cosmological scenario.

The presence of the scalar singlet facilitates strong first order phase
transition (SFOPT) which is otherwise impossible to achieve within the SM
particle content. In the presence of fast expansion, we observe, such SFOPT
gives rise to gravitational wave (GW) signal with intensities that are
detectable by several space-based and ground-based GW detectors. In a fast
expanding universe such strong GW signals can be generated for
$\lambda_{H\phi}\gtrsim 0.1$ unlike the standard cosmological scenario which
requires large enough portal coupling. The GW signal thus generated can
explain the 12.5-year pulsar timing data due to NANOGrav. However, SFOPT
supporting NANOGrav result is possible only for specific choices of the
parameters that control the fast expansion. Typically, we see $n\lesssim 4$
lies within the 2$\sigma$ posterior contours for amplitude and spectral slope
obtained by the NANOGrav collaboration. For all such $n\lesssim 4$ the cut-off
scale turns out to be $\Lambda\gtrsim 10^{10}\leavevmode\nobreak\ \rm GeV$
that can well explain observed DM abundance for suitable choice of the
temperature $T_{R}\sim\mathcal{O}(100)\leavevmode\nobreak\ \rm MeV$. Thus, a
simple framework of the kind discussed here connects reheat temperature of the
universe with a SFOPT in a fast expanding cosmological background. A possible
underlying UV complete theory which can give rise to such dimension-5
effective interactions have also been discussed briefly. It will however be
interesting, in the context of a UV complete theory, to look into interplay of
different parameters giving rise to DM relic abundance and a stochastic GW
signal in a fast expanding universe, together with more exotic physics. We
shall address such issues in forthcoming works.

## Acknowledgements

BB would like to acknowledge discussions with Abhijit Kumar Saha, Debasish
Borah and e-mail conversations Sunny Vagnozzi. AP would like to thank Biswajit
Banerjee for the implementation of the CosmoTransition package. Work of ADB is
supported in part by the National Science Foundation of China (11775093,
11422545, 11947235).

## References

  * Arzoumanian  _et al._ (2020) Z. Arzoumanian _et al._ (NANOGrav), (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Sesana  _et al._ (2004) A. Sesana, F. Haardt, P. Madau, and M. Volonteri, Astrophys. J. 611, 623 (2004), arXiv:astro-ph/0401543 .
  * Ellis and Lewicki (2020) J. Ellis and M. Lewicki, (2020), arXiv:2009.06555 [astro-ph.CO] .
  * Blasi  _et al._ (2020) S. Blasi, V. Brdar, and K. Schmitz, (2020), arXiv:2009.06607 [astro-ph.CO] .
  * Chigusa  _et al._ (2020) S. Chigusa, Y. Nakai, and J. Zheng, (2020), arXiv:2011.04090 [hep-ph] .
  * Neronov  _et al._ (2020) A. Neronov, A. Roper Pol, C. Caprini, and D. Semikoz, (2020), arXiv:2009.14174 [astro-ph.CO] .
  * Nakai  _et al._ (2020) Y. Nakai, M. Suzuki, F. Takahashi, and M. Yamada, (2020), arXiv:2009.09754 [astro-ph.CO] .
  * Addazi  _et al._ (2020) A. Addazi, Y.-F. Cai, Q. Gan, A. Marciano, and K. Zeng, (2020), arXiv:2009.10327 [hep-ph] .
  * Bhattacharya  _et al._ (2020a) S. Bhattacharya, S. Mohanty, and P. Parashari, (2020a), arXiv:2010.05071 [astro-ph.CO] .
  * Kolb and Turner (1990) E. W. Kolb and M. S. Turner, Front. Phys. 69, 1 (1990).
  * D’Eramo  _et al._ (2017) F. D’Eramo, N. Fernandez, and S. Profumo, JCAP 05, 012 (2017), arXiv:1703.04793 [hep-ph] .
  * D’Eramo  _et al._ (2018) F. D’Eramo, N. Fernandez, and S. Profumo, JCAP 02, 046 (2018), arXiv:1712.07453 [hep-ph] .
  * Elahi  _et al._ (2015) F. Elahi, C. Kolda, and J. Unwin, JHEP 03, 048 (2015), arXiv:1410.6157 [hep-ph] .
  * Chen and Kang (2018) S.-L. Chen and Z. Kang, JCAP 05, 036 (2018), arXiv:1711.02556 [hep-ph] .
  * Bernal  _et al._ (2019) N. Bernal, F. Elahi, C. Maldonado, and J. Unwin, JCAP 1911, 026 (2019), arXiv:1909.07992 [hep-ph] .
  * Ade  _et al._ (2014) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Aghanim  _et al._ (2018) N. Aghanim _et al._ (Planck), (2018), arXiv:1807.06209 [astro-ph.CO] .
  * Bernal and Hajkarim (2019) N. Bernal and F. Hajkarim, Phys. Rev. D 100, 063502 (2019), arXiv:1905.10410 [astro-ph.CO] .
  * Bhattacharya  _et al._ (2020b) S. Bhattacharya, S. Mohanty, and P. Parashari, Phys. Rev. D 102, 043522 (2020b), arXiv:1912.01653 [astro-ph.CO] .
  * Bernal  _et al._ (2020) N. Bernal, A. Ghoshal, F. Hajkarim, and G. Lambiase, (2020), arXiv:2008.04959 [gr-qc] .
  * Ratra and Peebles (1988) B. Ratra and P. J. E. Peebles, Phys. Rev. D 37, 3406 (1988).
  * Hall  _et al._ (2010) L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, JHEP 03, 080 (2010), arXiv:0911.1120 [hep-ph] .
  * Bernal  _et al._ (2017) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, Int. J. Mod. Phys. A32, 1730023 (2017), arXiv:1706.07442 [hep-ph] .
  * Edsjo and Gondolo (1997) J. Edsjo and P. Gondolo, Phys. Rev. D 56, 1879 (1997), arXiv:hep-ph/9704361 .
  * Barman  _et al._ (2020a) B. Barman, S. Bhattacharya, and B. Grzadkowski, (2020a), arXiv:2009.07438 [hep-ph] .
  * De Simone  _et al._ (2010) A. De Simone, V. Sanz, and H. P. Sato, Phys. Rev. Lett. 105, 121802 (2010), arXiv:1004.1567 [hep-ph] .
  * Chung  _et al._ (1999) D. J. Chung, E. W. Kolb, and A. Riotto, Phys. Rev. D 60, 063504 (1999), arXiv:hep-ph/9809453 .
  * Moroi  _et al._ (1993) T. Moroi, H. Murayama, and M. Yamaguchi, Phys. Lett. B 303, 289 (1993).
  * Kawasaki and Moroi (1995) M. Kawasaki and T. Moroi, Prog. Theor. Phys. 93, 879 (1995), arXiv:hep-ph/9403364 .
  * Kofman  _et al._ (1997) L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys. Rev. D 56, 3258 (1997), arXiv:hep-ph/9704452 .
  * Linde (1990) A. D. Linde,  _Particle physics and inflationary cosmology_ , Vol. 5 (1990) arXiv:hep-th/0503203 .
  * Barman  _et al._ (2020b) B. Barman, D. Borah, and R. Roshan, JCAP 11, 021 (2020b), arXiv:2007.08768 [hep-ph] .
  * Kakizaki  _et al._ (2015) M. Kakizaki, S. Kanemura, and T. Matsui, Phys. Rev. D 92, 115007 (2015), arXiv:1509.08394 [hep-ph] .
  * Wainwright (2012) C. L. Wainwright, Comput. Phys. Commun. 183, 2006 (2012), arXiv:1109.4189 [hep-ph] .
  * Coleman and Weinberg (1973) S. R. Coleman and E. J. Weinberg, Phys. Rev. D7, 1888 (1973).
  * Basler  _et al._ (2017) P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt, and A. Wlotzka, JHEP 02, 121 (2017), arXiv:1612.04086 [hep-ph] .
  * Arnold and Espinosa (1993) P. B. Arnold and O. Espinosa, Phys. Rev. D47, 3546 (1993), [Erratum: Phys. Rev.D50,6662(1994)], arXiv:hep-ph/9212235 [hep-ph] .
  * Kosowsky  _et al._ (1992a) A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. D45, 4514 (1992a).
  * Paul  _et al._ (2019) A. Paul, B. Banerjee, and D. Majumdar, JCAP 1910, 062 (2019), arXiv:1908.00829 [hep-ph] .
  * Kosowsky and Turner (1993) A. Kosowsky and M. S. Turner, Phys. Rev. D47, 4372 (1993), arXiv:astro-ph/9211004 [astro-ph] .
  * Huber and Konstandin (2008) S. J. Huber and T. Konstandin, JCAP 0809, 022 (2008), arXiv:0806.1828 [hep-ph] .
  * Kosowsky  _et al._ (1992b) A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. Lett. 69, 2026 (1992b).
  * Kamionkowski  _et al._ (1994) M. Kamionkowski, A. Kosowsky, and M. S. Turner, Phys. Rev. D49, 2837 (1994), arXiv:astro-ph/9310044 [astro-ph] .
  * Caprini  _et al._ (2008) C. Caprini, R. Durrer, and G. Servant, Phys. Rev. D77, 124015 (2008), arXiv:0711.2593 [astro-ph] .
  * Hindmarsh  _et al._ (2014) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014), arXiv:1304.2433 [hep-ph] .
  * Giblin and Mertens (2013) J. T. Giblin, Jr. and J. B. Mertens, JHEP 12, 042 (2013), arXiv:1310.2948 [hep-th] .
  * Giblin and Mertens (2014) J. T. Giblin and J. B. Mertens, Phys. Rev. D90, 023532 (2014), arXiv:1405.4005 [astro-ph.CO] .
  * Hindmarsh  _et al._ (2015) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D92, 123009 (2015), arXiv:1504.03291 [astro-ph.CO] .
  * Caprini and Durrer (2006) C. Caprini and R. Durrer, Phys. Rev. D74, 063521 (2006), arXiv:astro-ph/0603476 [astro-ph] .
  * Kahniashvili  _et al._ (2008a) T. Kahniashvili, A. Kosowsky, G. Gogoberidze, and Y. Maravin, Phys. Rev. D78, 043003 (2008a), arXiv:0806.0293 [astro-ph] .
  * Kahniashvili  _et al._ (2008b) T. Kahniashvili, L. Campanelli, G. Gogoberidze, Y. Maravin, and B. Ratra, Phys. Rev. D78, 123006 (2008b), [Erratum: Phys. Rev.D79,109901(2009)], arXiv:0809.1899 [astro-ph] .
  * Kahniashvili  _et al._ (2010) T. Kahniashvili, L. Kisslinger, and T. Stevens, Phys. Rev. D81, 023004 (2010), arXiv:0905.0643 [astro-ph.CO] .
  * Caprini  _et al._ (2009) C. Caprini, R. Durrer, and G. Servant, JCAP 0912, 024 (2009), arXiv:0909.0622 [astro-ph.CO] .
  * Barman  _et al._ (2020c) B. Barman, A. Dutta Banik, and A. Paul, Phys. Rev. D 101, 055028 (2020c), arXiv:1912.12899 [hep-ph] .
  * Paul  _et al._ (2020) A. Paul, U. Mukhopadhyay, and D. Majumdar, (2020), arXiv:2010.03439 [hep-ph] .
  * Steinhardt (1982) P. J. Steinhardt, Phys. Rev. D25, 2074 (1982).
  * Dev  _et al._ (2019) P. S. B. Dev, F. Ferrer, Y. Zhang, and Y. Zhang, JCAP 1911, 006 (2019), arXiv:1905.00891 [hep-ph] .
  * Shajiee and Tofighi (2019) V. R. Shajiee and A. Tofighi, Eur. Phys. J. C79, 360 (2019), arXiv:1811.09807 [hep-ph] .
  * Caprini  _et al._ (2016) C. Caprini _et al._ , JCAP 1604, 001 (2016), arXiv:1512.06239 [astro-ph.CO] .
  * Ellis  _et al._ (2018) J. Ellis, M. Lewicki, and J. M. No,  (2018), 10.1088/1475-7516/2019/04/003, [JCAP1904,003(2019)], arXiv:1809.08242 [hep-ph] .
  * Ellis  _et al._ (2019) J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, JCAP 1906, 024 (2019), arXiv:1903.09642 [hep-ph] .
  * Thrane and Romano (2013) E. Thrane and J. D. Romano, Phys. Rev. D88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Schmitz (2020) K. Schmitz, (2020), arXiv:2002.04615 [hep-ph] .
  * Moore  _et al._ (2015) C. Moore, R. Cole, and C. Berry, Class. Quant. Grav. 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * Alanne  _et al._ (2020) T. Alanne, T. Hugle, M. Platscher, and K. Schmitz, JHEP 03, 004 (2020), arXiv:1909.11356 [hep-ph] .
  * Vagnozzi (2020) S. Vagnozzi, (2020), arXiv:2009.13432 [astro-ph.CO] .
  * Dutta Banik  _et al._ (2018) A. Dutta Banik, A. K. Saha, and A. Sil, Phys. Rev. D 98, 075013 (2018), arXiv:1806.08080 [hep-ph] .
  * Konar  _et al._ (2020) P. Konar, A. Mukherjee, A. K. Saha, and S. Show, (2020), arXiv:2007.15608 [hep-ph] .
  * Fairbairn and Hogan (2013) M. Fairbairn and R. Hogan, JHEP 09, 022 (2013), arXiv:1305.3452 [hep-ph] .


# Sandwiching the Riemann hypothesis

R.C. McPhedran,
School of Physics,
University of Sydney

###### Abstract

We consider a system of three analytic functions, two of which are known to
have all their zeros on the critical line $\Re(s)=\sigma=1/2$. We construct
inequalities which constrain the third function, $\xi(s)$, on $\Im(s)=0$ to
lie between the other two functions, in a sandwich structure. We investigate
what can be said about the location of zeros and radius of convergence of
expansions of $\xi(s)$, with promising results.

##  1 Introduction

The Riemann hypothesis has a long history and important consequences in number
theory, which are treated in numerous textbooks e.g. [1, 2]. Riemann’s
hypothesis was that all the non-trivial zeros of the Riemann zeta function
$\zeta(s)$ lie on the critical line $\Re(s)=1/2$. There have been many methods
developed in order to try to prove this result by analytic reasoning and by
impressive computer evaluations (with the hypothesis verified numerically up
to $3\times 10^{12}$ [3]).

Among the many threads followed in the literature, we concentrate here on one
based on sums over inverse powers of the zeros of $\zeta(s)$ on the critical
line [4, 5, 6]. This is chosen because it permits a combination of numerical
exploration and analytic reasoning, with significant examples of both in the
literature. The focus of this work will be on the Keiper-Li theory, developed
independently by J. B. Keiper and X. J. Li. This is based on expansions of the
symmetrized counterpart of $\zeta(s)$, $\xi(s)$. The expansions involve two
sets of constants, $a_{n}$ and $\lambda(n)$, the first due to Li, and the
second to both authors. Li proved that the Riemann hypothesis holds if and
only if all $\lambda(n)$ are non-negative. The author and colleagues have
recently developed an expression for the positive constants $a_{n}$, and
evaluated to good accuracy the first 4000 of them [7, 8].

We will combine the Keiper-Li framework with results from the literature
concerning two functions, denoted here as $\xi_{+}(s)$ and $\xi_{-}(s)$. These
consist of combinations of $\xi(s+1/2)$ and $\xi(s-1/2)$, and it is known that
all their zeros lie on the critical line and are simple [9, 10, 11, 12]. We
will consider the properties of $\xi_{+}(s)$ and $\xi_{-}(s)$, and show that
along the positive real axis $\xi(s)$ is constrained to lie between
$\xi_{+}(s)$ and $\xi_{-}(s)$, in what we will call a sandwich. We will reason
from the properties of the sandwich that, under a mapping of the complex
variable $s=1/2+it$ onto the unit circle, the radius of convergence of the
expansion of $\xi(s)$ is unity. This result is highly significant with respect
to the Riemann hypothesis.

##  2 Basic equations and some results

We will be interested here in some functions related to the Riemann $\xi$
function:

$\xi(s)=\frac{1}{2}s(s-1)\frac{\Gamma(s/2)\zeta(s)}{\pi^{s/2}}.$ (1)

(Note that Li [6] uses a definition of $\xi(s)$ which is two times larger than
that used here.) This is an entire function, with $\log\xi(s)$ having
logarithmic singularities at the roots $\rho$ of $\xi(s)$ and no other
singularities [2]. The particular functions of interest are $\xi(s)$,
$\xi(s+1/2)$, $\xi(s-1/2)$ and two combinations of the latter two:

$\xi_{-}(s)=\xi(s+1/2)-\xi(s-1/2),~{}~{}\xi_{+}(s)=\xi(s+1/2)+\xi(s-1/2).$ (2)

The functions $\xi_{-}(s)$ and $\xi_{+}(s)$ are notable in that it has been
proved that these functions obey the Riemann hypothesis: all their zeros
$\rho$ lie on the critical line $\Re(s)=\sigma=1/2$. The functions were
studied by Taylor [9], Lagarias and Suzuki [11] and Ki[12]. The fact that all
the non-trivial zeros of the $\xi_{-}(s)$ combination lie on the critical line
was first established by P.R. Taylor, and published posthumously. Lagarias and
Suzuki considered the $\xi_{+}(s)$ combination, and showed that all its
complex zeros lie on the critical line, while Ki proved that all the complex
zeros were simple. A further useful property is that the complex zeros of the
$\xi_{+}(s)$ and $\xi_{-}(s)$ combinations strictly alternate on the critical
line, and have the same distribution function of zeros. The common
distribution function is indeed that corresponding to any prescribed argument
value of $\xi(s)$ on the line $\sigma=1$.

Lagarias and Suzuki [11] in their Theorem 4 present a general result relating
to functions consisting of a superposition of two identical parts, each of the
parts having zeros symmetrically placed about a midline in a critical strip,
and the two parts separated sufficiently so their critical strips do not
overlap. Then the symmetry of the superposition guarantees that the moduli of
the two parts can only be equal halfway between the two separated critical
regions. In turn, if the superposition of the two parts is a sum with a phase
factor of unit modulus, that sum will have all its zeros on the midline of the
total system. This result then includes Taylor’s result for $\xi_{-}(s)$ as a
special case, along with that for $\xi_{+}(s)$. Related previous work is [1,
13].

We investigate the properties of the five $\xi$ functions using power series
known to have coefficients all of which are known to be positive. The first
has been studied intensively by Pustyl’nikov [14, 15] and subsequent authors.
The set of coefficients $\xi_{r}$ occurs in the following expansion:

$\xi(s+1/2)=\sum_{r=0}^{\infty}\xi_{r}s^{2r},$ (3)

where $\xi_{0}=\xi(1/2)\approx 0.49712077818831410991$. Pustyl’nikov proved
that all the $\xi_{r}$ are positive and that this is a necessary condition for
the Riemann hypothesis to hold. Accurate numerical techniques for the
evaluation of the $\xi_{r}$ have been developed in subsequent work- see for
example [16]. A comprehensive and accurate tabulation of values due to Dr.
Rick Kreminski [17] was previously available on the Internet, but it appears
to be currently inaccessible.

The expansion (3) can be re-expressed to give

$\xi_{+}(s)=2\sum_{n=0}^{\infty}\xi_{n}s^{2n}-\sum_{n=0}^{\infty}\sum_{r=n+1}^{\infty}\left[s\binom{2r}{2n+1}-\binom{2r}{2n}\right]\xi_{r}s^{2n},$
(4)

and

$\xi_{-}(s)=\sum_{n=0}^{\infty}\sum_{r=n+1}^{\infty}\left[s\binom{2r}{2n+1}-\binom{2r}{2n}\right]\xi_{r}s^{2n}.$
(5)

These are respectively even and odd under the substitution $s\rightarrow 1-s$.
Adding and subtracting (4) and (5) gives the expansion (3) for $\xi(s+1/2)$
and the corresponding equation for $\xi(s-1/2)$.

The corresponding result for $2\xi(s)$ also contains both odd and even powers
of $s$ and can be written

$2\xi(s)=1+\sum_{l=1}^{\infty}{\cal E}_{l}s^{2l}-\sum_{l=1}^{\infty}{\cal
O}_{l}s^{2l-1},$ (6)

where

${\cal E}_{l}=2\sum_{r=l}^{\infty}{2r\choose
2l}\frac{\xi_{r}}{2^{2r-2l}},~{}{\cal O}_{l}=2\sum_{r=l}^{\infty}{2r\choose
2l-1}\frac{\xi_{r}}{2^{2r-2l+1}}.$ (7)

The first ten coefficients ${\cal E}$ of even order are

$\displaystyle
0.023343864534226183135,0.00025318173031652700506,1.7209870418615355778*10^{-6},$
$\displaystyle
8.3159682500277216307*10^{-9},3.0655602327633313510*10^{-11},9.0229664497612087603*10^{-14},$
$\displaystyle
2.1893251340686846583*10^{-16},4.4843405072454944930*10^{-19},7.8974339566658717737*10^{-22},$
$\displaystyle 1.2134779622875435114*10^{-24}.$ (8)

The first ten coefficients ${\cal O}$ of odd order are

$\displaystyle
0.023095708966121033814,0.00049798384992294867235,5.0502547922191741696*10^{-6},$
$\displaystyle 3.2378414618810769603*10^{-8},\
1.4852419214918940045*10^{-10},5.2238407222768796166*10^{-13},$ $\displaystyle
1.4729420831622495667*10^{-15},3.4352126539793423994*10^{-18},6.7821598786771781572*10^{-21},$
$\displaystyle 1.1540437076606000624*10^{-23}.$ (9)

The second set of power series is based around the Li coefficients $a_{n}$
[6]. These occur in the expansion

$2\xi\left(\frac{1}{1-z}\right)=\phi(z)=1+\sum_{j=1}^{\infty}a_{j}z^{j},$ (10)

for $z$ in the unit disc. An equivalent form for $w$ to the right of the
critical line is obtained using the substitution $1-z=1/w$:

$2\xi(w)=\phi\left(\frac{w-1}{w}\right)=1+\sum_{j=1}^{\infty}a_{j}\left(\frac{w-1}{w}\right)^{j}=\phi_{1}(w).$
(11)

To the left of the critical line, using the symmetry property
$\xi(s)=\xi(1-s)$, we find:

$2\xi(w)=\phi\left(\frac{w}{w-1}\right)=1+\sum_{j=1}^{\infty}a_{j}\left(\frac{w}{w-1}\right)^{j}=\phi_{2}(w).$
(12)

Note that $\phi_{1}(1)=1$ exactly, and $\phi_{2}(0)=1$ exactly.

Li’s paper [6] is important in that it establishes as a necessary condition
for the Riemann hypothesis to hold that the logarithmic derivative of the
function $\xi(1/(1-z))$ be analytic in the unit disc: see also Bombieri and
Lagarias [10] and Coffey [18]. Important related work by Keiper [5] predated
that of Li. Note that the following expression has been established for the
$a_{n}$ [7, 8]:

$a_{n}=2\sum_{p=1}^{n}{\cal C}_{n,p}\Sigma_{p}^{\xi},$ (13)

where the $\Sigma_{p}^{\xi}$ are positive-valued sums involving the $\xi_{r}$:

$\Sigma_{p}^{\xi}=\sum_{r=1}^{\infty}\frac{\xi_{r}}{2^{2r}}r^{p},$ (14)

and $\Sigma_{0}^{\xi}=1/2-\xi_{0}$. The $C_{n,p}$ arise in coefficients of
polynomials $a_{r}(n)$ generated by the expansion of a quotient function:

$\left(\frac{1+w}{1-w}\right)^{r}=1+\sum_{n=1}^{\infty}a_{r}(n)w^{n}.$ (15)

All the $C_{n,p}$ are non-negative, being zero if the pair $n,p$ mixes an even
and an odd integer, and being positive if both are even or odd. Note that
$C_{n,n}=4^{n}/n!$ with $C_{n,p}$ tending to zero as $n$ increases. The
$C_{n,p}$ obey a simple recurrence relation [7, 8], from which they may be
determined exactly.

The equation (13) has been used to obtain the first 4000 $a_{n}$ values to 200
decimal places accuracy [7, 8], and they may be downloaded. The first 500
$\xi_{r}$ values were sufficient to supply this accuracy. Note that

$a_{1}=\lambda_{1}=\sigma_{1}=1+\gamma/2-\log(4\pi)/2=8\Sigma_{1}^{\xi}.$ (16)

Here $a_{1}\approx 0.023095708966121033814$.

The ${\cal C}_{n,p}$ obey the useful relation:

$\sum_{p=1}^{n}{\cal C}_{n,p}=4n.$ (17)

Also, from equation (14) the $\Sigma_{p}^{\xi}$ increase monotonically with
$p$. Hence, the $a_{n}$ for $n>1$ satisfy the bounds

$na_{1}<a_{n}<8n\Sigma_{n}^{\xi}.$ (18)

From [7] $\log\Sigma_{n}^{\xi}$ has the leading terms for $n$ large

$n[\log n-\log\log n-2]$ (19)

so the lower and upper bounds in (18) differ considerably in their growth
rates as $n$ increases.

One important consequence of the inequalities (18) is that

$\sum_{n=1}^{N}a_{n}>a_{1}\sum_{n=1}^{N}n=a_{1}\frac{N(N+1)}{2},$ (20)

which diverges quadratically as $N\rightarrow\infty$. A second is that the
series $\phi_{1}(w)$ must diverge if $|(w-1)/w|>1$, i.e. if $Re(w)<1/2$ and
$|z|>1$, while the series $\phi_{2}(w)$ must diverge if $|w/(w-1)|>1$, i.e. if
$Re(w)>1/2$ and $|z|<1$.

##  3 The first two sandwiches

The equation (11) can be used to construct inequalities and monotonicity
properties among the key functions introduced in the previous section.
Consider the difference $\xi(u+a)-\xi(u+b)$ where $u$,$a$ and $b$ are real
with $a>b$ and $u+a$ and $u+b$ both exceeding 1/2 :

$2\xi(u+a)-2\xi(u+b)=\sum_{j=1}^{\infty}a_{j}\left[\left(1-\frac{1}{u+a}\right)^{j}-\left(1-\frac{1}{u+b}\right)^{j}\right]>0,$
(21)

since every element of the summand is positive. Hence, $\xi(u+1/2)>\xi(u)$ and
$\xi(u)>\xi(u-1/2)$. For $\xi(u-1/2)$ we have

$2\xi(u-1/2)=1+\sum_{j=1}^{\infty}a_{j}\left(1-\frac{1}{u-1/2}\right)^{j}>1,~{}{\rm
for}~{}u>3/2.$ (22)

The full set of inequalities for the first three functions of interest is then

$\xi\left(u+\frac{1}{2}\right)>\xi(u)>\xi\left(u-\frac{1}{2}\right)>0.$ (23)

Each of the three functions increases monotonically as the argument $u$
increases, and also if the full sum is replaced by a partial sum, the partial
sums increase monotonically with the upper limit on $j$.

We now consider mappings which move the location of lines along which zeros
are located onto circles in the complex plane. For the case of $\xi(s)$ a
convenient mapping from the critical line $\Re(s)=\sigma=1/2$ onto the unit
circle has already been given:

$w=u+iv=1-\frac{1}{s}=\frac{s-1}{s}.$ (24)

The inverse transformation is

$s=\frac{1}{1-w}.$ (25)

For the function $\xi\left(s+1/2\right)$ the corresponding forward
transformation is

$w_{h}=1-\frac{1}{s+1/2}=\frac{s-1/2}{s+1/2}.$ (26)

Its inverse transformation is

$s=-\frac{1}{2}+\frac{1}{1-w_{h}}.$ (27)

The known zeros of $\xi\left(s+1/2\right)$ lie on $\sigma=0$ and are mapped
onto the unit circle in the plane of complex $w_{h}$.

For the function $\xi\left(s-1/2\right)$ the corresponding forward
transformation is

$w_{m}=1-\frac{1}{s-1/2}=\frac{s-3/2}{s-1/2}.$ (28)

Its inverse transformation is

$s=\frac{1}{2}+\frac{1}{1-w_{m}}.$ (29)

The known zeros of $\xi\left(s-1/2\right)$ lie on $\sigma=1$ and are mapped
onto the unit circle in the plane of complex $w_{m}$.

We next connect $w_{h}$ and $w_{m}$ to the complex variable $w$. Eliminating
$s$ between (24) and (27) we find

$w_{h}=\frac{1+w}{3-w},~{}~{}w=\frac{3w_{h}-1}{w_{h}+1}.$ (30)

The equation for the fixed point of this transformation is

$w(3-w)=1+w,~{}{\rm or}~{}(w-1)^{2}=0.$ (31)

The fixed point is then of second order at $w=w_{h}=1$. The corresponding
equations relating to $w_{m}$ are

$w_{m}=\frac{3w-1}{w+1},~{}~{}w=\frac{1+w_{m}}{3-w_{m}},$ (32)

while again the fixed point for the transformation yields $w=1=w_{m}$ being of
second order.

Figure 1 illustrates curves of constant modulus in the complex $w$ plane
pertaining to these discussions. The black unit circle corresponds to $|w|=1$,
and thus to the mapping of the critical line through equation (24). The black
circle centred on $w=1/2$ of radius $1/2$ corresponds to the constraint
$|w_{m}|=1$, and the mapping of the line $\sigma=1$, with the black vertical
line being $|w_{m}|=3$, $u=-1/3$. The red lines are for $|w_{h}|=1$ and
$|w_{h}|=2$. The two fixed points occur where the two circles touch, with the
line $|w_{h}|=1$ being tangent to both.

Figure 1:  The transformations pertaining to $w$ (blue circle), $w_{m}$ (black
circle and black line) and $w_{h}$ (red lines).

The behaviour of the three functions $\xi(s+1/2)$, $\xi(s)$ and $\xi(s-1/2)$
along a part of the real axis is shown in Fig. 2, and is in accord with the
inequalities (23) . The figure shows the first of the two sandwiches we will
study. Our aim will be to understand the location of the singularities of the
three functions, as reflected in their power series and their behaviour along
the real axis.

Figure 2: The behaviour of the three functions $\xi(s+1/2)$ (red), $\xi(s)$
(blue) and $\xi(s-1/2)$ (black) along part of the real axis is shown.

Titchmarsh [19] gives two useful results concerning the radius of the circle
of convergence of power series with coefficients $a_{n}$ in the variable $z$
and their behaviour along the real axis. For convenience the radius of the
circle of convergence is normalised to be unity. The first is that if
$a_{n}\geq 0$ for all $n$, then $z=1$ is a singular point. The second is that
if $a_{n}$ is real for all values of $n$ and $\sum a_{n}$ is properly
divergent, then $z=1$ is a singular point. We note that both criteria apply to
power series with the coefficients which are the Li $a_{n}$ (see equation
(20)). If the coefficients are the $\xi_{n}$, the first will apply.

We thus deduce that the circle of convergence of $\phi(z)$ cuts the positive
real axis at a singular point of its power series (10). The same conclusion
applies to the circles of convergence pertaining to $\phi(z_{h})$ where
$z_{h}=z+1/2$ and to $\phi(z_{m})$ where $z_{m}=z-1/2$.

The behaviour of the three functions $\xi_{+}(s)$, $\xi(s)$ and $\xi_{-}(s)$
along an arc of the unit circle is shown in Fig. 3. Specifically, the function
plotted is $\log\xi[1/(1-\exp(i\theta)]$, where $\theta$ denotes the angular
position on the unit circle. The dips in the graph indicate zeros of the
function. For $\xi_{+}$ they occur where $\arg[\xi(1+(i/2)\cot(\theta/2))]$ is
an odd multiple of $\pi/2$, and for $\xi_{+}$ where it is an even multiple of
$\pi/2$. Note that zeros of $\xi_{+}(s)$ occur closer to those of $\xi(s)$
than is the case for zeros of $\xi_{-}(s)$, and that zeros of $\xi_{+}(s)$
strictly alternate with those of $\xi_{-}(s)$. (This is a consequence of the
monotonic variation of $\arg[\xi(1+(i/2)t]$ with $t$.)

Figure 3: The behaviour of the three functions $\xi_{+}(s)$ (red), $\xi(s)$
(blue) and $\xi_{-}(s)$ (black) along part of unit circle is shown.

The behaviour of $\xi_{+}(s)$, $\xi(s)$ and $\xi_{-}(s)$ along a part of the
real axis is shown in Fig. 4. The figure gives evidence for the following
inequalities:

$\xi_{+}(\sigma)>\xi(\sigma)>\xi_{-}(\sigma)>0~{}{\rm for}~{}\sigma\geq 1.$
(33)

These inequalities follow immediately from (23), with the exception of
$\xi(\sigma)>\xi_{-}(\sigma)$. This last follows from the asymptotic
expression:

$\displaystyle 2\xi(\sigma)-2[\xi(\sigma+1/2)-\xi(\sigma-1/2)]$
$\displaystyle=$ $\displaystyle
1+\sum_{j}a_{j}\left(1-\frac{1}{\sigma}\right)^{j}$
$\displaystyle\left\\{1+\exp\left[j\left(\log\left(1-\frac{1}{(\sigma-1/2)}\right)-\log\left(1-\frac{1}{\sigma}\right)\right)\right]\right.$
$\displaystyle\left.-\exp\left[j\left(\log\left(1-\frac{1}{(\sigma+1/2)}\right)-\log\left(1-\frac{1}{\sigma}\right)\right)\right]\right\\}$

giving for $j<<\sigma^{2}$:

$\displaystyle
2\xi(\sigma)-2[\xi(\sigma+1/2)-\xi(\sigma-1/2)]=1+\sum_{j}a_{j}\left(1-\frac{1}{\sigma}\right)^{j}\left[1-\frac{j}{\sigma^{2}}-\frac{j}{\sigma^{3}}+O\left(\left(\frac{1}{\sigma}\right)^{4}\right)\right]$
(35)

This is increasing with $\sigma$, diverging as $\sigma\rightarrow\infty$ by
virtue of the increasing values of $a_{j}$. The functions $\xi_{+}(\sigma)$,
$\xi(\sigma)$ and $\xi_{-}(\sigma)$ also increase monotonically with
$\sigma>1$ and are all positive for $\sigma>1/2$ . They must all diverge for
$\sigma\rightarrow\infty$.

Figure 4: The behaviour of the three functions $\xi_{+}(s)$ (red), $\xi(s)$
(blue) and $\xi_{-}(s)$ (black) along part of the real axis is shown.

##  4 The third sandwich

We next consider the logarithms of the three functions $\xi_{+}(s)$, $\xi(s)$
and $\xi_{-}(s)$, as shown in Fig. 5. The curves there illustrate a similar
behaviour to previous figures for $\xi_{+}(s)$ and $\xi(s)$. For $\xi_{-}(s)$
there is a notable influence of the first order zero at $s=1/2$. The monotonic
behaviour of the logarithm enables us to go from the inequalities (33) to

$\log\xi_{+}(\sigma)>\log\xi(\sigma)>\log\xi_{-}(\sigma)>0~{}{\rm
for}~{}\sigma\geq 1.$ (36)

The monotonic increase of the three functions of (36) with $\sigma$ is also
preserved.

Figure 5: The behaviour of the three functions $\log(\xi_{+}(s))$ (red),
$\log(\xi(s))$ (blue) and $\Re\log(\xi_{-}(s))$ (black) along part of the real
axis is shown.

A power series for $\log\xi(s)$ with coefficients depending on the equation
(6) can be constructed in two ways. The first is to use Mathematica or a
similar symbolic package. The second is to differentiate the series for
$\xi(s)$ and divide the result by $\xi(s)$, before integrating the result for
the quotient. The result is the same, to within slight numerical differences:

$\displaystyle\log(\xi(s))=-0.023095708966121033814s+0.023077158647902301379s^{2}+0.0000370527438173686409s^{3}$
$\displaystyle-0.00001840680531542237958s^{4}-1.43018671152521547*10^{-7}s^{5}+4.6906069489794377*10^{-8}s^{6}$
$\displaystyle+6.534558785292532*10^{-10}s^{7}-1.5860851386884509*10^{-10}s^{8}-3.141596838950986*10^{-12}s^{9}$
$\displaystyle+5.99771484715187*10^{-13}s^{10}+1.53698978205383*10^{-14}s^{11}-2.39484594174383*10^{-15}s^{12}$
$\displaystyle-7.5638205856654*10^{-17}s^{13}+9.851878115521*10^{-18}s^{14}+3.727558439514*10^{-19}s^{15}$
$\displaystyle-4.123828024852*10^{-20}s^{1}6-1.836227510375*10^{-21}s^{17}+1.743694268780*10^{-22}s^{18}$
$\displaystyle+9.03450466141*10^{-24}s^{19}-7.4119966493*10^{-25}s^{20}$ (37)

This is a power series with coefficients of mixed sign, unlike that for
$\xi(s)$. It has been obtained using Mathematica, by the first method, with 20
decimal places requested as the accuracy goal. Note that higher order
coefficients have fewer decimals than specified in the accuracy goal.

The main focus of this section will be a discussion of the radius of
convergence of the expansion (10), and its equivalent for
$\log[2\xi(1/(1-z))]$. From [19], the radius of convergence associated with
the series in (10) is

$R=\lim_{n\rightarrow\infty}R_{n}=\lim_{n\rightarrow\infty}\exp\left[\frac{-\log
a_{n}}{n}\right].$ (38)

In Fig. 6 we show the variation of $\log a_{n}/n$ with $n$, with a rapid rise
for values of $n$ up to around 300 being succeeded by a slow fall off for
higher values. Also shown are the associated values of $R_{n}$ from equation
(38) and a numerical fit for the range of $n$ from 1000 to 4000:

$\frac{\log a_{n}}{n}\approx 0.11652745618-0.01092578334\log n.$ (39)

The simple fit function gives a reasonably accurate representation of the
variation of $\log a_{n}/n$ in the range shown, and confirms that its leading
varying term is logarithmic. Values of $a_{n}$ for $n$ much larger than 4000
would be needed to get additional terms in the expression for $\log a_{n}/n$
and thus for $R_{n}$.

Figure 6: (Left) The values of $\log a_{n}/n$ for $n$ ranging up to 4000.
(Right) The $R_{n}$ values (blue) and the fit function (39) (red) for $n$
ranging from 1000 to 4000.

Titchmarsh [19] gives a necessary and sufficient condition that a power series
with coefficients $a_{n}$ (known here to be non-negative) should have a
singularity at the point $z=1$ lying on its circle of convergence. Let $b_{n}$
denote the quantities

$b_{n}=\sum_{m=0}^{n}{n\choose m}a_{m}.$ (40)

Then the criterion for a singularity to occur at $z=1$ is

$\lim_{n\rightarrow\infty}b_{n}^{-\frac{1}{n}}=\frac{1}{2}.$ (41)

The application of this criterion to the power series with coefficients
$a_{n}$ is shown in Fig. 7. While both the radius of convergence estimate and
the singularity estimate converge slowly with increasing $n$, they are
confirmatory of there being a singularity at $z=1$ lying on the unit circle of
convergence.

Figure 7: The values of the singularity estimate $b_{n}$ as a function of $n$.

Of more importance to the following discussion than the series for
$\log\xi(s)$ are those for $2\xi\left(\frac{1}{1-z}\right)=\phi(z)$ and
$\phi^{\prime}(z)/\phi(z)$. For the last, Li states that a necessary and
sufficient condition for the nontrivial zeros of $\zeta(s)$ to lie on the
critical line is that $\phi^{\prime}(z)/\phi(z)$ is analytic in the unit disc.
Furthermore, given that the series for $\phi^{\prime}(z)/\phi(z)$ around $z=0$
is

$\frac{\phi^{\prime}(z)}{\phi(z)}=\sum_{n=0}^{\infty}\lambda_{n+1}z^{n},$ (42)

the necessary and sufficient condition is equivalent to the requirement that
$\lambda_{n}\geq 0$ for every positive $n$. Keiper had previously shown that
non-negativity of the $\lambda_{n}$ was a necessary condition for the Riemann
hypothesis, and had calculated the $\lambda_{n}$ to high accuracy up to
$n=4000$. Note that Keiper based his definition of the $\lambda_{n}$ around
the expansion of $\log\phi(z)$, resulting in a factor of $n$ between his
$\lambda_{n}^{K}$ and those of Li ($\lambda_{n}^{L}$).

One way of going from equation (10) to the required series for
$\phi^{\prime}(z)/\phi(z)$ is to form the series for $1/\phi(z)$:

$\frac{1}{\phi(z)}=1+\sum_{j=1}^{\infty}A_{j}z^{j},$ (43)

via the recurrence relation

$A_{1}=-a_{1},~{}A_{j}=-a_{j}-\sum_{p=1}^{j-1}a_{p}A_{j-p}.$ (44)

The required series then comes from the product:

$[\sum_{j=1}^{J}ja_{j}z^{j-1}]\times[1+\sum_{j=1}^{J}A_{j}z^{j}],$ (45)

where $J$ specifies the number of values for $\lambda$’s occurring in the
series.

The series for $1/\phi(z)$ is to order 20:

$\displaystyle
1-0.0230957089661210338143102479065z-0.0459061617276994534365358813998z^{2}$
$\displaystyle-0.0681486316594069122599620826562z^{3}-0.089545433048398089752378459144z^{4}$
$\displaystyle-0.109826396050622458253486711461z^{5}-0.128731287368373527211506921969z^{6}$
$\displaystyle-0.146012159089427308368403139905z^{7}-0.161435608581258419635214039431z^{8}$
$\displaystyle-0.174784932895325800969831596126z^{9}-0.185862161803051549940505866001z^{10}$
$\displaystyle-0.194489954353028235590838150534z^{11}-0.200513344703746739263814823922z^{12}$
$\displaystyle-0.203801323942739062262628361287z^{13}-0.20424824564594307240714200573z^{14}$
$\displaystyle-0.20177504405427703081208059797z^{15}-0.19633025494134585958756091470z^{16}$
$\displaystyle-0.18789083050993567984657371103z^{17}-0.17646274097813564381563835719z^{18}$
$\displaystyle-0.16208135689084447159310105477z^{19}-0.14481160761104406228773575035z^{20}$
$\displaystyle+O\left(z^{21}\right).$ (46)

Here, all coefficients of $z$ are negative. This results in a monotonic
decreasing function, as shown in Fig. 8. Note that the truncated series (46)
works well until $z$ approaches the exponential region near unity.

Figure 8: The function $1/\phi(z)$ as a function of $z$ is compared with the
series (46).

The series for $\phi^{\prime}(z)/\phi(z)$ is to order 19:

$\displaystyle\frac{\phi^{\prime}(z)}{\phi(z)}=\sum_{n=0}^{\infty}\lambda_{n+1}^{L}z^{n}\hskip
142.26378pt$
$\displaystyle=0.0230957089661210338143102479065+0.0923457352280466703857284861921z$
$\displaystyle+0.207638920554324803791492046618z^{2}+0.368790479492241638590511489638z^{3}$
$\displaystyle+0.575542714461177452431106405493z^{4}+0.82756601228237929742500282202z^{5}$
$\displaystyle+1.12446011757095949058282010802z^{6}+1.46575567714706063265551454198z^{7}$
$\displaystyle+1.85091604838253415532604486792z^{8}+2.27933936319315774369303405737z^{9}$
$\displaystyle+2.75036083822019606035454709285z^{10}+3.26325532062461984807908598991z^{11}$
$\displaystyle+3.81724005784794598710436795129z^{12}+4.4114776786805985120806412969z^{13}$
$\displaystyle+5.0450793720267934585351114375z^{14}+5.7171082488687926394190666698z^{15}$
$\displaystyle+6.4265828721172029011455409609z^{16}+7.1724809382917229592529707263z^{17}$
$\displaystyle+7.9537430943119003048082250779z^{18}+8.7692768720932151295994613534z^{19}$
$\displaystyle+O\left(z^{20}\right)$ (47)

All coefficients of $z$ being positive, this is a monotonically increasing
function, positive within the radius of convergence of the series. We noted at
the beginning of this section that the monotonic increase of the functions
$\xi_{+}(s)$, $\xi(s)$ and $\xi_{-}(s)$ for $s$ real and larger than unity was
preserved for $\log\xi_{+}(s)$, $\log\xi(s)$ and $\log\xi_{-}(s)$, whose
series incorporate the same coefficients divided by their power as do their
derivatives.

More definite information comes from the knowledge that $\xi_{+}(s)$ and
$\xi_{-}(s)$ have all their zeros on the critical line, mapped by
$s\rightarrow 1/(1-z)$ onto the boundary of the unit circle. The point $z=1$
is then a limit point for a sequence of zeros (see Fig. 3), and is thus an
essential singularity of the logarithms of these functions. For each of
$\log[\xi_{+}(\sigma)]$ and $\log[\xi_{-}(\sigma)]$, the function tends to
infinity as $\sigma\rightarrow 1$. By virtue of the inequalities (36), we know
then that $\log[\xi(\sigma)]$ tends to infinity as $\sigma\rightarrow 1$ (in
keeping with equation (10)).

Consider what would be the situation if the Riemann hypothesis failed. There
would then be at least one zero lying properly inside the unit circle, and a
singularity of $\log[\phi(z_{*})]$ at a point $z_{*}$ with $|z_{*}|<1$ (with
naturally a corresponding singularity occurring outside the unit circle). Take
$\rho_{*}$ to be the minimum of all such $|z_{*}|$. Then the circle of radius
$\rho_{*}$ centred on $z=0$ is the circle of convergence of the logarithm of
the expansion in (10). This means that the series diverges on any ray as the
modulus of $z$ approaches $\rho_{*}$. Applying this to the real axis, the
modulus of $\log[\phi(\sigma)]$ has to diverge to positive infinity at
$\sigma=\rho_{*}<1$, and to exceed the modulus of $\log[\xi_{+}(\sigma)]$
there, in contradiction with (36).

Another way of looking at this situation is that the sandwiching shows that
$\log[\phi(\sigma)]$ cannot diverge before $\log[\xi_{+}(\sigma)]$ diverges.
It also cannot diverge after $\log[\xi_{-}(\sigma)]$ diverges.

Both arguments lead to the conclusion that indeed the radius of convergence of
the function $\log[\phi(\sigma)]$ is unity, and that $\sigma=1$ is an
essential singularity.

## References

  * [1] Titchmarsh, E.C. & Heath-Brown, D.R., 1986 The Theory of the Riemann Zeta-function, Oxford, Oxford pp. 254–291.
  * [2] Edwards, H.M., Riemann’s Zeta Function, Dover, Mineola (2001) pp. 132–134.
  * [3] Platt, D. and Trudgian,T. 2020, The Riemann hypothesis is true up to $3\times 10^{12}$, arXiv:2004.09765v1.
  * [4] Lehmer, D.H. 1988 The sum of like powers of the zeros of the Riemann zeta function Math. Comp. 50 265–273.
  * [5] Keiper, J.B. 1992 Power Series Expansions of Riemann’s $\xi$ Function Math. Comp. 58 765-773.
  * [6] Li, X.J. 1997 The positivity of a sequence of numbers and the Riemann hypothesis J. Number Th. 65 325–333.
  * [7] McPhedran, R.C., Scott, T.C. and Maignan, A., 2022 Comments on and Extensions to Criteria of Keiper and Li for the Riemann Hypothesis HAL Archive, hal-03579652v1 and ACM Communications in Computer Algebra, 57, pp. 85-110 (2023).
  * [8] McPhedran, R.C. 2024 Numerical investigations of the Keiper-Li Criterion for the Riemann hypothesis arXiv:2311.06294, pp. 27.
  * [9] Taylor, P.R. 1945 On the Riemann zeta-function, Q.J.O., 16, 1-21.
  * [10] Bombieri, E. and Lagarias, J.C. 1999 Complements to Li’s Criterion for the Riemann Hypothesis J. Number Th. 77 274–287.
  * [11] Lagarias, J.C. and Suzuki, M. 2006 The Riemann hypothesis for certain integrals of Eisenstein series J. Number Theory 118 98–122.
  * [12] Ki, H. 2006 Zeros of the constant term in the Chowla-Selberg formula Acta Arithmetica 124 197-204.
  * [13] Pólya, G. 1926 Bemerkung uber die Integraldarstellung der Riemannschen xi-Funktion Acta Mathematica,48, 305–317.
  * [14] Pustyl’nikov, L.D. 1999 On a property of the classical zeta-function associated with the Riemann hypothesis Russian Mathematical Surveys 54 262–263.
  * [15] Pustyl’nikov, L.D. 2000 On the asymptotic behaviour of the Taylor series coefficients of $\xi(s)$ Russian Mathematical Surveys 55 349–350.
  * [16] Griffin, M., Ono,K., Rolen,L. and Zagier, D. Jensen polynomials for the Riemann zeta function and other sequences 2019 Proceedings of the National Academy of Sciences, 116, 11103-11110.
  * [17] Kreminski, R. 2005
http://faculty.colostate-
pueblo.edu/rick.kreminski/stieltjesrelated/sigmavalues2005.txt

  * [18] Coffey, M.W. 2005 Toward verification of the Riemann hypothesis: application of the Li criterion Mathematical Physics, Analysis and Geometry 8 211-255.
  * [19] Titchmarsh, E.C. 1939 The Theory of Functions, Oxford, Oxford, Chapter 7.