# Imprints of MeV Scale Hidden Dark Sector at Planck Data

Sougata Ganguly<EMAIL_ADDRESS>School of Physical Sciences, Indian
Association for the Cultivation of Science,
2A $\&$ 2B Raja S.C. Mullick Road, Kolkata 700032, India  Sourov Roy
<EMAIL_ADDRESS>School of Physical Sciences, Indian Association for the
Cultivation of Science,
2A $\&$ 2B Raja S.C. Mullick Road, Kolkata 700032, India  Abhijit Kumar Saha
<EMAIL_ADDRESS>School of Physical Sciences, Indian Association for the
Cultivation of Science,
2A $\&$ 2B Raja S.C. Mullick Road, Kolkata 700032, India

###### Abstract

New light species can contribute to the number of effective relativistic
degrees of freedom ($N_{\rm eff}$) at Cosmic Microwave Background (CMB) which
is precisely measured by Planck. In this work, we consider an MeV scale
thermally decoupled non-minimal dark sector and study the imprint of the dark
sector dynamics on the measurement of $N_{\rm eff}$ at the time of CMB
formation. We have predicted the allowed region of model parameter space in
the light of constraints arising from the measurements of both $N_{\rm eff}$
and dark matter relic density by Planck. It turns out that the impact of the
dark sector dynamics on $N_{\rm eff}$ is significant in case of a non-
hierarchical mass spectrum of the dark sector particles.

##  I Introduction

The increasing tension due to non-detection of dark matter at several direct
and indirect detection experiments [1, 2, 3, 4, 5] motivate to explore
alternative scenarios beyond weakly interacting massive particle (WIMP)
paradigm [6, 7, 8, 9, 10]. One of the attractive proposals is freeze in
mechanism [11, 12] where the DM abundance results from decays or annihilations
of visible sector particles. In this scenario, the DM never reaches thermal
equilibrium with the Standard Model (SM) bath due to its feeble interaction
strength. Such DM candidates are popularly dubbed as feebly interacting
massive particle (FIMP). Interestingly, the final abundance of the FIMP dark
matter is sensitive to the initial production mechanism in contrast to WIMP. A
different kind of scenario has been discussed in [13, 14, 15, 16] where the
dark matter is a part of a secluded dark sector and the annihilation of the
dark matter into the other dark sector particles set the relic abundance. This
is known as the secluded dark sector freeze-out.

The presence of any beyond Standard Model (BSM) light species can leave
detectable imprints at Planck satellite experiment. The Planck data provides
the number of relativistic degrees of freedom at the time of Cosmic Microwave
Background (CMB) to be $N_{\rm eff}^{\rm CMB}=2.99^{+0.34}_{-0.33}$ with
95$\%$ confidence limit (C.L) [17]. In SM alone, an accurate analysis of
neutrino decoupling by including exact collision terms, effect of neutrino
oscillations, and finite temperature QED estimates $N_{\rm eff}^{\rm SM}=3.04$
at the time of (CMB) [18, 19, 20, 21]. Earlier efforts of constraining light
BSM degrees of freedom can be found in [22, 23, 24, 25, 26, 27, 28, 29, 30].
These studies are based on the assumptions that SM neutrino decouples
instantaneously and the new BSM particle is in thermal equilibrium either with
the neutrino bath or with the photon bath. A more detailed analysis in [31,
32, 33, 34] by incorporating non-instantaneous SM neutrino decoupling imposes
lower bound on masses of BSM degrees of freedom from the measurement of
$N_{\rm eff}^{\rm CMB}$. In addition to that, during Big Bang Nucleosynthesis
(BBN), non-negligible abundances of the additional light species increase the
Hubble parameter which in turn may alter the abundances of Helium-4 and
Deuterium. The measured abundances of Helium-4 and Deuterium give rise to an
upper bound on number of relativistic degrees of freedom at the time of BBN
($N_{\rm eff}^{\rm BBN}$) to be $2.878\pm 0.278$ with $68.3\%$ C.L [35]. Thus
light species with mass $\lesssim\mathcal{O}(1)$MeV can also contribute to the
$N_{\rm eff}^{\rm BBN}$ [24, 36, 37, 38].

The measurement of effective relativistic degrees of freedom can be used as a
probe of light dark matter as well. A concise discussion on other possible
modes of probing a light dark sector can be found in [39, 40, 41]. In [42, 43,
44, 45, 31, 46, 47], it has been argued that a light dark matter can
contribute sizeably to $N_{\rm eff}^{\rm CMB}$ if it remains in thermal
contact with either neutrino or electron bath till late time. This poses a
lower bound on the MeV scale DM which differs depending on the spin of the
dark matter candidate.

In this work, we consider an MeV scale non-minimal dark sector consisting of
two real gauge singlet scalar fields. The lightest scalar is assumed to be
neutrinophilic. The heavier scalar, being stable in the Universe lifetime can
be identified as a viable dark matter candidate. Both the scalars interact
very feebly with the visible sector and hence never attain thermal equilibrium
with the SM bath, thus form a secluded dark sector. First, the non-thermal
production of lighter neutrinophilic scalar takes place from the inverse decay
of SM neutrinos and subsequently its annihilation yields the dark matter. The
production rate of DM depends on two factors. One is the abundance of the
mediator particle and secondly DM coupling strength to the mediator also
determines the efficiency of DM yield. After the production, the dark matter
annihilates back to the lighter scalar and in fact a stronger coupling leads
to formation of internal dark sector equilibrium with a different temperature
other than that of SM. Finally, when the ratio between dark sector interaction
strength and the expansion rate of the Universe comes down below unity, the
dark matter decouples and freezes out.

We anticipate that such a non-trivial dynamics of MeV scale dark sector can
have two fold impacts on the observation of Planck measurement of $N_{\rm
eff}^{\rm CMB}$. Firstly, if the lighter scalar decays to SM neutrinos during
or after neutrino decoupling, it raises the neutrino bath temperature and in
turn increases $N_{\rm eff}^{\rm CMB}$. Secondly, when mass hierarchy between
dark matter and the lighter scalar is small, it may increase the $N_{\rm
eff}^{\rm CMB}$ further. In order to address such possibilities, we have
performed a detailed study by solving system of coupled differential equations
for the evolutions of number densities and the temperatures for the relevant
species. We also include constraints on DM relic density in our present
analysis. The bound on the measurement of $\Delta N_{\rm eff}^{\rm CMB}$ are
expected to be more stringent with the next generation CMB experiments with
improved sensitivity. In fact the proposed CMB Stage IV (CMB-S4) experiment
[48] has the ability to strengthen the upper bound of $\Delta N_{\rm eff}^{\rm
CMB}$ upto 0.06 at 95% C.L and can potentially probe our proposed setup.

##  II Basic set up

In this section, we describe the particle contents of our proposed set up and
their interaction pattern. We consider an extension of SM with two real gauge
singlet scalar fields $S$ and $\phi$. The effective Lagrangian of our interest
is given by,

$\displaystyle-\mathcal{L}\supset\frac{1}{2}m_{S}^{2}S^{2}+\frac{1}{2}m_{\phi}^{2}\phi^{2}+\frac{\lambda}{4}S^{2}\phi^{2}+\dfrac{\kappa}{2}\phi\overline{\nu_{i}}\nu_{i},$
(1)

where we have considered the relevant mass scales ($m_{S}$ and $m_{\phi}$) and
the coupling coefficients ($\lambda$ and $\kappa$) as real and positive. In
our framework, the scalar field $\phi$ is purely neutrinophilic whereas $S$ is
the DM candidate. The stability of $S$ can be ensured by imposing a
$\mathbb{Z}_{2}$ symmetry under which $S$ has odd charge. The last term of Eq.
1 describes the interaction of the dark sector with the neutrino bath. Here
the couplings of $\phi$ with the neutrinos are flavor universal and $i$
represents the SM neutrino flavors. In our analysis, we work in the regime
$m_{\phi}<m_{S}$. In the opposite limit where the DM is lighter than the
mediator particle, the dark matter $S$ has negligible impact on $N_{\rm
eff}^{\rm CMB}$. This is because in this case we need smaller values of the
portal coupling to satisfy the relic density constraint, which in turn implies
that the impact on $N_{\rm eff}^{\rm CMB}$ is negligible.

The effective vertex $\phi\overline{\nu}\nu$ is clearly not invariant under
$SU(2)_{L}\otimes U(1)_{Y}$. This can be originated from an $SU(2)_{L}\otimes
U(1)_{Y}$ gauge invariant set up which is valid at higher energy scale. In
appendix A, we have discussed two such realisations, guided by the
introduction of a discrete symmetry $\mathbb{Z}_{4}$ in both cases, under
which $\phi$ transforms non-trivially. This choice automatically forbids the
interaction terms involving $\phi$ and other SM particles at linear order in
$\phi$. Additionally, there exist other renormalisable terms such as
$\phi^{2}|\Phi|^{2},\,S^{2}|\Phi|^{2}$ in the Lagrangian, where $\Phi$ is the
SM Higgs doublet. These may lead to new production channels for the dark
sector (DS) particles. However, in the present work we are interested in a
specific scenario where the DS gets populated dominantly from the SM neutrinos
(for example in [49]). This can be realized by either choosing the Higgs
portal coupling ($\phi^{2}|\Phi|^{2}$ etc.) parameters to be tiny enough or
assuming the reheating temperature of the Universe much lower than the Higgs
boson mass scale [50, 51], making the Higgs abundance Boltzmann suppressed in
the early Universe. In both the cases, the DS does not equilibrate with the SM
bath through Higgs portal.

Figure 1: Feynman diagrams indicating productions, annihilations of particles
as well as decay mode of $\phi$.

In the present study, we aim to probe the secluded non-minimal dark sector
(DS) with field contents $S$ and $\phi$ at Planck experiment. The dark sector
is secluded when $S$ and $\phi$ never reach thermal equilibrium with the SM
bath. The non-thermalisation of dark sector with SM bath implies,

$\displaystyle\sum_{i}n_{\nu_{i}}^{\rm eq}(T_{\nu_{i}})\langle\sigma
v\rangle_{\overline{\nu}_{i}\nu_{i}\rightarrow\phi}<H,$ (2)

where $\langle\sigma v\rangle_{\overline{\nu}_{i}\nu_{i}\rightarrow\phi}$ is
the thermally averaged cross section for the inverse decay of $\phi$ and $H$
represents the Hubble parameter of the Universe. A conservative check for Eq.
2 can be performed at the freeze in temperature $T_{\nu}\sim m_{\phi}$ of
$\phi$ yield. Then using the standard analytical expressions for
$\langle\sigma v\rangle_{\overline{\nu}_{i}\nu_{i}\rightarrow\phi}$,
$n_{\nu_{i}}^{\rm eq}$ and Hubble parameter, the Eq. 2 can be converted to,

$\displaystyle\left(\dfrac{m_{\phi}}{1\rm
MeV}\right)>0.3\times\left(\dfrac{y}{10^{-10}}\right)^{2}\,\,,$ (3)

where $y^{2}=\sum_{i}\kappa^{2}$. We have confirmed numerically that Eq. 3 is
a reasonable check to keep the DS and SM neutrino out of thermal equilibrium
for any temperature of the SM neutrino bath.

We start with the assumption that both $\phi$ and $S$ have zero initial
abundance. The inverse decay of SM neutrino yields $\phi$ and conversion of
$\phi$ takes place through $\phi\phi\rightarrow SS$ process. At a later stage
when number density of $S$ is sufficient, the reverse process
$SS\rightarrow\phi\phi$ turns effective (see Fig. 1). If the rate of such
forward and reverse interactions in dark sector are fast enough, a dark
equilibrium can be formed with a new temperature $T_{D}=T_{S}\neq
T_{\gamma},\,T_{\nu}$. Assuming the dark sector is internally thermalised, we
have estimated the dark sector temperature $T_{D}$ by solving the Boltzmann
equation of total dark sector energy density $\rho_{D}=\rho_{S}+\rho_{\phi}$.
The equilibration of the dark sector can be ensured by the following
condition.

$\displaystyle n_{S}^{\rm eq}(T_{D})\langle\sigma
v\rangle_{SS\rightarrow\phi\phi}>H,$ (4)

where we have considered $\langle\sigma v\rangle_{SS\rightarrow\phi\phi}$ to
be $s$-wave dominated. Simultaneous requirements of dark sector equilibrium
and non-thermalisation of dark and visible sectors put constraints on the
model parameters and the consequences of these constraints will be discussed
in section III.

We are particularly interested in the dynamics of a MeV scale secluded dark
sector in presence of dark $S-\phi$ equilibrium. The impact of such MeV scale
dark sector on Planck can be parameterized by the effective number of
relativistic species at CMB as defined by,

$\displaystyle N_{\rm eff}^{\rm
CMB}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\left(\frac{\rho_{\nu}}{\rho_{\gamma}}\right),$
(5)

where $\rho_{\gamma}=2\times\frac{\pi^{2}}{30}T_{\gamma}^{4}$,
$\rho_{\nu}=2\times 3\times\frac{7}{8}\times\frac{\pi^{2}}{30}T_{\nu}^{4}$ are
the energy densities of photon and neutrino baths respectively. Before
neutrino decoupling, the temperatures of the neutrino bath and the photon bath
are equal. In the post neutrino decoupling era, $T_{\gamma}$ and $T_{\nu}$
evolve separately. If all the BSM particles are having mass much larger than
eV scale, they do not contribute to the energy budget during CMB and hence in
Eq. 5, one can write,

$\displaystyle N_{\rm eff}^{\rm
CMB}=3\times\left(\frac{11}{4}\right)^{4/3}\left(\frac{T_{\nu}}{T_{\gamma}}\right)^{4}.$
(6)

Similarly one can define the number of relativistic degrees of freedom at the
time of BBN and it is given by ( considering all BSM particles are non-
relativistic during BBN)

$\displaystyle N_{\rm eff}^{\rm
BBN}=\dfrac{8}{7}\left(\dfrac{\rho_{\nu}}{\rho_{\gamma}}\right).$ (7)

##  III Boltzmann equations and Numerical analysis

To study the effect of the dark sector dynamics on $\Delta N_{\rm eff}$, we
need to solve the Boltzmann equations for comoving number densities $Y_{i}$
($i=\phi$, $S$), $\xi=T_{\nu}/T_{\gamma}$, and $\xi_{D}=T_{D}/T_{\gamma}$. In
deriving the above Boltzmann equations, we have used the following
assumptions.

  * •

The temperatures of all the neutrino flavors are same i.e.
$T_{\nu_{e}}=T_{\nu_{\mu}}=T_{\nu_{\tau}}\equiv T_{\nu}$.

  * •

The entropy of the photon bath is conserved since at late time $\phi$ can only
decay into a pair of neutrinos. Using the conservation of entropy in the
photon bath, we have defined $Y_{i}=n_{i}/s_{\gamma}$ and the relic abundance
of the DM is calculated accordingly111 We estimate the relic abundance of DM
using $\Omega_{S}h^{2}=\frac{g_{*s,\gamma}}{g_{*s}^{\rm today}}\times
2.745\times 10^{8}\times
Y_{S}^{(x\rightarrow\infty)}\times\left(\frac{m_{S}}{1\,\text{GeV}}\right)$
where $g_{*s,\gamma}=2$ and $g_{*s}^{\rm today}=3.91$..

  * •

The interaction rate between the dark sector particles such as $S$ and $\phi$
is sufficient enough to keep them in thermal equilibrium and they share a
common temperature $T_{D}$ (see discussion above Eq. 4).

Using the above assumptions, first we write the evolution of $\xi$ as a
function of $x\equiv m_{S}/T_{\gamma}$.

$\displaystyle\dfrac{d\xi}{dx}$
$\displaystyle=\dfrac{1}{x}\left[\xi-\dfrac{\mathcal{B}_{T_{\gamma}}(T_{\gamma},T_{\nu})}{\mathcal{B}_{T_{\nu}}(T_{\gamma},T_{\nu})}\right]\,\,\,,$
(8)

where

$\displaystyle\mathcal{B}_{T_{\nu}}(T_{\gamma},T_{\nu})=\dfrac{\mathcal{C}_{\rm
el}(T_{\gamma},T_{\nu})+\mathcal{C}_{\rm
inel}(T_{\gamma},T_{\nu})-\mathcal{C}_{\nu\nu\leftrightarrow\phi}-4H\rho_{\nu}}{\frac{\partial\rho_{\nu}}{\partial
T_{\nu}}}\,\,,$
$\displaystyle\mathcal{B}_{T_{\gamma}}(T_{\gamma},T_{\nu})=\dfrac{-\mathcal{C}_{\rm
el}(T_{\gamma},T_{\nu})-\mathcal{C}_{\rm
inel}(T_{\gamma},T_{\nu})-4H\rho_{\gamma}-3H(\rho_{e}+p_{e})}{\frac{\partial(\rho_{\gamma}+\rho_{e})}{\partial
T_{\gamma}}}\,\,.$

Here, $\mathcal{C}_{\rm el}(T_{\gamma},T_{\nu})$ and $\mathcal{C}_{\rm
inel}(T_{\gamma},T_{\nu})$ are the relevant collision terms for the elastic
and inelastic processes such as
$e^{+}e^{-}\leftrightarrow\bar{\nu}_{i}\nu_{i}$,
$e^{-}(e^{+})\nu_{i}(\bar{\nu}_{i})\leftrightarrow
e^{-}(e^{+})\nu_{i}(\bar{\nu}_{i})$,
$e^{+}(e^{-})\nu_{i}(\bar{\nu}_{i})\leftrightarrow
e^{+}(e^{-})\nu_{i}(\bar{\nu}_{i})$. $\mathcal{C}_{\nu\nu\leftrightarrow\phi}$
is the collision term for $\phi\leftrightarrow\bar{\nu}_{i}\nu_{i}$ where
$i=e,\mu,\tau$. The matrix amplitude square of these processes and the generic
form of collision terms for elastic, inelastic, and decay processes are given
in decay processes are given in table 3 of appendix B and appendix C
respectively. The energy densities of electron, photon, and neutrino are
respectively denoted by $\rho_{e}$, $\rho_{\gamma}$, and $\rho_{\nu}$. $p_{e}$
denotes the pressure of electron.

The evolution of the co-moving number densities of $\phi$ ($Y_{\phi}$), and
$S$ ($Y_{S}$) are governed by the following two coupled Boltzmann equations.

$\displaystyle\dfrac{dY_{\phi}}{dx}$ $\displaystyle=\dfrac{h_{\rm
eff}s_{\gamma}}{xH}\langle\sigma v_{\rm
rel}\rangle_{SS\to\phi\phi}^{T_{D}}\left[Y_{S}^{2}-\left(\dfrac{Y_{S}^{\rm
eq}(T_{D})}{Y_{\phi}^{\rm eq}(T_{D})}\right)^{2}Y_{\phi}^{2}\right]$
$\displaystyle+\dfrac{h_{\rm
eff}}{xH}\left(\langle\Gamma_{\phi}\rangle_{T_{\nu}}Y_{\phi}^{\rm
eq}(T_{\nu})-\langle\Gamma_{\phi}\rangle_{T_{D}}Y_{\phi}\right)\,\,\,$ (10)

$\displaystyle\dfrac{dY_{S}}{dx}$ $\displaystyle=-\dfrac{h_{\rm
eff}s_{\gamma}}{xH}\langle\sigma v_{\rm
rel}\rangle_{SS\to\phi\phi}^{T_{D}}\left[Y_{S}^{2}-\left(\dfrac{Y_{S}^{\rm
eq}(T_{D})}{Y_{\phi}^{\rm eq}(T_{D})}\right)^{2}Y_{\phi}^{2}\right]\,\,\,.$
(11)

Here, $s_{\gamma}$ is the entropy density of the photon bath, and $Y^{\rm
eq}_{\phi(S)}(T_{D})$ denotes the equilibrium co-moving number density of
$\phi(S)$ at hidden sector temperature $T_{D}$. The Hubble parameter $H$ and
$h_{\rm eff}$ are defined as follows.

$\displaystyle H$ $\displaystyle=\sqrt{\dfrac{8\pi}{3M^{2}_{\rm
Pl}}\left(\rho_{\gamma}+\rho_{\nu}+\rho_{S}+\rho_{\phi}\right)}\,\,,$
$\displaystyle h_{\rm eff}$ $\displaystyle=1+\dfrac{1}{3}\dfrac{d\ln
g_{*s,\gamma}(T_{\gamma})}{d\ln T_{\gamma}}\,\,,$ (12)

where $M_{\rm Pl}=1.22\times 10^{19}\,\rm GeV$ and $g_{*s,\gamma}(T_{\gamma})$
is the relativistic degrees of freedom contributing to the entropy density of
the photon bath. The energy densities of photon and neutrino baths are
indicated by $\rho_{\gamma}$, $\rho_{\nu}$ respectively. Note that, in the
expression of the Hubble parameter, we have included the energy densities of
dark sector particles in addition to the contributions from SM fields. On the
right hand side (RHS) of Eq. 10 and Eq. 11, $\langle\sigma v_{\rm
rel}\rangle^{T_{D}}_{SS\to\phi\phi}$ stands for the thermally averaged
annihilation cross section of $SS\to\phi\phi$ at temperature $T_{D}$. The
analytical form of $\sigma_{SS\to\phi\phi}$ is given by,

$\displaystyle\sigma_{SS\to\phi\phi}=\dfrac{\lambda^{2}}{32\pi
s}\sqrt{\dfrac{s-4m_{\phi}^{2}}{s-4m_{S}^{2}}}\,\,,$ (13)

where $\sqrt{s}$ is the centre of mass energy of the collision.
$\langle\Gamma_{\phi}\rangle_{T_{D}}$ and
$\langle\Gamma_{\phi}\rangle_{T_{\nu}}$ on the RHS of Eq. 10 are the thermally
averaged total decay width of $\phi$ at temperature $T_{D}$ and $T_{\nu}$
respectively. A general expression of the thermally averaged decay width at
temperature $T_{i}$ is given by
$\langle\Gamma_{\phi}\rangle_{T_{i}}=\Gamma_{\phi}\dfrac{K_{1}(m_{\phi}/T_{i})}{K_{2}(m_{\phi}/T_{i})}$,
where $\Gamma_{\phi}=\frac{y^{2}m_{\phi}}{16\pi}$ is the total decay width of
$\phi$. In deriving the Boltzmann equations for $Y_{\phi}$ and $Y_{\chi}$, we
have also taken into account the effect of non-zero chemical potential for
dark sector particles as evident from the presence of square of the ratio
$(Y_{S}^{\rm eq}(T_{D})/Y_{\phi}^{\rm eq}(T_{D}))$ in the RHS of Eq. 10.

As we discussed earlier, the dark sector is connected to SM by the portal
coupling $y$ and our choice of $y$ keeps dark and the visible sector out of
thermal equilibrium. Initially inverse decay of $\phi$ populates the dark
sector and subsequently $\phi\phi\to SS$ process occurs. The Boltzmann
equation that governs the evolution of total dark sector energy density
($\rho_{D}=\rho_{\phi}+\rho_{S}$) is read as

$\displaystyle\frac{d\rho_{D}}{dt}+3H(\rho_{D}+p_{D})=\mathcal{C}_{\nu\nu\leftrightarrow\phi},$
(14)

where $p_{D}$ stands for the pressure and
$\mathcal{C}_{\nu\nu\leftrightarrow\phi}$ includes the collision integrals for
both decay and inverse decay processes of $\phi$. We use the assumption of
Maxwell- Boltzmann statistics and write $\rho_{D}$ and $p_{D}$ as function of
$T_{D}$ as given by

$\displaystyle\rho_{D}=\left(\frac{Y_{\phi}}{Y_{\phi}^{\rm
eq}(T_{D})}\right)\rho_{\phi}^{\rm eq}(T_{D})+\left(\frac{Y_{S}}{Y_{S}^{\rm
eq}(T_{D})}\right)\rho_{S}^{\rm eq}(T_{D}),$ (15) $\displaystyle
p_{D}=\left(\frac{Y_{\phi}}{Y_{\phi}^{\rm eq}(T_{D})}\right)p_{\phi}^{\rm
eq}(T_{D})+\left(\frac{Y_{S}}{Y_{S}^{\rm eq}(T_{D})}\right)p_{S}^{\rm
eq}(T_{D}),$ (16)

In the RHS of the above equation, the quantities, $\rho_{\phi}^{\rm eq}$ and
$p_{\phi}^{\rm eq}$ represent the energy density and pressure at equilibrium
as given by

$\displaystyle\rho_{i}^{\rm
eq}(T_{D})=\frac{g_{i}}{2\pi^{2}}T_{D}^{4}\left[3\left(\frac{m_{i}}{T_{D}}\right)^{2}K_{2}\left(\frac{m_{i}}{T_{D}}\right)+\left(\frac{m_{i}}{T_{D}}\right)^{3}K_{1}\left(\frac{m_{i}}{T_{D}}\right)\right],$
$\displaystyle p_{i}^{\rm
eq}=\frac{g_{i}}{2\pi^{2}}T_{D}^{4}\left(\frac{m_{i}}{T_{D}}\right)^{2}K_{2}\left(\frac{m_{i}}{T_{D}}\right),$
(17)

where once again, the ratio $\left(\frac{Y_{i}}{Y_{i}^{\rm eq}(T_{D})}\right)$
in Eqs.15-16 takes care of the non-zero chemical potential during the
temperature evolution for the ith species. Utilizing these, we obtain,

$\displaystyle\frac{dT_{D}}{dx}=\frac{s_{\gamma}Y_{\phi}^{\rm
eq}(T_{D})Y_{S}^{\rm
eq}(T_{D})\left\\{\frac{1}{s_{\gamma}}\mathcal{C}_{\nu\nu\leftrightarrow\phi}-3H(Y_{\phi}+Y_{S})T_{D}\right\\}+F_{1}+F_{2}}{xH\sum_{i\neq
j}Y_{i}Y_{j}^{\rm eq}\left\\{\frac{d\rho_{i}^{\rm
eq}}{dT_{D}}-\frac{1}{n_{i}^{\rm eq}(T_{D})}\frac{dn_{i}^{\rm
eq}(T_{D})}{dT_{D}}\rho_{i}^{\rm eq}(T_{D})\right\\}},$ (18)

where,

$\displaystyle F_{1}=s_{\gamma}\langle\sigma
v\rangle_{SS\to\phi\phi}^{T_{D}}\left[Y_{S}^{2}-\left(\frac{Y_{S}^{\rm
eq}(T_{D})}{Y_{\phi}^{\rm
eq}(T_{D})}\right)^{2}Y_{\phi}^{2}\right]\left(\rho_{S}^{\rm eq}Y_{\phi}^{\rm
eq}-\rho_{\phi}^{\rm eq}Y_{S}^{\rm eq}\right),$ $\displaystyle
F_{2}=-\rho_{\phi}^{\rm eq}(T_{D})Y_{\chi}^{\rm
eq}(T_{D})\left[\langle\Gamma_{\phi\to\nu\nu}\rangle_{T_{\nu}}Y_{\phi}^{\rm
eq}(T_{\nu})-\langle\Gamma_{\phi\to\nu\nu}\rangle_{T_{D}}Y_{\phi}\right].$

and the indices $i$ and $j$ stand for $\phi$ and $S$.

(a) (b) (c)

Figure 2: (a) Variation of ${T_{D}}$ as a function of $T_{\gamma}$ for
$m_{S}=5.1\rm MeV$, $\lambda=6.61\times 10^{-6}$ (black line) and $m_{S}=12\rm
MeV$, $\lambda=5.18\times 10^{-6}$ (blue line) and these parameters yield the
correct DM relic abundance [17]. (b) Evolution of co-moving number density of
$\phi$ ($Y_{\phi}$) (dashed lines) and $S$ ($Y_{S}$) (solid lines) as a
function of the SM bath temperature $T_{\gamma}$ for the same choice of
parameters as in 2(a). (c) The variation of $\Delta N^{\rm CMB}_{\rm eff}$ as
a function of $T_{\gamma}$. Here the parameter choices and color codes are
same as Fig. 2(a). The red horizontal solid line indicates the upper limit on
$\Delta N_{\rm eff}$ at CMB by Planck $2\sigma$ data. In all plots, we have
considered $m_{\phi}=5\rm MeV$ and $y=2\times 10^{-10}$.

(a) (b) (c)

Figure 3: $\Delta N^{\rm CMB}_{\rm eff}$ as a function of $\lambda$ for (a)
$m_{\phi}=3$MeV, (b) $m_{\phi}=7$MeV and (c) $m_{\phi}=10$MeV. The green
points are for $\delta=(m_{S}-m_{\phi})=0.5$MeV whereas blue points indicate
the $\Delta N_{\rm eff}^{\rm CMB}$ values for $\delta=2$MeV. Here we have
considered $y=2\times 10^{-10}$. The symbol ‘$\star$’ stands for the $\lambda$
value that yields correct relic abundance [17]. In all three figures, the grey
region is disallowed from the measurement of $\Delta N^{\rm CMB}_{\rm eff}$ by
Planck within $2\sigma$ limit.

At first, we would like to calculate the value of $\Delta N_{\rm eff}^{\rm
CMB}$ in SM. In order to do so, we have numerically solved Eq. 8 in the
absence of new physics. Following [52, 53, 54], we have computed the collision
terms for SM $\nu$ decoupling, considering Fermi-Dirac (FD) distribution
function for all the SM fermions and appropriate Pauli blocking factors. For
simplification, we have neglected the mass of the electron ($m_{e}$) since at
the time of $\nu$ decoupling, $T_{\gamma}\gg m_{e}$. We have obtained the
value of $T_{\nu}/T_{\gamma}=0.714$ at CMB and the correspondingly $N_{\rm
eff}^{\rm CMB,SM}\simeq 3.02$. Now we proceed to estimate $N_{\rm eff}^{\rm
CMB}$ in presence of the MeV scale secluded dark sector. We solve the coupled
system of Boltzmann equations (Eqs.8-11,18) with the assumption of initial
zero abundance for the dark sector particles and $T_{D}^{\rm initial}=0$ to
study the dark sector dynamics along with the estimate of $\Delta{\rm N}_{\rm
eff}^{\rm CMB}$. The obtained numerical results in presence of the new physics
are presented as follows.

Figure 4: The variation of $\Delta N^{\rm CMB}_{\rm eff}$ as a function of
$T_{\gamma}$ for $y=2\times 10^{-10}$ (solid line), $y=1.5\times 10^{-10}$
(dashed line), and $y=1.25\times 10^{-10}$ (dotted line). Here we choose
$m_{S}=5.1\rm MeV$, $m_{\phi}=5\rm MeV$. The red horizontal solid line
indicates the upper limit on $\Delta N_{\rm eff}$ at CMB by Planck $2\sigma$
data.

In Fig. 2, we demonstrate the evolution patterns of $T_{D}$, $Y_{\phi}$,
$Y_{S}$, and $\Delta N_{\rm eff}^{\rm CMB}$ as functions of the SM temperature
$T_{\gamma}$. We have fixed $m_{\phi}$ and $y$ at $5$MeV and $2\times
10^{-10}$ respectively and considered two specific DM masses which are 5.1MeV
and 12MeV. The corresponding $\lambda$ value for each of the DM mass is fixed
by the relic density requirement ($\Omega_{S}h^{2}\sim 0.12$). The dark sector
temperature starts from zero initial value and increases at the early stage of
its evolution due to the production of $\phi$ from
$\bar{\nu}_{i}\nu_{i}\to\phi$ process. The temperature $T_{D}$ continues to
grow until it reaches a maximum value, after which it starts to redshift like
radiation upto $T_{D}\lesssim m_{\phi}$ due to the expansion of the Universe.
Around $T_{\gamma}\sim m_{\phi}$, we notice that the redshift for $T_{D}$ gets
slower. This occurs since the production rate of $\nu\nu\rightarrow\phi$
reaches maximum around this temperature. For $m_{S}=5.1$MeV we also observe
the impact of dark matter freeze out in the evolution of $T_{D}$, where late
time conversion process $SS\to\phi\phi$ slows down the redshift of $T_{D}$ for
the second time. This pattern is not visible for $m_{S}=12$MeV as in this case
the dark matter freeze out occurs at a higher temperature when the
$\nu\nu\to\phi$ process is still active. Finally after the freeze out of the
dark matter and decay of $\phi$, the dark sector temperature redishifts like a
non-relativistic matter field with usual scale factor ($a$) dependence
$T_{D}\propto a^{-2}$.

In Fig. 2(b), we show the evolutions of $Y_{S}$ (solid lines) and $Y_{\phi}$
(dashed lines) as functions of $T_{\gamma}$ considering two different choices
of $(m_{S},\lambda)$ as indicated by black and blue curves for $m_{S}=5.1$MeV
and 12MeV, respectively. The corresponding values of $\lambda$ are chosen
accordingly to yield correct relic abundance with a larger DM mass requiring
lower value of $\lambda$ to obey the relic density bound. On the other hand,
in Fig. 2(c), we show the variation of $\Delta N_{\rm eff}^{\rm CMB}$ as a
function of $T_{\gamma}$ for the same choice of parameters as considered in
Fig. 2(b). It is observed from Fig. 2 that the dark matter $S$ of mass $12\rm
MeV$ freezes-out much earlier and at the time of DM freeze-out, the production
of $\phi$ is still in progress. As a result, when $\phi$ starts decaying to a
pair of SM neutrinos, the DM has already decoupled from $\phi$ bath and the
final value of $\Delta N_{\rm eff}$ is governed mostly by $m_{\phi}$ and $y$.
However the situation turns a bit different as we decrease the mass hierarchy
of $S$ and $\phi$. In case of $m_{S}=5.1\rm MeV$, decoupling of $S$ is delayed
in comparison to the earlier scenario and at the time of decay of $\phi$,
$SS\to\phi\phi$ is still active with $S$ having equivalent comoving abundance
of $\phi$ around $T_{D}\sim m_{\phi}$. As a result, the depletion of $\phi$
has been counterbalanced by the production $SS\to\phi\phi$, which in turn
enhances $\Delta N_{\rm eff}^{\rm CMB}$. Thus, $\Delta N_{\rm eff}^{\rm CMB}$
increases as we reduce the mass gap between $S$ and $\phi$.

Next, the dependence of the parameter $\lambda$ on $\Delta N_{\rm eff}^{\rm
CMB}$ is shown in Fig. 3 for $y=2\times 10^{-10}$ and three different values
of $m_{\phi}$. The lowest chosen value of $m_{\phi}$ is 3MeV in Fig. 3 and for
this choice, $\phi$ remains out of equilibrium from SM neutrino bath as can be
confirmed from Eq. 3. For larger $m_{\phi}$ with same $y$, the Eq. 3 is easier
to satisfy. We consider two different hierarchical patterns of dark sector
particles which are $(m_{S}-m_{\phi})$= 0.5MeV and 2MeV. Let us define
$\delta=(m_{S}-m_{\phi})$ for convenience. The dark sector reaches internal
equilibrium for the chosen ranges of $\lambda$ in Figs. 3(a)-3(c) and it has
been numerically checked. We find that $\Delta N_{\rm eff}^{\rm CMB}$ is
insensitive to $\lambda$. This is because when the dark sector is internally
thermalized, the abundance of $\phi$ at the onset of its decay is mostly
governed by $y,m_{\phi}$ and $m_{S}$. We also notice $\Delta N_{\rm eff}$
increases with the decrease in mass gap between $S$ and $\phi$ as earlier
observed in Fig. 2(c). In all three subfigures of Fig. 3, the point
highlighted by ‘$\star$’ indicates the required $\lambda$ to yield correct
relic density for DM [17]. Another important outcome of Fig. 3 is related to
the dependence of $m_{\phi}$ on $\Delta N_{\rm eff}^{\rm CMB}$. When
$m_{\phi}$ is fixed at 3MeV, the predictions for $\Delta N_{\rm eff}^{\rm
CMB}$ considering both the mass differences $\delta=$ 0.5MeV and 2MeV are
below the Planck limit. With the increase of $m_{\phi}\sim 7$MeV with
$\delta=0.5$MeV, the prediction for $\Delta N_{\rm eff}^{\rm CMB}$ enters into
the disallowed region. This happens since larger $m_{\phi}$ increase the decay
width of $\phi$ and causes enhanced rate of entropy injection into the SM
neutrino bath. Further increase of $m_{\phi}$ to 10MeV shows different pattern
and $\delta=0.5$MeV again makes a comeback to the allowed region of $\Delta
N_{\rm eff}^{\rm CMB}$. This is due to the fact that a sufficient large
$m_{\phi}$ would make $\phi$ to decay early with less impact on evolution of
$T_{\nu}$. On the other hand for $\delta=2$MeV, the predictions for $\Delta
N_{\rm eff}^{\rm CMB}$ always remain inside the Planck favored region.

As we find from Fig. 2(c), the predictions on $\Delta N_{\rm eff}^{\rm CMB}$
for $m_{\phi}=5$MeV and $m_{S}=5.1$MeV are disfavored when $y=2\times
10^{-10}$. Now, let us examine the effects of varying $y$ for the same set of
$m_{\phi}$ and $m_{S}$. In Fig. 4, we show the variation of $\Delta N^{\rm
CMB}_{\rm eff}$ as a function of $T_{\gamma}$ for three different choices of
the portal coupling $y$ by fixing $m_{\phi}$ and $m_{S}$ at 5MeV and 5.1MeV
respectively. The choices for $\lambda$ are made to obtain the best fit value
of observed relic abundance $\Omega h^{2}=0.1200\pm 0.0012$ [17] after
formation of internal dark equilibrium. Earlier we found that $y=2\times
10^{-10}$ is ruled out by present CMB data which is observed once again in
Fig. 4. However reducing $y$ suppresses the energy injection rate and thus
decreases $\Delta N^{\rm CMB}_{\rm eff}$ and hence would be allowed by the
Planck data.

In Fig. 5, we have shown the evolutions of $N_{\rm eff}^{\rm BBN}$ following
the definition of Eq. 6 considering three different values of $m_{\phi}$ with
the $m_{S}$ is fixed by $m_{S}=m_{\phi}+0.5$MeV and $y=2\times 10^{-10}$. The
magnitude of $\lambda$ is determined by the requirement of obeying the relic
density bound. It turns out that $N_{\rm eff}^{\rm BBN}$ at $T_{\rm BBN}\sim
1$MeV increases with $m_{\phi}$ however stays inside the allowed $2\sigma$
limit [35]. Importantly, for lighter $m_{\phi}$, $N_{\rm eff}^{\rm BBN}$ is
smaller than its SM value which is approximately three. This owes to the fact
that for the lightest $m_{\phi}=3$MeV, as considered here, the production of
$\phi$ continues even during BBN from the inverse decay of $\phi$. This
reduces the temperature of the SM neutrino bath at BBN and we observe a dip in
the value of $N_{\rm eff}^{\rm BBN}$ at $T_{\rm BBN}\sim 1$MeV for
$m_{\phi}=3$MeV. Overall, we find that BBN requirements do not impose any
serious constraints to our model parameter space. A detailed analysis in order
to study the impact of a secluded non-minimal dark sector on the abundances of
the primordial light elements (e.g. Helium-4 and Deuterium [55, 56]) is beyond
the scope of the present article and will be pursued elsewhere.

Finally in Fig. 6, we show the predictions for $\Delta N_{\rm eff}^{\rm CMB}$
in $m_{S}-\lambda$ plane. We have fixed $m_{\phi}$ and $y$ at 5MeV and
$2\times 10^{-10}$ respectively. These particular choices ensure that $\phi$
does not equilibrate with SM neutrinos as followed from Eq. 3. The
overabundant region of DM is highlighted with light grey color. We have
numerically found that the condition for dark sector equilibrium following Eq.
4 with $m_{S}\lesssim 20$MeV poses much weaker constraint compared to the one
from observed DM relic in the $m_{S}-\lambda$ plane. Hence we do not show it
in Fig. 6.

Figure 5:  The evolution of $N_{\rm eff}^{\rm BBN}$ with $T_{\gamma}$ is shown
for three different sets of $(m_{\phi},m_{s})$ with $y=2\times 10^{-10}$. The
corresponding $\lambda$ value for each set is fixed by the relic density
requirement. The magenta and blue horizontal solid lines denote the upper and
lower limit of $\Delta N_{\rm eff}$ at the time of BBN at $2\sigma$ [35]
respectively. The perpendicular dot-dashed line represents $T_{\rm BBN}=1$MeV.

The $\Delta N_{\rm eff}^{\rm CMB}$ lines for different $m_{S}$ are also
depicted in Fig. 6. As explained earlier the predictions for $\Delta N_{\rm
eff}^{\rm CMB}$ remain insensitive to variation in $\lambda$. Clearly, one can
see that the value of $\Delta N_{\rm eff}^{\rm CMB}$ is maximum when $m_{S}$
turns closer to $m_{\phi}$. We also observed that as $m_{S}$ turns larger, its
impact on $\Delta N_{\rm eff}^{\rm CMB}$ gets reduced. The latest Planck
2$\sigma$ bound on $\Delta N_{\rm eff}^{\rm CMB}$ restricts $m_{S}$ to be
larger than around 7 MeV for $m_{\phi}=5$MeV and $y=2\times 10^{-10}$.

Figure 6:  Allowed parameter space in $m_{S}-\lambda$ plane from the relic
density bound along with the $\Delta N^{\rm CMB}_{\rm eff}$ predictions for
$m_{\phi}=5\rm MeV$ and $y=2\times 10^{-10}$. The overabundant region of DM is
shown by grey shade. In the top x-axis, we show the predictions for $\Delta
N_{\rm eff}^{\rm CMB}$ as function of DM mass $m_{S}$.

##  IV Summary and Conclusion

In this work, we examine the scope of probing a MeV scale secluded dark sector
from the observation of $\Delta N_{\rm eff}^{\rm CMB}$ by Planck. For the
purpose, we have considered a simple structure of secluded dark sector
comprising of two SM gauge singlet scalars $S$ and $\phi$. The dark sector
never reaches thermal equilibrium with the SM sector. The lighter scalar
$\phi$ is neutrinophilic and we identify the heavier one $S$ as our DM
candidate. Initially the dark sector particle $\phi$ gets populated from the
SM sector by the $\overline{\nu}_{i}\nu_{i}\to\phi$ process and afterwards
$\phi\phi\rightarrow SS$ process produces $S$. Depending on the interaction
strength between $S$ and $\phi$, dark thermalisation can occur which causes
freeze-out of $S$ via $SS\to\phi\phi$ process. We constrain our parameter
space from the observations of both dark matter relic density and measurement
of $N_{\rm eff}^{\rm CMB}$ by Planck ensuring internal dark thermal
equilibrium. We particularly emphasize the impact of the dark sector
parameters on the prediction of $N_{\rm eff}^{\rm CMB}$. Below we summarize
few important observations that came out from our analysis.

$\bullet$ We obtain unique predictions for the $\Delta N_{\rm eff}^{\rm CMB}$
for a particular DM mass that satisfies the relic density bound with suitable
value of $\lambda$, provided mass of the other scalar and its coupling
strength with the SM neutrinos are fixed.

$\bullet$ The dark matter mass has a nontrivial role on the prediction of
$\Delta N_{\rm eff}^{\rm CMB}$. The predictions for $\Delta{N_{\rm eff}^{\rm
CMB}}$ turn out to be maximum in the nearly degenerate spectrum of dark sector
particles and consideration of an increased hierarchy between the dark sector
particles reduces the impact of dark matter mass on $\Delta N_{\rm eff}^{\rm
CMB}$. We also find that prediction for $\Delta N_{\rm eff}^{\rm CMB}$ is
insensitive to the dark sector interaction rate, provided the dark thermal
equilibrium is reached.

$\bullet$ Increasing $m_{\phi}$ upto a certain value keeping other parameters
($y,~{}\lambda$ and $\delta$) fixed, enhances $\Delta{N_{\rm eff}^{\rm CMB}}$
hence restricted by the present Planck limit. When $m_{\phi}$ is too large,
the effect on $\Delta{N_{\rm eff}^{\rm CMB}}$ of $\phi$ diminishes.

$\bullet$ For a particular $m_{\phi}$ with a constant $y$, we are able to
impose a lower bound on the DM mass from the Planck $2\sigma$ data on $\Delta
N_{\rm eff}^{\rm CMB}$. As an example we find $m_{S}\gtrsim 7$ MeV when
$m_{\phi}=5$MeV and $y=2\times 10^{-10}$. Additionally, a lower bound on
$\lambda$ can be derived such that the DM relic abundance remain below 0.12 as
experimentally favored.

In summary, an MeV scale secluded non-minimal dark sector is difficult to
probe at dedicated dark matter search experiments. In view of this, we propose
that the measurement of $\Delta N_{\rm eff}^{\rm CMB}$ can be used as a tool
to test such a decoupled dark sector. Our work indeed shows that some of the
parameter space yield observable $\Delta N_{\rm eff}^{\rm CMB}$ and are
already within the reach of the Planck sensitivity. Upcoming CMB stage IV
experiments with improved sensitivity should be able to probe/refute the
allowed region of our model parameter space further.

##  V Acknowledgements

SG would like to thank Anirban Biswas for some insightful discussions. SG
would also like to thank the University Grants Commission (UGC), Government of
India, for providing financial support as a senior research fellow. AKS is
supported by NPDF grant PDF/2020/000797 from Science and Engineering Research
Board (SERB), Government of India.

##  Appendix A Two examples of $SU(2)_{L}\otimes U(1)_{Y}$ gauge invariant
model

In this appendix, we discuss two $SU(2)_{L}\otimes U(1)_{Y}$ invariant
frameworks that give rise to $\phi\overline{\nu}\nu$ vertex effectively at MeV
scale.

Case I: We consider a variant of type-I seesaw framework with three gauge
singlet right handed (RH) neutrinos. We impose $\mathbb{Z}_{4}$ discrete
symmetry to write the desired Lagrangian. The charge assignments of the
relevant fields under the new $\mathbb{Z}_{4}$ are presented in table 1. We
propose the following Lagrangian for generating the $\phi\overline{\nu}\nu$
vertex.

$\displaystyle\mathcal{L}\supset-\dfrac{c_{\phi
1}}{2}\phi\,\overline{N_{R_{1}}^{C}}N_{R_{1}}-\dfrac{c_{\phi
2}}{2}\phi\,\overline{N_{R_{2}}^{C}}N_{R_{2}}-Y_{e1}\overline{L}_{e}\widetilde{\Phi}N_{R_{1}}-Y_{\mu
1}\overline{L}_{\mu}\widetilde{\Phi}N_{R_{2}}-Y_{\tau
3}\overline{L}_{\tau}\widetilde{\Phi}N_{R_{3}}-M_{12}\overline{N_{R_{1}}^{C}}N_{R_{2}}-\frac{M_{33}}{2}\overline{N_{R_{3}}^{C}}{N_{R_{3}}}+\text{
h.c. },$ (19)

where $c_{\phi 1}$, $c_{\phi 2}$ are dimensionless $\mathcal{O}(1)$ coupling
coefficients and $\widetilde{\Phi}=i\sigma_{2}\Phi^{*}$.

Fields | $\mathbb{Z}_{4}$
---|---
$N_{R_{1}}$ | $i$
$N_{R_{2}}$ | $-i$
$N_{R_{3}}$ | $1$
$\phi$ | $-1$
$\\{L_{e},L_{\mu},L_{\tau}\\}$ | $\\{i,-i,1\\}$
$\Phi$ | 1
Table 1: Charge assignment for the relevant fields as proposed in case I. The
dark matter field $S$ transforms trivially under the imposed $\mathbb{Z}_{4}$.

After electroweak symmetry breaking, the SM neutrino mass matrix is given by,

$\displaystyle m_{\nu}\simeq-m_{D}M_{R}^{-1}m_{D}^{T},$ (20)

where

$\displaystyle m_{D}=\frac{v}{\sqrt{2}}\begin{pmatrix}Y_{e1}&0&0\\\ 0&Y_{\mu
2}&0\\\ 0&0&Y_{\tau
3}\end{pmatrix}\text{~{}~{}~{}and~{}~{}~{}}M_{R}=\begin{pmatrix}0&M_{12}&0\\\
M_{12}&0&0\\\ 0&0&M_{33}\end{pmatrix},$ (21)

The active-sterile neutrino mixing is $\theta_{\alpha i}$, quantified as
$\theta_{\alpha i}\simeq\left(m_{D}M_{R}^{-1}\right)_{\alpha i}$ [57]. Due to
this active-sterile mixing, the first two terms in Eq. 19 can be translated to
$\phi\overline{\nu}\nu$ vertex with the coupling coefficient being
proportional to $\sim(\theta_{\mu 1}^{2}+\theta_{e2}^{2})$. This set up is
able to provide correct order of active neutrino mass ($\sim 0.1$eV) for
$M_{33},M_{12}\sim\mathcal{O}(1)$GeV and $Y\sim\mathcal{O}(10^{-7})$ with the
desired value of $\phi\overline{\nu}\nu$ coupling coefficient
$\sim\mathcal{O}(10^{-10})$. A detailed study on neutrino masses and mixing is
beyond the scope of the present paper. To prevent the production of DS
particles from tree level RH neutrino decay ($N\to\phi\nu$), we have assumed
that the reheating temperature ($T_{\rm RH}$) of the Universe is less than $1$
GeV which is consistent with the lower limit on the $T_{\rm RH}$ ($\gtrsim
4$MeV) from big bang nucleosynthesis [50].

Case II: We may also start with the following $SU(2)_{L}\otimes U(1)_{Y}$
invariant Lagrangian at dimension six level to effectively generate the
$\phi\overline{\nu}\nu$ vertex [58],

$\displaystyle\mathcal{L}\supset-\frac{c_{2}}{\Lambda^{2}}\phi\,\overline{\widetilde{L_{l}}}\,\Phi\,\widetilde{\Phi}^{\dagger}\,L_{l}+\text{~{}h.c.~{}},$
(22)

where $c_{2}$ is the dimensionless $\mathcal{O}(1)$ coupling coefficient,
$\Lambda$ is the cut off scale and $\tilde{L}_{l}=i\sigma_{2}L_{l}^{C}$ with
$l=\\{e,\mu.\tau\\}$.

Fields | $\mathbb{Z}_{4}$
---|---
$\phi$ | $-1$
$\\{L_{e},L_{\mu},L_{\tau}\\}$ | $\\{i,i,i\\}$
$\Phi$ | 1
Table 2: Charge assignments for the relevant fields as proposed in case II.
The dark matter field $S$ transforms trivially under the imposed
$\mathbb{Z}_{4}$.

In this case also, we consider an additional discrete symmetry
$\mathbb{Z}_{4}$ and tabulate the charge assignments of the relevant fields in
table 2. After electroweak symmetry breaking one obtains the vertex
$\phi\overline{\nu}\nu$ with the coupling coefficient being proportional to
$\frac{c_{2}v^{2}}{\Lambda^{2}}$ where $v=246$GeV.

##  Appendix B Matrix amplitude squares

Process | $|\mathcal{M}|^{2}/G_{F}^{2}$
---|---
$\nu_{e}\bar{\nu}_{e}\to e^{-}e^{+}$ | $-32(g_{L}^{2}+g_{R}^{2})su$
$\nu_{\mu(\tau)}\bar{\nu}_{\mu(\tau)}\to e^{-}e^{+}$ | $-32\left[(g_{L}-1)^{2}+g_{R}^{2}\right]su$
$\nu_{e}e^{-}\to\nu_{e}e^{-}$ | $32\left[(g_{L}^{2}+g_{R}^{2})s^{2}+\\{(g_{L}-1)^{2}+g_{R}^{2}\\}u^{2}\right]$
$\nu_{\mu(\tau)}e^{-}\to\nu_{\mu(\tau)}e^{-}$ | $32\left[(g_{L}-1)^{2}+g_{R}^{2})\right](s^{2}+u^{2})$
$\nu_{e}e^{+}\to\nu_{e}e^{+}$ | $32\left[(g_{L}^{2}+g_{R}^{2})u^{2}+\\{(g_{L}-1)^{2}+g_{R}^{2}\\}s^{2}\right]$
$\nu_{\mu(\tau)}e^{+}\to\nu_{\mu(\tau)}e^{+}$ | $32\left[(g_{L}-1)^{2}+g_{R}^{2})\right](s^{2}+u^{2})$
$\phi\to\overline{\nu}_{i}\nu_{i}$ | $\frac{2\kappa^{2}m_{\phi}^{2}}{G_{F}^{2}}$
Table 3: Matrix amplitude squares of the relevant processes for neutrino
decoupling. Here $G_{F}$ is the Fermi’s coupling constant, $s,u$ are the
Mandelstam variables, $g_{L}=\dfrac{1}{2}+\sin^{2}\theta_{W}$,
$g_{R}=\sin^{2}\theta_{W}$, and $\theta_{W}$ is the Weinberg Angle. In
deriving these amplitudes, we have neglected the mass of the electron [59,
31].

##  Appendix C Collision Terms

$\bullet$ The collision term for the evolution of the energy density of $A$
for a generic inelastic process such as $A(p_{1})+B(p_{2})\to
C(p_{3})+D(p_{4})$ is given by

$\displaystyle\mathcal{C}^{A}_{\rm inel}=\int
d\Pi_{1}d\Pi_{2}d\Pi_{3}d\Pi_{4}E_{1}(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4})|\mathcal{M}^{2}|_{AB\to
CD}\Lambda_{\rm inel}(E_{1},E_{2},E_{3},E_{4})\,\,,$ (23)

where

$\displaystyle\Lambda_{\rm inel}=f_{C}(E_{3})f_{D}(E_{4})(1\pm
f_{A}(E_{1}))(1\pm f_{B}(E_{2}))-f_{A}(E_{1})f_{B}(E_{2})(1\pm
f_{C}(E_{3}))(1\pm f_{D}(E_{4})).$ (24)

Here the distribution function of the $i^{th}$ species is denoted by $f_{i}$.
$p_{1}$, $p_{2}$, $p_{3}$, $p_{4}$ are four momenta of $A$, $B$, $C$, $D$
respectively. $|\mathcal{M}^{2}|_{AB\to CD}$ is the matrix amplitude square of
$AB\to CD$ process and $d\Pi_{i}=\dfrac{d^{3}\vec{p}_{i}}{2E_{i}(2\pi)^{3}}$.

$\bullet$ The collision term for the evolution of the energy density of $A$
for a generic elastic process such as $A(p_{1})+B(p_{2})\to A(p_{3})+B(p_{4})$
is given by

$\displaystyle\mathcal{C}^{A}_{\rm el}=\dfrac{1}{2}\int
d\Pi_{1}d\Pi_{2}d\Pi_{3}d\Pi_{4}(E_{1}-E_{3})(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4})|\mathcal{M}^{2}|_{AB\to
AB}\Lambda_{\rm el}(E_{1},E_{2},E_{3},E_{4}),$ (25)

where

$\displaystyle\Lambda_{\rm el}=f_{A}(E_{3})f_{B}(E_{4})(1\pm
f_{A}(E_{1}))(1\pm f_{B}(E_{2}))-f_{A}(E_{1})f_{B}(E_{2})(1\pm
f_{A}(E_{3}))(1\pm f_{B}(E_{4}))\,\,.$ (26)

Here $|\mathcal{M}^{2}|_{AB\to AB}$ is the matrix amplitude square of $AB\to
AB$ elastic process.

$\bullet$ For $A(p_{1})\to B(p_{2})+C(p_{3})$ process, the collision term for
the evolution of the energy density of $A$ is as follows:

$\displaystyle\mathcal{C}^{A}_{A\leftrightarrow BC}=\int
d\Pi_{1}d\Pi_{2}d\Pi_{3}E_{1}(2\pi)^{4}\delta^{4}(p_{1}-p_{2}-p_{3})|\mathcal{M}^{2}|_{A\to
BC}\Lambda_{A\leftrightarrow BC}(E_{1},E_{2},E_{3})\,\,,$ (27)

where $|\mathcal{M}^{2}|_{A\to BC}$ is the matrix amplitude square for $A\to
BC$ process and

$\displaystyle\Lambda_{A\leftrightarrow BC}=f_{B}(E_{2})f_{C}(E_{3})(1\pm
f_{A}(E_{1}))-f_{A}(E_{1})(1\pm f_{B}(E_{2}))(1\pm f_{C}(E_{3}))\,\,.$ (28)

## References

  * Aprile  _et al._ [2012] E. Aprile _et al._ (XENON100), Phys. Rev. Lett. 109, 181301 (2012), arXiv:1207.5988 [astro-ph.CO] .
  * Angloher  _et al._ [2016] G. Angloher _et al._ (CRESST), Eur. Phys. J. C 76, 25 (2016), arXiv:1509.01515 [astro-ph.CO] .
  * Agnese  _et al._ [2016] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 116, 071301 (2016), arXiv:1509.02448 [astro-ph.CO] .
  * Klasen  _et al._ [2015] M. Klasen, M. Pohl, and G. Sigl, Prog. Part. Nucl. Phys. 85, 1 (2015), arXiv:1507.03800 [hep-ph] .
  * Tan  _et al._ [2016] A. Tan _et al._ (PandaX-II), Phys. Rev. Lett. 117, 121303 (2016), arXiv:1607.07400 [hep-ex] .
  * Feng [2010] J. L. Feng, Ann. Rev. Astron. Astrophys. 48, 495 (2010), arXiv:1003.0904 [astro-ph.CO] .
  * Roszkowski  _et al._ [2018] L. Roszkowski, E. M. Sessolo, and S. Trojanowski, Rept. Prog. Phys. 81, 066201 (2018), arXiv:1707.06277 [hep-ph] .
  * Schumann [2019] M. Schumann, J. Phys. G 46, 103003 (2019), arXiv:1903.03026 [astro-ph.CO] .
  * Lin [2019] T. Lin, PoS 333, 009 (2019), arXiv:1904.07915 [hep-ph] .
  * Leane [2020] R. K. Leane (2020) pp. 203–228, arXiv:2006.00513 [hep-ph] .
  * 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. A 32, 1730023 (2017), arXiv:1706.07442 [hep-ph] .
  * Pospelov  _et al._ [2008] M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B 662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * Feng  _et al._ [2008] J. L. Feng, H. Tu, and H.-B. Yu, JCAP 10, 043 (2008), arXiv:0808.2318 [hep-ph] .
  * Cheung  _et al._ [2011] C. Cheung, G. Elor, L. J. Hall, and P. Kumar, JHEP 03, 042 (2011), arXiv:1010.0022 [hep-ph] .
  * Chu  _et al._ [2012] X. Chu, T. Hambye, and M. H. G. Tytgat, JCAP 05, 034 (2012), arXiv:1112.0493 [hep-ph] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Mangano  _et al._ [2005] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico, Nucl. Phys. B 729, 221 (2005), arXiv:hep-ph/0506164 .
  * de Salas and Pastor [2016] P. F. de Salas and S. Pastor, JCAP 07, 051 (2016), arXiv:1606.06986 [hep-ph] .
  * Akita and Yamaguchi [2020] K. Akita and M. Yamaguchi, JCAP 08, 012 (2020), arXiv:2005.07047 [hep-ph] .
  * Bennett  _et al._ [2021] J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor, and Y. Y. Y. Wong, JCAP 04, 073 (2021), arXiv:2012.02726 [hep-ph] .
  * Boehm  _et al._ [2012] C. Boehm, M. J. Dolan, and C. McCabe, JCAP 12, 027 (2012), arXiv:1207.0497 [astro-ph.CO] .
  * Brust  _et al._ [2013] C. Brust, D. E. Kaplan, and M. T. Walters, JHEP 12, 058 (2013), arXiv:1303.5379 [hep-ph] .
  * Nollett and Steigman [2015] K. M. Nollett and G. Steigman, Phys. Rev. D 91, 083505 (2015), arXiv:1411.6005 [astro-ph.CO] .
  * Wilkinson  _et al._ [2016] R. J. Wilkinson, A. C. Vincent, C. Bœhm, and C. McCabe, Phys. Rev. D 94, 103525 (2016), arXiv:1602.01114 [astro-ph.CO] .
  * Abazajian and Heeck [2019] K. N. Abazajian and J. Heeck, Phys. Rev. D 100, 075027 (2019), arXiv:1908.03286 [hep-ph] .
  * D’Eramo  _et al._ [2018] F. D’Eramo, R. Z. Ferreira, A. Notari, and J. L. Bernal, JCAP 11, 014 (2018), arXiv:1808.07430 [hep-ph] .
  * Luo  _et al._ [2021] X. Luo, W. Rodejohann, and X.-J. Xu, JCAP 03, 082 (2021), arXiv:2011.13059 [hep-ph] .
  * Luo  _et al._ [2020] X. Luo, W. Rodejohann, and X.-J. Xu, JCAP 06, 058 (2020), arXiv:2005.01629 [hep-ph] .
  * Biswas  _et al._ [2021] A. Biswas, D. Borah, and D. Nanda, JCAP 10, 002 (2021), arXiv:2103.05648 [hep-ph] .
  * Escudero [2019] M. Escudero, JCAP 02, 007 (2019), arXiv:1812.05605 [hep-ph] .
  * Escudero  _et al._ [2019] M. Escudero, D. Hooper, G. Krnjaic, and M. Pierre, JHEP 03, 071 (2019), arXiv:1901.02010 [hep-ph] .
  * Ibe  _et al._ [2020] M. Ibe, S. Kobayashi, Y. Nakayama, and S. Shirai, JHEP 04, 009 (2020), arXiv:1912.12152 [hep-ph] .
  * Escudero Abenza [2020] M. Escudero Abenza, JCAP 05, 048 (2020), arXiv:2001.04466 [hep-ph] .
  * Fields  _et al._ [2020] B. D. Fields, K. A. Olive, T.-H. Yeh, and C. Young, JCAP 03, 010 (2020), [Erratum: JCAP 11, E02 (2020)], arXiv:1912.01132 [astro-ph.CO] .
  * Depta  _et al._ [2019] P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, JCAP 04, 029 (2019), arXiv:1901.06944 [hep-ph] .
  * Ghosh and Sachdeva [2020] D. Ghosh and D. Sachdeva, JCAP 10, 060 (2020), arXiv:2007.01873 [hep-ph] .
  * Giovanetti  _et al._ [2021] C. Giovanetti, M. Lisanti, H. Liu, and J. T. Ruderman, (2021), arXiv:2109.03246 [hep-ph] .
  * Knapen  _et al._ [2017] S. Knapen, T. Lin, and K. M. Zurek, Phys. Rev. D 96, 115021 (2017), arXiv:1709.07882 [hep-ph] .
  * Choudhury and Sachdeva [2019a] D. Choudhury and D. Sachdeva, Phys. Rev. D 100, 035007 (2019a), arXiv:1903.06049 [hep-ph] .
  * Choudhury and Sachdeva [2019b] D. Choudhury and D. Sachdeva, Phys. Rev. D 100, 075007 (2019b), arXiv:1906.06364 [hep-ph] .
  * Boehm  _et al._ [2013] C. Boehm, M. J. Dolan, and C. McCabe, JCAP 08, 041 (2013), arXiv:1303.6270 [hep-ph] .
  * Kitabayashi and Kurosawa [2016] T. Kitabayashi and Y. Kurosawa, Phys. Rev. D 93, 033002 (2016), arXiv:1509.05564 [hep-ph] .
  * Green and Rajendran [2017] D. Green and S. Rajendran, JHEP 10, 013 (2017), arXiv:1701.08750 [hep-ph] .
  * Berlin and Blinov [2018] A. Berlin and N. Blinov, Phys. Rev. Lett. 120, 021801 (2018), arXiv:1706.07046 [hep-ph] .
  * Berlin and Blinov [2019] A. Berlin and N. Blinov, Phys. Rev. D 99, 095030 (2019), arXiv:1807.04282 [hep-ph] .
  * Sabti  _et al._ [2020] N. Sabti, J. Alvey, M. Escudero, M. Fairbairn, and D. Blas, JCAP 01, 004 (2020), arXiv:1910.01649 [hep-ph] .
  * Abazajian  _et al._ [2019] K. Abazajian _et al._ , (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Du  _et al._ [2020] Y. Du, F. Huang, H.-L. Li, and J.-H. Yu, JHEP 12, 207 (2020), arXiv:2005.01717 [hep-ph] .
  * Hannestad [2004] S. Hannestad, Phys. Rev. D 70, 043506 (2004), arXiv:astro-ph/0403291 .
  * Hasegawa  _et al._ [2019] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram, and S. Hannestad, JCAP 12, 012 (2019), arXiv:1908.10189 [hep-ph] .
  * Hannestad and Madsen [1995] S. Hannestad and J. Madsen, Phys. Rev. D 52, 1764 (1995), arXiv:astro-ph/9506015 .
  * Dolgov [2002] A. D. Dolgov, Phys. Rept. 370, 333 (2002), arXiv:hep-ph/0202122 .
  * Kreisch  _et al._ [2020] C. D. Kreisch, F.-Y. Cyr-Racine, and O. Doré, Phys. Rev. D 101, 123505 (2020), arXiv:1902.00534 [astro-ph.CO] .
  * Arbey  _et al._ [2020] A. Arbey, J. Auffinger, K. P. Hickerson, and E. S. Jenssen, Comput. Phys. Commun. 248, 106982 (2020), arXiv:1806.11095 [astro-ph.CO] .
  * An  _et al._ [2022] R. An, V. Gluscevic, E. Calabrese, and J. C. Hill, JCAP 07, 002 (2022), arXiv:2202.03515 [astro-ph.CO] .
  * Kersten and Smirnov [2007] J. Kersten and A. Y. Smirnov, Phys. Rev. D 76, 073005 (2007), arXiv:0705.3221 [hep-ph] .
  * Kouvaris  _et al._ [2015] C. Kouvaris, I. M. Shoemaker, and K. Tuominen, Phys. Rev. D 91, 043519 (2015), arXiv:1411.3730 [hep-ph] .
  * Kawasaki  _et al._ [2000] M. Kawasaki, K. Kohri, and N. Sugiyama, Phys. Rev. D 62, 023506 (2000), arXiv:astro-ph/0002127 .


# On complexity constants of linear and quadratic models for derivative-free
trust-region algorithms111This study was financed in part by the Coordenação
de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance
Code 001

A. E. Schwertner222Graduate Program in Mathematics, State University of
Maringá, Paraná, Brazil F. N. C. Sobral333Department of Mathematics, State
University of Maringá, Paraná, Brazil 444 Corresponding author,
<EMAIL_ADDRESS>
###### Abstract

Complexity analysis has become an important tool in the convergence analysis
of optimization algorithms. For derivative-free optimization algorithms, it is
not different. Interestingly, several constants that appear when developing
complexity results hide the dimensions of the problem. This work organizes
several results in literature about bounds that appear in derivative-free
trust-region algorithms based on linear and quadratic models. All the
constants are given explicitly by the quality of the sample set, dimension of
the problem and number of sample points. We extend some results to allow
“inexact” interpolation sets. We also provide a clearer proof than those
already existing in literature for the underdetermined case.

Keywords: Polynomial interpolation, Derivative-free trust-region algorithms,
Minimum Frobenius norm, Error bounds

##  1 Introduction

Polynomial interpolation is one of the main aspects of model-based derivative-
free algorithms. Given a function $f:\mathbb{R}^{n}\to\mathbb{R}$ and a set of
interpolation points $\mathcal{Y}=\\{y^{0},\dots,y^{p}\\}$, the interpolating
polynomial $\mathrm{m}\in\mathcal{P}_{n}^{a}(\mathbb{R})$ is such that

$\mathrm{m}(y^{i})=f(y^{i}),\quad i=0,\dots,p,$ (1)

where $\mathcal{P}_{n}^{a}(\mathbb{R})$ is the space of polynomials of degree
at most $a$ in $n$ variables. It is well understood that quadratic polynomials
($a=2$) provides a good approximation, being able to capture the curvature of
$f$ using a reasonable amount of interpolation points [5, p.35].

More recently, the wort-case complexity analysis of optimization algorithms
became popular. Although most of the bounds obtained are very pessimistic, the
analysis provide further insights into the parameters and the dependence on
the problems’ dimensions.

When studying worst-case complexity of derivative-free trust-region algorithms
using linear or quadratic polynomials, one has to deal with bounds related to
the geometric quality of $\mathcal{Y}$ and to the Hessian of the model (in the
quadratic case). In [5], several bounds were provided, which explicitly show
the dependence on the problems’ dimensions. Such bounds were used, for
example, in the excellent complexity analysis of [2, 6]. The complexity bounds
generated the desire of performing derivative-free optimization in smaller
subspaces, to further decrease the impact of dimensionality in the
construction of models [3].

In [9], some error bounds were obtained for a relaxed version of (1)

$|\mathrm{m}(y^{i})-f(y^{i})|\leq\kappa\delta^{2},$ (2)

where $\kappa$ is a constant and $\delta$ can be viewed as the precision of
the model [8, 9] or just as the trust-region radius [5]. Condition (2) is
related to linear and quadratic models build from support vector regression
strategies. Unfortunately, only the determined case is analyzed in [9], that
is, the case where $p=dim(\mathcal{P}_{n}^{a}(\mathbb{R}))-1$. In the case of
quadratic models, that would require $O(n^{2})$ interpolation points in
$\mathcal{Y}$. This work intends to extend the bounds obtained in [9] to the
underdetermined case, which is known to have better performance in practical
implementations.

The main contributions of the present work are

  * •

the bounds associated with linear and quadratic polynomial models are
organized, making clear their dependence on the dimension and quality of the
interpolation set;

  * •

the bounds from [9] are extended to the underdetermined case.

This work is organized as follows.

Section 2 is concerned with bounds to linear and quadratic polynomial models
which are uniquely determined by the sample set. Section 3 is the main
contribution of this work and discusses error bounds for underdetermined
quadratic models, since they are the most efficient models in practice [7].
Finally, in Section 4 we state and organize all the bounds and discuss
possible open questions.

Throughout this work, we will refer to the interpolating polynomial
$\mathrm{m}$ as the one satisfying the relaxed condition (2),
$\overline{B}(x,\delta)$ is the closed ball centered at $x$ with radius
$\delta$, $\left\lVert\cdot\right\lVert$ is the Euclidean vector norm or its
respective induced norm for matrices, $\mathcal{Y}=\\{y^{0},\dots,y^{p}\\}$ is
the set of $p+1$ sample points, for which $f$ is assumed to be known,
$\mathcal{P}_{n}^{a}(\mathbb{R})$ is the space of polynomials of degree at
most $a$ in $n$ variables, and $q=dim(\mathcal{P}_{n}^{a}(\mathbb{R}))-1$.

##  2 Determined models

Let us assume that $\mathcal{Y}\subset\overline{B}(y^{0},\delta)$. Determined
interpolation models need that the number of sample points in $\mathcal{Y}$
equals the dimension of $\mathcal{P}_{n}^{a}(\mathbb{R})$, that is, $p=q$ in
our notation. In the linear case, we assume that $p=n=q$ and define matrices

$\mathbf{L}_{L}=\left[\begin{array}[]{c}(y^{1}-y^{0})^{T}\\\ \vdots\\\
(y^{n}-y^{0})^{T}\end{array}\right]=\left[\begin{array}[]{ccc}y_{1}^{1}-y_{1}^{0}&\dots&y_{n}^{1}-y_{n}^{0}\\\
\vdots&\ddots&\vdots\\\
y_{1}^{n}-y_{1}^{0}&\dots&y_{n}^{n}-y_{n}^{0}\end{array}\right],$ (3)

and

$\hat{\mathbf{L}}_{L}=\frac{1}{\delta}\mathbf{L}_{L}.$ (4)

In the quadratic case, we assume $p=(n^{2}+3n)/2=q$ and define the matrix

$\mathbf{L}_{Q}=\left[\begin{array}[]{c}\big{(}\overline{\varphi}(y^{1}-y^{0})\big{)}^{T}\\\
\vdots\\\
\big{(}\overline{\varphi}(y^{q}-y^{0})\big{)}^{T}\end{array}\right],$ (5)

where
$\overline{\varphi}(x)=\bigg{[}x_{1},x_{2},\dots,x_{n},\frac{1}{2}x_{1}^{2},x_{1}x_{2},x_{1}x_{3},\dots,x_{1}x_{n},\frac{1}{2}x_{2}^{2},\dots,x_{n-1}x_{n},\frac{1}{2}x_{n}^{2}\bigg{]}^{T}$
is the vector whose elements form the natural basis of monomials in
$\mathcal{P}_{n}^{2}(\mathbb{R})$, as defined in [5, Section 3.1]. We also
consider the matrix

$\hat{\mathbf{L}}_{Q}=\mathbf{L}_{Q}\left[\begin{array}[]{cc}\mathbf{D}_{L}^{-1}&0\\\
0&\mathbf{D}_{Q}^{-1},\end{array}\right]$ (6)

where $\mathbf{D}_{L}=\delta\mathbf{I}_{n\times n}$ e
$\mathbf{D}_{Q}=\delta^{2}\mathbf{I}_{(q-n)\times(q-n)}$. Matrices
$\mathbf{L}_{L}$ and $\mathbf{L}_{Q}$ are related with the construction of
determined interpolation models, while (4) and (6) are their respective scaled
versions, mainly used for theoretical purposes.

We will assume that the following assumptions are valid.

######  Assumption 2.1.

$\nabla f$ is Lipschitz continuous with constant $L$ in a sufficiently large
open bounded domain $\mathcal{X}\subset\mathbb{R}^{n}$.

######  Assumption 2.2.

$f$ is bounded below.

######  Assumption 2.3.

If the model to be built is linear, then the matrix $\mathbf{L}_{L}$ is non-
singular and there is a constant $\kappa_{L}>0$ such that
$\left\lVert\hat{\mathbf{L}}_{L}^{-1}\right\lVert\leq\kappa_{L}$. If the model
is quadratic, then the matrix $\mathbf{L}_{Q}$ is nonsingular and there is a
constant $\kappa_{Q}>0$ such that
$\left\lVert\hat{\mathbf{L}}_{Q}^{-1}\right\lVert\leq\kappa_{Q}$.

The next assumption is related to the relaxed concept of interpolation models
considered in this work.

######  Assumption 2.4.

There is a constant $\kappa\geq 0$ such that for every
$y^{j}\in\mathcal{Y}\subset\overline{B}(y^{0},\delta)$

$\left\lvert\mathrm{m}(y^{j})-f{(y^{j})}\right\lvert\leq\kappa\delta^{2}.$

Note that Assumption 2.4 includes the creation of “exact” interpolation models
and also models under uncertainty [9], under numerical errors or even when we
can control the precision in which $f$ is calculated [1]. The next theorem
provides the error bounds for determined polynomial interpolation under
Assumption 2.4. Its proof is provided in [9].

######  Theorem 2.5.

Suppose that Assumptions 2.1, 2.2, 2.3 and 2.4 hold. Then, for all
$x\in\mathcal{X}\cap\overline{B}(y^{0},\delta)$, if the model is linear

$\displaystyle\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert$
$\displaystyle=0,$ $\displaystyle\left\lVert\nabla
f(x)-\nabla\mathrm{m}(x)\right\lVert$
$\displaystyle\leq\left(L+\left(\frac{1}{2}L+2\kappa\right)\kappa_{L}\sqrt{n}\right)\delta,$
$\displaystyle\left\lvert f(x)-\mathrm{m}(x)\right\lvert$
$\displaystyle\leq\left(\frac{1}{2}L+\kappa+\left(\frac{1}{2}L+2\kappa\right)\kappa_{L}\sqrt{n}\right)\delta^{2},$

and, if the model is quadratic,

$\displaystyle\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert$
$\displaystyle\leq 2\kappa_{Q}\sqrt{2q}(\kappa+L),$
$\displaystyle\left\lVert\nabla f(x)-\nabla\mathrm{m}(x)\right\lVert$
$\displaystyle\leq\left(2\kappa_{Q}\sqrt{q}\left(1+\sqrt{2}\right)(\kappa+L)\right)\delta,$
$\displaystyle\left\lvert f(x)-\mathrm{m}(x)\right\lvert$
$\displaystyle\leq\left(\frac{1}{2}L+\kappa+\kappa_{Q}\sqrt{q}\left(2+3\sqrt{2}\right)(\kappa+L)\right)\delta^{2}.$

In order to ensure Assumption 2.3, we need to make additional hypotheses about
the geometry of the sample set. The existence and uniqueness of model m, also
known as _poisedness_ , is not enough. In linear interpolation, for example,
this can be achieved by assuming that the points in $\mathcal{Y}$ are affinely
linear independent [10]. Here, we assume that the set $\mathcal{Y}$ is
$\Lambda$-poised on the ball $\overline{B}(y^{0},\delta)$, for some
$\Lambda>0$. This definition is taken from [5] and will be better discussed in
Section 3 for underdetermined polynomials. The following lemma, adapted from
[9], shows that Assumption 2.3 is valid when the models to be built are linear
or quadratic and $\mathcal{Y}$ is $\Lambda$-poised.

######  Lemma 2.6.

Let $\Lambda>0$. If $\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{n}\\}$ is
$\Lambda$-poised in $\overline{B}(y^{0},\delta)$ with respect to the basis
$\phi$ of $\mathcal{P}^{1}_{n}$, then
$\left\lVert\hat{\mathbf{L}}_{L}^{-1}\right\lVert\leq\Lambda\sqrt{n}$. If
$\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{q}\\}$ is a $\Lambda$-poised set in
$\overline{B}(y^{0},\delta)$ with respect to the basis $\phi$ in
$\mathcal{P}^{2}_{n}$, then,
$\left\lVert\hat{\mathbf{L}}_{Q}^{-1}\right\lVert\leq
4\Lambda\sqrt{(q+1)^{3}}$.

##  3 Underdetermined models

The main drawback of determined polynomials is the high number of
interpolation points needed in the quadratic case. One could use only linear
polynomials, but quadratic polynomials are known to better explore the
curvature of $f$. Therefore, we are interested in building quadratic
polynomials using less than $q$ interpolation points. Now we assume that
$n<p<q$ and $q=dim(\mathcal{P}_{n}^{2}(\mathbb{R}))-1=(n^{2}+3n)/2$.

When dealing with classical interpolation theory, in the sense of (1),
building an underdetermined interpolation quadratic model
$\mathrm{m}\in\mathcal{P}_{n}^{2}(\mathbb{R})$ can be viewed as finding
$\alpha\in\mathbb{R}^{q+1}$ such that

$\sum_{j=0}^{q}\alpha_{j}\phi_{j}(y^{i})=f(y^{i})\text{, for }i=0,\dots,p,$

or, equivalently,

$\textbf{M}(\phi,\mathcal{Y})\alpha=f(\mathcal{Y}),$ (7)

where $\phi=\\{\phi_{0}(x),\phi_{1}(x),\dots,\phi_{q}(x)\\}$ is a basis for
$\mathcal{P}_{n}^{2}(\mathbb{R})$,

$\textbf{M}(\phi,\mathcal{Y})=\left[\begin{array}[]{cccc}\phi_{0}(y^{0})&\phi_{1}(y^{0})&\cdots&\phi_{q}(y^{0})\\\
\phi_{0}(y^{1})&\phi_{1}(y^{1})&\cdots&\phi_{q}(y^{1})\\\
\vdots&\vdots&\ddots&\vdots\\\
\phi_{0}(y^{p})&\phi_{1}(y^{p})&\cdots&\phi_{q}(y^{p})\end{array}\right],\text{
and }f(\mathcal{Y})=\left[\begin{array}[]{c}f(y^{0})\\\ f(y^{1})\\\ \vdots\\\
f(y^{p})\end{array}\right].$ (8)

Under the above conditions, the number of interpolation points is less than
the dimension of the space of the polynomials
$\mathcal{P}_{n}^{2}(\mathbb{R})$, that is,
$\left\lvert\mathcal{Y}\right\lvert=p+1<q+1=dim(\mathcal{P}_{n}^{2}(\mathbb{R}))$.
Therefore, the interpolation polynomials defined by (7) are no longer unique.

In this section, our goal is to build an underdetermined quadratic model,
which approximately interpolates the function $f$ over
$\mathcal{Y}\subset\overline{B}(y^{0},\delta)$, in the sense of Assumption
2.4. First, let us consider matrices

$\mathbf{L}_{s}=\left[\begin{array}[]{c}(y^{1}-y^{0})^{T}\\\ \vdots\\\
(y^{p}-y^{0})^{T}\end{array}\right]=\left[\begin{array}[]{ccc}y_{1}^{1}-y_{1}^{0}&\dots&y_{n}^{1}-y_{n}^{0}\\\
\vdots&\ddots&\vdots\\\
y_{1}^{p}-y_{1}^{0}&\dots&y_{n}^{p}-y_{n}^{0}\end{array}\right],$ (9)

and

$\hat{\mathbf{L}}_{s}=\frac{1}{\delta}\mathbf{L}_{s}.$ (10)

In order to obtain the main error bounds for the underdetermined case, we need
to assume some geometric conditions. In other words, we need the sample set
$\mathcal{Y}$ to be poised for linear regression and the Hessian of the model
to be bounded. Such properties are given by assumptions 3.1 and 3.2.

######  Assumption 3.1.

The matrix $\mathbf{L}_{s}\in\mathbb{R}^{p\times n}$ has full column rank,
that is, $rank(\mathbf{L}_{s})=n$, and there is a constant $\kappa_{s}>0$ such
that $\left\lVert\hat{\mathbf{L}}^{\dagger}_{s}\right\lVert\leq\kappa_{s}$,
where $\hat{\mathbf{L}}^{\dagger}_{s}$ denotes the Moore-Penrose pseudo-
inverse of the matrix $\hat{\mathbf{L}}_{s}$.

######  Assumption 3.2.

There is a constant $\kappa_{H}\geq 0$ such that
$\left\lVert\mathbf{H}\right\lVert\leq\kappa_{H}$.

The following theorem establishes error bounds for the underdetermined model
and its derivatives in $\mathcal{X}\cap\overline{B}(y^{0},\delta)$.

######  Theorem 3.3.

Suppose that Assumptions 2.1, 2.2, 3.1 and 3.2 hold. Then, for all
$x\in\mathcal{X}\cap\overline{B}(y^{0},\delta)$,

$\displaystyle\left\lvert f(x)-\mathrm{m}(x)\right\lvert$
$\displaystyle\leq\bigg{(}\frac{1}{2}\big{(}L+\kappa_{H}\big{)}+\kappa+2\kappa_{s}\sqrt{p}\bigg{(}L+\kappa+\frac{3}{4}\kappa_{H}\bigg{)}\bigg{)}\delta^{2},$
(11) $\displaystyle\left\lVert\nabla f(x)-\nabla\mathrm{m}(x)\right\lVert$
$\displaystyle\leq
2\kappa_{s}\sqrt{p}\bigg{(}L+\kappa+\frac{3}{4}\kappa_{H}\bigg{)}\delta.$ (12)

###### Proof.

The proof is strongly based in the ideas of [5, Theorem 5.4] and [9, Theorem
2]. In addition $\left\lVert\mathbf{H}\right\lVert$ and
$\left\lVert\hat{\mathbf{L}}^{\dagger}_{s}\right\lVert$ were replaced by their
respective bounds in equations (11) and (12). ∎

###  3.1 Minimum Frobenius norm models

To establish a bound for the norm of the Hessian of the model, we must specify
its construction. In this sense, we will reduce the degree of freedom in
choosing the model, assuming that it is the minimum Frobenius norm model. We
will also need to define what is a relaxed minimum Frobenius norm model.

Let $\overline{\phi}$ be the natural basis of
$\mathcal{P}^{2}_{n}(\mathbb{R})$ and consider its split into the sets
$\overline{\phi}_{L}=\\{1,x_{1},\dots,x_{n}\\}$ and
$\overline{\phi}_{Q}=\big{\\{}\frac{1}{2}x_{1}^{2},x_{1}x_{2},\dots,\frac{1}{2}x_{n}^{2}\big{\\}}$.
Then we can write the underdetermined interpolation model as

$\mathrm{m}(x)=\alpha_{L}^{T}\overline{\phi}_{L}(x)+\alpha_{Q}^{T}\overline{\phi}_{Q}(x),$

where $\alpha_{L}$ and $\alpha_{Q}$ are appropriate parts of the vector of
coefficients $\alpha$. As in [5, p. 80-81], we can define the minimum
Frobenius norm solution $\alpha^{\textrm{mfn}}$ as the solution of the
optimization problem,

$\begin{array}[]{cc}\min&\frac{1}{2}\left\lVert\alpha_{Q}\right\lVert^{2}\\\
\text{s.t.}&\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})\alpha_{L}+\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})\alpha_{Q}=f(\mathcal{Y}),\end{array}$
(13)

in the variables $\alpha_{L}$ and $\alpha_{Q}$, where $\mathcal{Y}$ is the set
of interpolation points and the matrices
$\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})$ and
$\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})$ are submatrices of
$\mathbf{M}(\overline{\phi},\mathcal{Y})$, given in (8), whose columns
correspond to the elements of $\overline{\phi}_{L}$ and $\overline{\phi}_{Q}$,
respectively.

It is important to point out that the condition of existence and uniqueness of
the minimum Frobenius norm model is that the matrix defined by

$\mathbf{F}(\overline{\phi},\mathcal{Y})=\left[\begin{array}[]{cc}\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})^{T}&\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})\\\
\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})^{T}&0\end{array}\right]$

is nonsingular. Note that $\mathbf{F}(\overline{\phi},\mathcal{Y})$ is
nonsingular if and only if $\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})$ has a
full column rank and
$\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})^{T}$
is positive definite on the nullspace of
$\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})^{T}$, this last condition being
guaranteed if $\mathbf{M}(\overline{\phi}_{Q},\mathcal{Y})$ has full row rank
[5, p. 81].

We say that $\mathcal{Y}$ is poised in the minimum Frobenius norm sense if the
matrix $\mathbf{F}(\overline{\phi},\mathcal{Y})$ is nonsingular. Note that
poisedness in the minimum Frobenius norm sense implies poisedness in the
linear interpolation or regression senses and, as a result, poisedness for
underdetermined quadratic interpolation in the minimum norm sense [5, p. 81].

######  Remark 3.4.

As we saw earlier, if $\mathcal{Y}$ is poised in the minimum Frobenius norm
sense, then $\mathbf{F}(\overline{\phi},\mathcal{Y})$ is nonsingular and hence
$\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})$ has a full column rank. So, as

$\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})=\begin{bmatrix}1&0\\\
e&\mathbf{Id}\end{bmatrix}\left[\begin{array}[]{cc}1&{y^{0}}^{T}\\\
0&\mathbf{L}_{s}\end{array}\right]=\mathbf{E}^{-1}\left[\begin{array}[]{cc}1&{y^{0}}^{T}\\\
0&\mathbf{L}_{s}\end{array}\right]$

where $e=[1,1,\dots,1]^{T}\in\mathbb{R}^{p}$, and
$\mathbf{Id}\in\mathbb{R}^{p\times p}$ is the identity matrix. Since
$\mathbf{E}$ is an elementary nonsingular matrix and
$\mathbf{M}(\overline{\phi}_{L},\mathcal{Y})$ has full column rank, it follows
that $\mathbf{L}_{s}$ has full column rank, which partially satisfies
Assumption 3.1.

Therefore, we will assume that the following hypothesis holds.

######  Assumption 3.5.

The set of sample points
$\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{p}\\}\subset\mathbb{R}^{n}$ is poised in
the minimum Frobenius norm sense on $\overline{B}(y^{0},\delta)$.

Consider the following auxiliary lemma.

######  Lemma 3.6.

There is a number $\sigma_{\infty}>0$ such that for any choice of $v$
satisfying $\left\lVert v\right\lVert_{\infty}=1$, there is
$z\in\overline{B}(0,1)$ such that $\left\lvert
v^{T}\phi(z)\right\lvert\geq\sigma_{\infty}$. If, in addition,
$v^{T}\overline{\phi}(x)$ is a quadratic polynomial and $\overline{\phi}$ is
the natural basis in $\mathcal{P}_{n}^{2}(\mathbb{R})$, then

$\max_{x\in\overline{B}(0,1)}\left\lvert
v^{T}\overline{\phi}(x)\right\lvert\geq\frac{1}{4}.$

If $\overline{\phi}$ is the natural basis of
$\mathcal{P}_{n}^{1}(\mathbb{R})$, then,

$\max_{x\in\overline{B}(0,1)}\left\lvert
v^{T}\overline{\phi}(x)\right\lvert\geq 1.$

###### Proof.

See lemmas 3.10, 3.11, and 6.7 in [5]. ∎

######  Remark 3.7.

Let $\overline{\phi}$ be the natural basis of
$\mathcal{P}_{n}^{1}(\mathbb{R})$. Given $\overline{v}\in\mathbb{R}^{n+1}$
with $\left\lVert\overline{v}\right\lVert=1$, there are
$\beta\in(0,\sqrt{n+1})$ and $v\in\mathbb{R}^{n+1}$ such that $\left\lVert
v\right\lVert_{\infty}=1$ and $v=\beta\overline{v}$. Thus, by Lemma 3.6, we
have that

$\max_{x\in\overline{B}(0,1)}\left\lvert\overline{v}^{T}\overline{\phi}(x)\right\lvert=\max_{x\in\overline{B}(0,1)}\frac{1}{\beta}\left\lvert
v^{T}\overline{\phi}(x)\right\lvert\geq\frac{1}{\sqrt{n+1}}\max_{x\in\overline{B}(0,1)}\left\lvert
v^{T}\overline{\phi}(x)\right\lvert\geq\frac{1}{\sqrt{n+1}}.$

Equation (13) is strongly related with the exact interpolation condition (1).
In order to relax such condition we will first define minimum Frobenius norm
Lagrange polynomials.

######  Definition 3.8.

Given a set $\mathcal{Y}$ of $p+1$ interpolation points, the set of $p+1$
polynomials $\ell_{j}(x)=\alpha_{j}^{T}\overline{\phi}(x)$, $j=0,\dots,p$, is
called set of minimum Frobenius norm Lagrange polynomials for
$\overline{\phi}$, if $\alpha_{j}=((\alpha_{j})_{L},(\alpha_{j})_{Q})$ is the
solution of

$\begin{array}[]{ll}\min&\frac{1}{2}\left\lVert\alpha_{Q}\right\lVert^{2}\\\
\mbox{s.
t.}&M(\overline{\phi}_{L},\mathcal{Y})\alpha_{L}+M(\overline{\phi}_{Q},\mathcal{Y})\alpha_{Q}=e_{j+1},\end{array}$

where $e_{j}$ is the $j$-th canonical vector in $\mathbb{R}^{p+1}$.

By comparing the KKT conditions of Definition 3.8 and equation (13) it is
possible to show that the minimum Frobenius norm polynomial that interpolates
$f$ in $\mathcal{Y}$ can be defined as
$\mathrm{m}(x)=\sum_{j=1}^{p+1}\ell_{j}(x)f(y^{j}).$ So, in order to relax
(13), we will say that our relaxed minimum Frobenius norm polynomials are
given by

$\mathrm{m}(x)=\sum_{j=1}^{p+1}\ell_{j}(x)\gamma_{j}$ (14)

where $\gamma_{j}$, $j=1,\dots,p+1$, are chosen such that Assumption 2.4
holds. One can think that $\mathrm{m}$ was obtained by solving (13) with
$f(y^{j})$ replaced by $\gamma_{j}$ for all $j=0,\dots,p$.

The theorems that follow are the main contribution of this work. Their
establish bounds for $\left\lVert\nabla^{2}m(x)\right\lVert$ and
$\left\lVert\hat{L}^{\dagger}_{s}\right\lVert$ under the hypothesis of well
poised interpolation sets and, therefore, assumptions 3.1 and 3.2 can be
removed from Theorem 3.3. Using a partial definition from [5], we first define
what we mean as well poised interpolation sets for minimum Frobenius norm
quadratic models.

######  Definition 3.9.

Let $\Lambda>0$ and $B\subseteq\mathbb{R}^{n}$ be given. Let $\overline{\phi}$
be the natural basis of monomials of $\mathcal{P}_{n}^{2}(\mathbb{R})$. A
poised set $\mathcal{Y}$ is said to be $\Lambda$-poised in $B$ in the minimum
Frobenius norm sense if and only if for any $x\in B$, the solution
$\lambda(x)\in\mathbb{R}^{p+1}$ of

$\begin{array}[]{ll}\min&\frac{1}{2}\left\lVert
M(\overline{\phi}_{Q},\mathcal{Y})^{T}\lambda(x)-\overline{\phi}_{Q}(x)\right\lVert^{2}\\\
\mbox{s.t.}&M(\overline{\phi}_{L},\mathcal{Y})^{T}\lambda(x)=\overline{\phi}_{L}(x)\end{array}$

is such that $\left\lVert\lambda(x)\right\lVert_{\infty}\leq\Lambda$.

Note that Definition 3.9 does not need to be relaxed. We are now ready to
provide the necessary bounds in terms of the quality of the interpolation set,
the Lipschitz constant and the dimensions of the problem.

######  Theorem 3.10.

Suppose that Assumptions 2.1 and 3.5 hold. Assume that the set
$\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{p}\\}$ is $\Lambda$-poised in
$\overline{B}(y^{0},\delta)$ in the minimum Frobenius norm sense and
$\delta_{\max}>0$ is an upper bound for $\delta$. Then,

$\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert\leq\left(\kappa+\frac{L}{2}\right)\frac{4\Lambda(p+1)\sqrt{2(q+1)}}{c(\delta_{\max})^{2}},$

where $c(\delta_{\max})=\min\\{1,1/\delta_{\max},1/\delta_{\max}^{2}\\}$.

###### Proof.

This proof follows the proof of [5, Theorem 5.7]. Assume, without loss of
generality, that $y^{0}=0$. By Lemma 3.6, the definition of
$\Lambda$-poisedness and arguments very similar to [5, Theorem 5.7] we have
that

$\left\lVert\nabla^{2}\ell_{j}(x)\right\lVert\leq\sqrt{2(q+1)}\frac{4\Lambda}{\delta^{2}c(\delta_{\max})^{2}},\quad
j=0,\dots,p,$ (15)

where $\delta_{\max}$ and $c(\delta_{\max})$ were defined by the theorem.

Now, let us consider the function

$\hat{f}(x)=f(x)-f(y^{0})-\nabla f(y^{0})^{T}(x-y^{0}).$

Note that $\hat{f}(y^{0})=0$, $\nabla\hat{f}(y^{0})=0$ and the Hessian remains
unchanged. In addition, if $\mathrm{m}(x)$ is a relaxed minimum Frobenius norm
model for $f$ over the set $\mathcal{Y}$, then it is easy to see that

$\hat{\mathrm{m}}(x)=\mathrm{m}(x)-f(y^{0})-\nabla f(y^{0})^{T}(x-y^{0})$

also satisfies Assumption 2.4 for $\hat{f}$ and the points in $\mathcal{Y}$.
Therefore, we can assume without loss of generality that $f(y^{0})=0$ and
$\nabla f(y^{0})=0$. Thus, we have that $\left\lvert
f(x)\right\lvert=\left\lvert f(x)-f(y^{0})-\nabla
f(y^{0})^{T}(x-y^{0})\right\lvert\leq(L/2)\delta^{2}$, for all $x\in
B(y^{0},\delta)$ and, by (14), we can prove that

$\left\lvert\gamma_{j}\right\lvert\leq(\kappa+L/2)\delta^{2},\quad
j=0,\dots,p.$ (16)

By using the definition of $\mathrm{m}(x)$ in (14) and bounds (15) and (16),
we finally get that

$\begin{split}\left\lVert\mathbf{H}\right\lVert&=\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert=\left\lVert\sum_{j=0}^{p}\gamma_{j}\nabla^{2}\ell_{j}(x)\right\lVert\leq\sum_{j=0}^{p}\left\lvert\gamma_{j}\right\lvert\left\lVert\nabla^{2}\ell_{j}(x)\right\lVert\\\
&\leq\left(\kappa+\frac{L}{2}\right)\frac{4\Lambda(p+1)\sqrt{2(q+1)}}{c(\delta_{\max})^{2}},\end{split}$

and the proof is complete. ∎

######  Theorem 3.11.

If $\mathcal{Y}=\\{y^{0},y^{1},\dots,y^{p}\\}$ is a $\Lambda$-poised set in
$\overline{B}(y^{0},\delta)$ in the sense of the minimum Frobenius norm, then

$\left\lVert\hat{\mathbf{L}}_{s}^{\dagger}\right\lVert\leq\Lambda\sqrt{2(n+1)}(p+1).$
(17)

###### Proof.

It is known that $\Lambda$-poisedness does not depend on the scale of the
sample set and it is invariant with respect to shifts [5]. Therefore, let us
consider the set

$\hat{\mathcal{Y}}=\bigg{\\{}0,\frac{y^{1}-y^{0}}{\delta},\dots,\frac{y^{p}-y^{0}}{\delta}\bigg{\\}},$

which is $\Lambda$-poised in $\overline{B}(0,1)$, and the matrices given by

$\hat{\mathbf{M}}_{L}=\mathbf{M}(\overline{\phi}_{L},\hat{\mathcal{Y}})=\left[\begin{array}[]{cc}1&0\\\
e&\hat{\mathbf{L}}_{s}\end{array}\right],\mathbf{E}=\left[\begin{array}[]{cc}1&0\\\
-e&\mathbf{Id}\end{array}\right],\mathbf{E}^{-1}=\left[\begin{array}[]{cc}1&0\\\
e&\mathbf{Id}\end{array}\right],\text{ and }$
$\mathbf{Q}=\mathbf{E}\hat{\mathbf{M}}_{L}=\left[\begin{array}[]{cc}1&0\\\
0&\hat{\mathbf{L}}_{s}\end{array}\right],$

where $\hat{\mathbf{M}}_{L}\in\mathbb{R}^{(p+1)\times(n+1)}$,
$\mathbf{E}\in\mathbb{R}^{(p+1)\times(p+1)}$,
$\mathbf{E}^{-1}\in\mathbb{R}^{(p+1)\times(p+1)}$,
$\mathbf{Q}\in\mathbb{R}^{(p+1)\times(n+1)}$,
$e=[1,1,\dots,1]^{T}\in\mathbb{R}^{p}$, $\mathbf{Id}\in\mathbb{R}^{p\times p}$
is the identity matrix and $\hat{\mathbf{L}}_{s}$ was defined in (10). Note
that the Moore-Penrose pseudoinverse of $\mathbf{Q}$ is given by

$\mathbf{Q}^{\dagger}=\left[\begin{array}[]{cc}1&0\\\
0&\hat{\mathbf{L}}_{s}^{\dagger}\end{array}\right],$

where $\mathbf{Q}^{\dagger}\in\mathbb{R}^{(n+1)\times(p+1)}$. Thus,

$\left\lVert\mathbf{Q}^{\dagger}\right\lVert=\max\Big{\\{}1,\left\lVert\hat{\mathbf{L}}_{s}^{\dagger}\right\lVert\Big{\\}}\geq\left\lVert\hat{\mathbf{L}}_{s}^{\dagger}\right\lVert.$
(18)

By Definition 3.9, for every $x\in\overline{B}(0,1)\subset\mathbb{R}^{n}$
there exists $\lambda(x)\in\mathbb{R}^{p+1}$, with
$\left\lVert\lambda(x)\right\lVert_{\infty}\leq\Lambda$, such that
$\hat{\mathbf{M}}_{L}^{T}\lambda(x)=\overline{\phi}_{L}(x)$. Since
$\mathbf{E}$ is a non-singular matrix, we have that for each
$x\in\overline{B}(0,1)$,

$\displaystyle\hat{\mathbf{M}}_{L}^{T}\lambda(x)=\overline{\phi}_{L}(x)$
$\displaystyle\Longleftrightarrow\hat{\mathbf{M}}_{L}^{T}\big{(}\mathbf{E}^{T}\mathbf{E}^{-T}\big{)}\lambda(x)=\overline{\phi}_{L}(x)$
(19)
$\displaystyle\Longleftrightarrow\big{(}\hat{\mathbf{M}}_{L}^{T}\mathbf{E}^{T}\big{)}\big{(}\mathbf{E}^{-T}\lambda(x)\big{)}=\overline{\phi}_{L}(x)$
$\displaystyle\Longleftrightarrow\big{(}\mathbf{E}\hat{\mathbf{M}}_{L}\big{)}^{T}\big{(}\mathbf{E}^{-T}\lambda(x)\big{)}=\overline{\phi}_{L}(x)$
$\displaystyle\Longleftrightarrow\mathbf{Q}^{T}\big{(}\mathbf{E}^{-T}\lambda(x)\big{)}=\overline{\phi}_{L}(x)$
$\displaystyle\Longleftrightarrow\mathbf{E}^{-T}\lambda(x)=\mathbf{Q}^{\dagger^{T}}\overline{\phi}_{L}(x).$

Note that $\left\lVert\lambda(x)\right\lVert_{\infty}\leq\Lambda$ and
$\mathbf{E}^{-T}\lambda(x)=\Bigg{[}\sum_{j=0}^{p}[\lambda(x)]_{j},\lambda(x)^{T}\Bigg{]}^{T}$.
Thus, from (19), we get

$\displaystyle\left\lVert\mathbf{Q}^{\dagger^{T}}\overline{\phi}_{L}(x)\right\lVert$
$\displaystyle=\left\lVert\mathbf{E}^{-T}\lambda(x)\right\lVert=\sqrt{\Bigg{(}\sum_{j=0}^{p}[\lambda(x)]_{j}\Bigg{)}^{2}+\sum_{j=1}^{p}[\lambda(x)]_{j}^{2}}$
(20)
$\displaystyle\leq\sqrt{\Bigg{(}\sum_{j=0}^{p}\left\lvert[\lambda(x)]_{j}\right\lvert\Bigg{)}^{2}+\sum_{j=1}^{p}\left\lvert[\lambda(x)]_{j}\right\lvert^{2}}$
$\displaystyle\leq\sqrt{\big{(}(p+1)\Lambda\big{)}^{2}+p\Lambda^{2}}<\sqrt{2(p+1)^{2}\Lambda^{2}}$
$\displaystyle=\Lambda\sqrt{2}(p+1).$

Let $\overline{v}\in\mathbb{R}^{n+1}$ be a singular right unit vector of
$\mathbf{Q}^{\dagger^{T}}$ associated with the largest singular value
$\sigma_{1}$, and $\overline{x}\in\mathbb{R}^{n}$ the maximizer of
$\left\lvert\overline{v}^{T}\overline{\phi}_{L}(x)\right\lvert$ in
$\overline{B}(0,1)$. Using the SVD decomposition and the fact that
$\left\lVert\mathbf{Q}^{\dagger}\right\lVert=\sigma_{1}$ it is not hard to
show (see [5, Lemma 3.13]) that

$\left\lVert\mathbf{Q}^{\dagger^{T}}\overline{\phi}_{L}(\overline{x})\right\lVert\geq|\overline{v}^{T}\overline{\phi}_{L}(\overline{x})|\left\lVert\mathbf{Q}^{\dagger^{T}}\right\lVert.$
(21)

Then, from Lemma 3.6, Remark 3.7, (20), and (21),

$\displaystyle\Lambda\sqrt{2}(p+1)$
$\displaystyle\geq\left\lVert\mathbf{Q}^{\dagger^{T}}\overline{\phi}_{L}(\overline{x})\right\lVert\geq\left\lvert\overline{v}^{T}\overline{\phi}_{L}(\overline{x})\right\lvert\left\lVert\mathbf{Q}^{\dagger^{T}}\right\lVert$
$\displaystyle=\left\lvert\overline{v}^{T}\overline{\phi}_{L}(\overline{x})\right\lvert\left\lVert\mathbf{Q}^{\dagger}\right\lVert=\max_{x\in\overline{B}(0,1)}\left\lvert\overline{v}^{T}\overline{\phi}_{L}(x)\right\lvert\left\lVert\mathbf{Q}^{\dagger}\right\lVert$
$\displaystyle\geq\frac{1}{\sqrt{n+1}}\left\lVert\mathbf{Q}^{\dagger}\right\lVert.$

Thus, from the previous inequality and from (18), it follows that

$\left\lVert\hat{\mathbf{L}}_{s}^{\dagger}\right\lVert\leq\left\lVert\mathbf{Q}^{\dagger}\right\lVert\leq\Lambda\sqrt{2(n+1)}(p+1),$

and we complete the proof. ∎

##  4 Conclusions

In this article, we organized several results from the literature about error
bounds for linear and quadratic models related to derivative-free trust-region
algorithms. We also extended the results relaxed results of [9] to
underdetermined models and provided a clearer proof than [5] for the bound on
$\left\lVert\mathbf{L_{s}}^{\dagger}\right\lVert$ in the minimum Frobenius
norm case. Table 1 provides a compilation of the results presented here and is
useful for future works in worst case complexity.

Future work may include the estimation of bounds for the Hessian of minimum
norm underdetermined quadratic models. In [4], the authors obtain bounds for
the projection of errors onto a specific linear subspace, which is not very
useful in practical terms. On the other hand, [5, p. 79], suggests that by
using an overall poisedness constant for the sample set, it is possible to
establish bounds for the model’s Hessian. Another interesting question is
whether models constructed by support vector regression [9] using
underdetermined quadratic polynomials are able to satisfy Assumption 2.4 and
their practical benefits under noisy blackbox functions.

Model Type | Error | Error bound | Ref.
---|---|---|---
Linear determined | $\left\lvert\mathrm{m}(x)-f(x)\right\lvert$ | $\left(\frac{1}{2}L+\kappa+\left(\frac{1}{2}L+2\kappa\right)\Lambda n\right)\delta^{2}$ | [9]
| $\left\lVert\nabla\mathrm{m}(x)-\nabla f(x)\right\lVert$ | $\left(L+\left(\frac{1}{2}L+2\kappa\right)\Lambda n\right)\delta$ | [9]
| $\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert$ | $0$ | [9]
Quadratic determined | $\left\lvert\mathrm{m}(x)-f(x)\right\lvert$ | $\left(\frac{1}{2}L+\kappa+4\Lambda\sqrt{q(q+1)^{3}}\left(2+3\sqrt{2}\right)(\kappa+L)\right)\delta^{2}$ | [9]
| $\left\lVert\nabla\mathrm{m}(x)-\nabla f(x)\right\lVert$ | $\left(8\Lambda\sqrt{q(q+1)^{3}}\left(1+\sqrt{2}\right)(\kappa+L)\right)\delta$ | [9]
| $\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert$ | $8\Lambda\sqrt{2q(q+1)^{3}}(\kappa+L)$ | [9]
Quadratic underdetermined | $\left\lvert\mathrm{m}(x)-f(x)\right\lvert$ | $\left(\frac{1}{2}\left(L+\kappa_{H}\right)+\kappa+2\kappa_{s}\sqrt{p}\left(L+\kappa+\frac{3}{4}\kappa_{H}\right)\right)\delta^{2}$ | This work
| $\left\lVert\nabla\mathrm{m}(x)-\nabla f(x)\right\lVert$ | $2\kappa_{s}\sqrt{p}\left(L+\kappa+\frac{3}{4}\kappa_{H}\right)\delta$ | This work
Min. Frobenius norm | $\left\lvert\mathrm{m}(x)-f(x)\right\lvert$ | $\left(\frac{1}{2}\left(L+\left(\kappa+\frac{L}{2}\right)\Lambda\frac{4(p+1)\sqrt{2(q+1)}}{c(\delta_{\max})^{2}}\right)+\kappa+\right.$ | This work
|  | $\quad+\left.2\Lambda\sqrt{2p(n+1)}(p+1)\left(L+\kappa+\left(\kappa+\frac{L}{2}\right)\Lambda\frac{3(p+1)\sqrt{2(q+1)}}{c(\delta_{\max})^{2}}\right)\right)\delta^{2}$ | This work
| $\left\lVert\nabla\mathrm{m}(x)-\nabla f(x)\right\lVert$ | $2\Lambda\sqrt{2p(n+1)}(p+1)\left(L+\kappa+\left(\kappa+\frac{L}{2}\right)\Lambda\frac{3(p+1)\sqrt{2(q+1)}}{c(\delta_{\max})^{2}}\right)\delta$ | This work
| $\left\lVert\nabla^{2}\mathrm{m}(x)\right\lVert$ | $\left(\kappa+\frac{L}{2}\right)\Lambda\frac{4(p+1)\sqrt{2(q+1)}}{c(\delta_{\max})^{2}}$ | This work
Table 1: Error bounds for linear and quadratic interpolation models under
Assumption 2.4. In the table, $L$ is the Lipschitz constant associated with
$\nabla f$, $\kappa$ is the relaxed interpolation condition, $\Lambda$ is
associated with the geometry of the sample set, $n$ is the dimension of the
domain, $q+1$ is the dimension of $\mathcal{P}_{n}^{a}(\mathbb{R})$, $p+1$ is
the number of sample points, and $\delta\leq\delta_{\max}$ are associated with
the neighborhood of sampling points.

## References

  * [1] E. Birgin, N. Krejić, and J. Martínez.  Inexact restoration for derivative-free expensive function minimization and applications.  J. Comp. Appl. Math., 410:114193, 2022.
  * [2] C. Cartis and L. Roberts.  A derivative-free Gauss-Newton method.  Math. Prog. Comp., 11(4):631–674, 2019.
  * [3] C. Cartis and L. Roberts.  Scalable subspace methods for derivative-free nonlinear least-squares optimization.  arXiv:2102.12016, 2021.
  * [4] A. R. Conn, K. Scheinberg, and L. N. Vicente.  Geometry of sample sets in derivative-free optimization: polynomial regression and underdetermined interpolation.  IMA J. Numer. Anal., 28(4):721–748, 2008.
  * [5] A. R. Conn, K. Scheinberg, and L. N. Vicente.  Introduction to derivative-free optimization.  SIAM, Philadelphia, 2009.
  * [6] R. Garmanjani, D. Júdice, and L. N. Vicente.  Trust-region methods without using derivatives: worst case complexity and the nonsmooth case.  SIAM J. Optim., 26(4):1987–2011, 2016.
  * [7] M. J. D. Powell.  The NEWUOA software for unconstrained optimization without derivatives.  In Large-scale nonlinear optimization, pages 255–297. Springer, Boston, 2006.
  * [8] M. J. D. Powell.  On fast trust region methods for quadratic models with linear constraints.  Math. Prog. Comp., 7(3):237–267, 2015.
  * [9] A. Verdério, E. W. Karas, L. G. Pedroso, and K. Scheinberg.  On the construction of quadratic models for derivative-free trust-region algorithms.  EURO J. Comput. Optim., 5(4):501–527, 2017.
  * [10] S. M. Wild.  Derivative-free optimization algorithms for computationally expensive functions.  PhD thesis, Cornell University, Ithaca, 2008.