# Gemini Dark Matter

Andrew Cheek<EMAIL_ADDRESS>Tsung-Dao Lee Institute & School of Physics
and Astronomy, Shanghai Jiao Tong University, China  Key Laboratory for
Particle Astrophysics and Cosmology (MOE) & Shanghai Key Laboratory for
Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai
200240, China  Yu-Cheng Qiu<EMAIL_ADDRESS>Tsung-Dao Lee Institute &
School of Physics and Astronomy, Shanghai Jiao Tong University, China  Key
Laboratory for Particle Astrophysics and Cosmology (MOE) & Shanghai Key
Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong University,
Shanghai 200240, China  Liang Tan<EMAIL_ADDRESS>Tsung-Dao Lee
Institute & School of Physics and Astronomy, Shanghai Jiao Tong University,
China  Key Laboratory for Particle Astrophysics and Cosmology (MOE) & Shanghai
Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong
University, Shanghai 200240, China

###### Abstract

The $S_{8}/\sigma_{8}$ tension in the large scale structure can be explained
by decaying dark matter with an almost degenerate spectrum and small enough
decay width. Here we propose the Gemini dark matter model, which contains a
heavy mother particle $\chi_{3}$ and two twins $\chi_{1/2}$ which are almost
degenerate in mass and are produced at the same time. The dark sector is
charged under the same Froggatt-Nielsen symmetry that can explain the
hierarchy of the Standard model Yukawa couplings. The slightly heavier
$\chi_{2}$ decays into $\chi_{1}$ and the axionic component of the flavon,
which washes out the small scale structure and resolves $S_{8}/\sigma_{8}$
tension. We present the production mechanism of Gemini dark matter and viable
parameter regions. We find that despite the preferred dark matter mass being
$\mathcal{O}(1)$–$\mathcal{O}(100)$ keV, they constitute cold dark matter. The
Gemini dark matter model predicts an abundance of dark radiation that will be
probed in future measurements of the CMB.

##  I Introduction

The standard model of modern cosmology, $\Lambda$CDM model, currently provides
a consistent picture of the observable Universe to a satisfying level of
precision [1]. In recent years however, some observations have put tension on
the $\Lambda$CDM [2] such as the Hubble tension and the $S_{8}/\sigma_{8}$
tension [3]. If either of these tensions persist, cosmologists will have to go
beyond $\Lambda$CDM. This paper focuses on a new physics solution to the
$S_{8}/\sigma_{8}$ tension.

Decaying dark matter (DDM) [4] has been proposed to resolve the
$S_{8}/\sigma_{8}$ tension [5, 6, 7, 8, 9, 10, 11, 12]. The basic solution
states that a cold dark matter (CDM) decays into a slightly lighter dark
particle. The mass difference gives the product particle some velocity so it
acts as warm dark matter (WDM), washing out structure a little. The decay
lifetime should be $\tau\sim\mathcal{O}(10)$–$\mathcal{O}(100)\,{\rm Gyr}$ and
the mass ratio between the decay product and the CDM quantified as
$\epsilon\equiv(m_{\rm CDM}^{2}-m_{\rm WDM}^{2})/2m_{\rm CDM}^{2}$, must be
$\sim 0.01$–$0.1$ [12]. This almost degenerate dark sector and highly
suppressed decay must also evade indirect detection bounds. This could be
achieved by assigning new quantum numbers and discrete symmetries.

The standard model (SM) of particle physics also contains several mysteries,
and one of them is the mass hierarchy between three generations of quarks and
leptons. One solution is the well-developed Froggatt-Nielsen (FN) mechanism
[13, 14, 15], which introduces a chiral $U(1)_{\rm FN}$ (global or gauge)
symmetry to the SM particles. This new symmetry forbids Yukawa terms at the
renormalizable level. The dimension at which the Yukawa term appears at an
effective nonrenormalizable level depends on the FN quantum numbers of the
given particle. To retrieve the SM description, the $U(1)_{\rm FN}$ symmetry
is spontaneously broken where the Yukawa couplings are now proportional to a
common parameter, $\lambda<1$ to some positive generation-dependent power, and
hierarchy emerges. This new symmetry is associated with a new scalar field,
named flavon $\Phi$, and it mediates flavor-changing current, whose coupling
strengths with fermions are proportional to the respective fermion masses. In
other words, flavon couplings are generation-dependent.

The FN framework can be easily extended to the dark sector [16]. In the
minimal setup, the dark sector interacts with the SM particles by mediating
the flavon $\Phi$. Cosmologically, dark matter can be produced through either
thermal freeze-out or non-thermal freeze-in [16, 17, 18, 19]. Since the
coupling between the flavon and the dark matter particle depends on dark
matter mass, the direct and indirect detection may have interesting
signatures.

Reference [17] considered DDM under the FN framework and suggested a potential
model to resolve the $S_{8}$, albeit very roughly. In this paper, we
systematically study the mass spectrum required to solve the $S_{8}$ tension.
We have identified a more consistent solution that results in different
cosmology and phenomenological signals. In the dark sector, three dark
fermions are charged under the FN symmetry. Two of them $\chi_{1/2}$ are
almost degenerate and the other one $\chi_{3}$ is much heavier, $m_{3}\gg
m_{1/2}$. The decay channel $\chi_{2}\to\chi_{1}+a$, where $a$ is the axionic
component of the complex FN scalar $\Phi$, is responsible for resolving the
$S_{8}/\sigma_{8}$. The $\chi_{1/2}$ are produced from $\chi_{3}$ decay,
$\chi_{3}\to\chi_{1/2}+a$. We name this model ‘Gemini dark matter’ (Gemini DM
hereafter), where particles $\chi_{1/2}$ are the twins and $\chi_{3}$ is the
mother particle. The mother particle is produced through freeze-in from the SM
thermal bath. The preferred twin masses in this model are roughly
$\mathcal{O}(1)$–$\mathcal{O}(100)$ keV, and the average velocity of them is
$\sim 10^{-6}$, which is cold despite its small mass.

The Gemini DM model predicts the existence of the dark radiation relic. During
the parturition $\chi_{3}\to\chi_{1/2}+a$, the light axionic flavon $a$ is
also produced as dark radiation, which shall contribute to the deviation from
neutrino effective number of species $\Delta N_{\rm eff}$. It may be probed in
future observations by CMB-S4 [20] and CMB-HD [21].

The paper goes as follows. Sec. II constructs the Gemini DM model. The
cosmological production of Gemini DM is explained in Sec. III. Summary and
discussions are in Sec. IV. Appendix A reviews the Froggatt-Nielsen mechanism
in the SM sector.

##  II Gemini dark matter

In this section, we construct the Gemini DM model by extending the Froggatt-
Nielsen symmetry to the dark sector, which contains Weyl fermions $\chi_{i}$
carries FN charge $n_{i}$, which is also Standard Model gauge group singlet.
There are three generations of $\chi$s. We focus on constructing a particle
spectrum that resolves the $S_{8}/\sigma_{8}$ tension, i.e. we require two
particles with almost degenerate masses.

With our setup, below some UV cutoff scale $\Lambda$, the leading order
nonrenormalizable operators that respect the global $U(1)_{\rm FN}$ appearing
in the dark sector Lagrangian density are

$\mathcal{L}\supset
i\bar{\chi}_{j}\bar{\sigma}^{\mu}\partial_{\mu}\chi_{j}-\frac{\beta_{ij}}{2}\frac{\Phi^{n_{i}+n_{j}}}{\Lambda^{n_{i}+n_{j}-1}}\chi_{i}\chi_{j}+{\rm
h.c.}\;,$ (1)

where $\beta_{ij}\sim\mathcal{O}(1)$ are the coupling constants and $\Phi$ is
the flavon field is with FN charge $-1$. Here bar indicates the conjugate and
$\bar{\sigma}^{\mu}=(1,-\vec{\sigma})$, where $\vec{\sigma}$ is an array of
Pauli matrices. The summation of dummy indices is employed here. The complex
FN scalar $\Phi$ can acquire a nontrivial vacuum expectation value (VEV),
$\langle\Phi\rangle=f_{a}$, under its potential, leading to spontaneous
symmetry breaking of $U(1)_{\rm FN}$. Since we are considering the global
$U(1)_{\rm FN}$, and quantum gravity cannot accommodate such symmetries, it
must be broken at least by the Planck scale [22, 23]. Therefore, the goldstone
boson $a$ (we call it axion for familiarity) acquires a mass through an
explicit breaking term. After spontaneous symmetry breaking, the field can be
written as $\Phi=e^{ia/f_{a}}(f_{a}+\phi)/\sqrt{2}$. A Taylor expansion gives

$\left(\frac{\Phi}{\Lambda}\right)^{n}\to\lambda^{n}e^{ina/f_{a}}\left(1+\frac{n\phi}{f_{a}}+\cdots\right)\;,$
(2)

where the FN parameter is defined as $\lambda\equiv f_{a}/\sqrt{2}\Lambda$.
The leading interaction terms in Eq.(1) are the Yukawa-type interactions that
determine the $\chi$ masses. After performing the phase rotation for each
fermion, $\chi_{j}\to e^{-in_{j}a/f_{a}}\chi_{j}$, one rotates the axion away.
This phase rotation generates couplings between $\partial a$ and the $\chi$’s
through the $\Phi$ fields kinetic term. Using the equation of motion for the
$\chi$’s, one obtains the interaction between $a$ and fermions, due to the
mismatch between mass eigenstates and coupling states, flavour changing
currents that couple to flavon are still present. After diagonalizing the mass
matrix, interactions in Eq. (1) are expressed as

$-\mathcal{L}\supset\frac{1}{2}m_{k}\chi_{k}\chi_{k}+g_{ij}^{\phi}\phi\chi_{i}\chi_{j}+g_{ij}^{a}a\chi_{i}\chi_{j}+{\rm
h.c.}\;,$ (3)

where $m_{k}$ is the $k$-th diagonal component of $D$ and

$\displaystyle D$ $\displaystyle={\rm diag}(m_{1},m_{2},m_{3})=U^{\top}MU$
$\displaystyle g_{ij}^{\phi}$ $\displaystyle=\frac{1}{2}\left[{\rm
sym}\left(U^{\top}\frac{1}{f_{a}}\frac{\partial(\lambda
M)}{\partial\lambda}U\right)\right]_{ij}$ $\displaystyle g_{ij}^{a}$
$\displaystyle=\left[{\rm sym}\left(\frac{ND}{f_{a}}\right)\right]_{ij}$
$\displaystyle M_{ij}$
$\displaystyle=\frac{1}{\sqrt{2}}f_{a}\beta_{ij}\lambda^{n_{i}+n_{j}-1}$
$\displaystyle N_{ij}$ $\displaystyle=(U^{\dagger})_{ik}n_{k}U_{kj}\;.$

The $U$ is a unitary matrix, $UU^{\dagger}=1$. Here $[{\rm sym}(A)]_{ij}\equiv
A_{ij}+A_{ji}-A_{ii}\delta_{ij}$, for any matrix $A$. The presence of this
symmetric sum for the off-diagonal couplings in $g_{ij}^{\phi,a}$ is due to
the spinor identity $\chi_{i}\chi_{j}=\chi_{j}\chi_{i}$. The nonvanishing
flavor-changing couplings between flavons and $\chi$ fields provide decay
channels.

The Gemini DM model requires at least three generations of $\chi$. The reason
is the following. In the minimal setup, we require the decay channel
$\chi_{2}\to\chi_{1}+a$ to resolve $S_{8}$ tension. Meanwhile, the mass split
between $\chi_{1}$ and $\chi_{2}$ is small, which is described by

$\epsilon=\frac{1}{2}\left(1-\frac{m_{1}^{2}}{m_{2}^{2}}\right)\in(0.01,0.1)\;.$
(4)

Recall that the mass matrix under the FN framework is approximately rank one.
If only two generations are present, then one must have $m_{2}\gg m_{1}$,
which is not desirable. With three generations, one could have $m_{3}\gg
m_{2}\approx m_{1}$ under a reasonable choice of $\beta_{ij}$. In the Gemini
DM model, we call $\chi_{3}$ the mother particle and $\chi_{1/2}$ the twins.

Here we present a specific parameter choice that gives us a desirable
spectrum. We fix the FN parameter as $\lambda=0.171$ [24] as a benchmark
model. The value of FN parameter $\lambda$ has degeneracy with FN charge
choices [25]. Choose the FN charge as

$n_{1}=4.5+n_{3}\;,\quad n_{2}=2.5+n_{3}\;.$ (5)

The charge $n_{3}$ only determines the overall scale, which is irrelevant when
determining couplings and mass ratios. The half charge can be absorbed by
redefining the FN parameter. We parametrize $\beta$ as

$\beta=\begin{pmatrix}1&1&1+c\\\ 1&1&1\\\ 1+c&1&1\end{pmatrix}\;.$ (6)

in order to obtain various spectra by varying the parameter $c$ and determine
how to produce the Gemini DM scenario. As shown in the upper panel of Fig. 1,
for $-1.7\lesssim c\lesssim 2.0$, one achieves the Gemini DM spectrum with an
$\epsilon$ that resolves the $S_{8}$ tension using Eq. (3). Furthermore, the
corresponding mass ratio of the much heavier mother particle $\chi_{3}$ to the
Gemini is shown in lower panel of Fig. 1.

Figure 1: The upper panel is the mass difference $\epsilon$ between the twins
$\chi_{1/2}$ plotted against the single parameter $c$ describes the
$\beta_{ij}$ in Eq. (6). The lower panel shows the mass ratio between the
mother particle and the twins as a function of $c$.

We parametrize the flavon couplings in Eq. (3) as

$\displaystyle g_{ij}^{\phi}$
$\displaystyle=\frac{1}{f_{a}}\left(m_{i}-m_{j}+m_{i}\delta_{ij}\right)\mathcal{A}_{ij}$
$\displaystyle g_{ij}^{a}$
$\displaystyle=\frac{1}{f_{a}}\left(m_{i}-m_{j}+m_{i}\delta_{ij}\right)\mathcal{B}_{ij}\;,$
(7)

where $\mathcal{A}$ and $\mathcal{B}$ are numerical matrices that are
calculated explicitly with our choice of $\beta$ (6) and FN charges (5). By
varying $c$, one obtains a one-to-one correspondence between couplings and
$\epsilon$. To determine the values and variance of the matrices $\mathcal{A}$
and $\mathcal{B}$ as shown in Table 1, we perform the following method. We
randomly sample $c$, adopting the uniform distribution as the prior, and
select out those couplings associated with $\epsilon\in(0.01,0.1)$. From these
values we collect the corresponding $\mathcal{A}$ and $\mathcal{B}$. We
calculate the average and their deviations as shown in Table 1. The result
here is that the ranges of values in the matrices are not large compared to
their central value. We can therefore justifiably treat them as constants for
the following calculations.

Table 1: Statistics of matrix elements $|\mathcal{A}_{ij}|^{2}$ and $|\mathcal{B}_{ij}|^{2}$ under the parametrization Eq. (5), (6) and Eq. (7) with $\epsilon\in(0.01,0.1)$, and randomly sampled $c$ under uniform distribution. Central values are averages. The upper and the lower uncertainties indicate the maximum and the minimum. $(ij)$ | $|\mathcal{A}_{ij}|^{2}$ | $|\mathcal{B}_{ij}|^{2}$
---|---|---
$(11)$ | $6.64^{+0.49}_{-0.47}\times 10^{3}$ | $48.6^{+0.4}_{-0.4}$
$(22)$ | $6.76^{+0.48}_{-0.48}\times 10^{3}$ | $49.4^{+0.3}_{-0.4}$
$(33)$ | $1.57\times 10^{3}\pm 10^{-3}$ | $12.3\pm 10^{-5}$
$(21)$ | $1.37\times 10^{2}\pm 10^{-1}$ | $0.999^{+0.001}_{-0.002}$
$(31)$ | $5.77^{+0.95}_{-0.88}\times 10^{-2}$ | $4.21^{+0.70}_{-0.65}\times 10^{-4}$
$(32)$ | $6.79^{+1.11}_{-0.99}\times 10^{-2}$ | $4.98^{+0.79}_{-0.75}\times 10^{-4}$

Since $\phi$ obtains its mass from the FN symmetry breaking, $m_{\phi}$ is
typically associated with the breaking scale $f_{a}$. The axion, $a$, is
typically very light, obtaining its mass from some explicit breaking, which we
treat as a free parameter and $m_{a}\ll m_{1/2}$. In line with this thinking,
we consider the case that $\chi$ dominantly decays through $a$. With
$\chi_{2}$ slightly heavier than $\chi_{1}$ ($m_{2}\gtrsim m_{1}$), $\chi_{3}$
decays via $\chi_{1/2}+a$ while $\chi_{2}$ dominantly decays to $\chi_{1}+a$.
Heavy $\phi$ decays into $\chi$s and $a$. These decay widths are given by

$\displaystyle\Gamma_{\chi_{i}\to\chi_{j}a}$
$\displaystyle=\frac{m_{i}|g_{ij}^{a}|^{2}}{16\pi}\gamma_{+}\left(\frac{m_{j}}{m_{i}},\frac{m_{a}}{m_{i}}\right)\;,$
$\displaystyle\Gamma_{\phi\to\chi_{k}\bar{\chi}_{l}}$
$\displaystyle=\frac{m_{\phi}|g_{ij}^{\phi}|^{2}}{8\pi(1+\delta_{kl})}\gamma_{-}\left(\frac{m_{k}}{m_{\phi}},\frac{m_{l}}{m_{\phi}}\right)\;.$
(8)

Kinematics requires that $m_{i}>m_{j}+m_{a}$ and $m_{\phi}>m_{k}+m_{l}$. The
two-body decay phase space factor $\gamma_{\pm}$ is

$\displaystyle\gamma_{\pm}(y,z)$ $\displaystyle=(1+y+z)^{3/2}(1\pm y-z)^{3/2}$
$\displaystyle\qquad\qquad\times\sqrt{(1-y+z)(1\mp y-z)}\;,$

for any real number $y$ and $z$. There are sub-dominant dark sector decays
into SM particles mediated by $\phi$ and $a$, $\chi_{i}\to\chi_{j}+{\rm
SM}+{\rm SM}$. However, it is a three-body decay, suppressed by at least
another coupling squared $\sim(m_{\rm SM}/f_{a})^{2}$. So the three-body decay
branching ratio (Br) is negligible.

Our twins are dark matter and the decay $\chi_{2}\to\chi_{1}+a$ would give the
lighter $\chi_{1}$ enough kinetic energy to wash out some structures, which
results in smaller $S_{8}$ for the late universe. As long as the decay width
is approximately

$\Gamma_{\chi_{2}\to\chi_{1}a}=\frac{1}{\tau_{8}}\sim\mathcal{O}(10^{-43})\,{\rm
GeV}\;,$ (9)

this decay resolves $S_{8}/\sigma_{8}$ tension, equivalently,
$\tau_{8}\sim\mathcal{O}(10)$–$\mathcal{O}(100)$ Gyr [12]. It is projected
into a region on the $f_{a}$-$m_{1/2}$ parameter space111Note that for
$\tau>100\,{\rm Gyr}$, the model still works as a dark matter model without
resolving the $S_{8}$. However, for $\tau<10\,{\rm Gyr}$, it has conflicts
with other structure formation [12] constraints.. Choosing $\epsilon=0.05$ and
$m_{a}=10^{-6}\,{\rm eV}$ as a benchmark model, the preferred parameter region
could be parametrized as

$m_{1/2}\approx 37\times\left(\frac{f_{a}^{2}}{\tau_{8}}\right)^{1/3}\;,$ (10)

which is the gray stripe shown in Fig. 2.

The axion mass $m_{a}$ puts an upper bound on $f_{a}$ via its overproduction
from the misalignment mechanism. For $m_{a}\sim 10^{-6}$ eV, this gives
$f_{a}<1.1\times 10^{13}$ GeV [26]. There is a lower bound on $f_{a}$ in this
figure. Our axion $a$ couples to nucleons with a coupling of order $\sim
1/f_{a}$. This means that supernovas can emit a large axion flux, which is
constrained as $f_{a}>10^{9}$ GeV by the Kamiokande-II neutrino detector at
the time of SN 1987A [27]. Furthermore, the axion will couple to the photon
via SM fermion loops, this coupling is highly suppressed $g_{a\gamma}\sim
m_{a}^{2}/f_{a}^{2}$, avoids constraints by orders of magnitude [28] and
leaves the axion stable on cosmological timescales. Note that with $m_{2}\sim
100\,{\rm keV}$ the three-body $\chi_{2}$ decay into charged fermions is
kinematically forbidden not just suppressed. However, the
$\chi_{2}\to\chi_{1}+\gamma+\gamma$ decay occurs through a loop but has a
decay width of
$\Gamma_{\chi_{1}\to\chi_{1}\gamma\gamma}\lesssim(m_{2}/f_{a})^{2}\Gamma_{\chi_{2}\to\chi_{1}a}\sim\mathcal{O}(10^{-48})\,{\rm
s}^{-1}$, well below current limits [29].

Figure 2: Constraints on $f_{a}$ and $m_{2}$ parameter space. Here we choose
the $\epsilon=0.05$ and $m_{a}=10^{-6}\,{\rm eV}$ as the benchmark model. The
gray stripe is the region that satisfies Eq. (9) and resolves $S_{8}$ tension.
The red region is excluded by the supernova axions [27].

The Gemini DM model has several parameters, among which $\epsilon$ fixes the
mass ratio between the twins, and $\tau_{8}$ restricts $m_{1/2}$ in terms of
the breaking scale $f_{a}$. For now, we still have the mass of the mother
particle $m_{3}$ and the scale $f_{a}$ undetermined. As long as $m_{a}$ is
small, it does not affect Gemini DM phenomenology. In the following section,
we show that the production of Gemini DM in the early Universe will allow us
to limit the $f_{a}$ and $m_{3}$ values that solve the $S_{8}$ tension. This
means that there is essentially only one free parameter, which is the $f_{a}$.

##  III Cosmological production and constraints

We propose a consistent production mechanism for the Gemini DM. The mother
particle $\chi_{3}$ freezes in through thermal $\phi$ decay, and subsequently
gives birth to the twins $\chi_{1/2}$ via $\chi_{3}\to\chi_{1/2}+a$. This
mechanism could be divided into three stages, which are (1) $\phi$ stays in
thermal; (2) $\phi$ decays into $\chi_{3}$; and (3) $\chi_{3}$ decays into
$\chi_{1}$ or $\chi_{2}$.

###  III.1 $\phi$ in thermal equilibrium with the Standard Model bath

From Fig. 2 one sees that current experiments require that at least
$f_{a}>10^{9}\,{\rm GeV}$. Meanwhile, the $m_{\phi}$ is generated from the
potential that breaks the $U(1)_{\rm FN}$. This means that $m_{\phi}$ is much
larger than the Higgs mass scale $\sim\mathcal{O}(100)\,{\rm GeV}$.
Furthermore, the SM thermal bath temperature required to produce appreciable
quantities of $\phi$ is much higher than the electroweak breaking scale. As
reviewed in appendix A, at this scale, both $\phi$ and $a$ couple to SM
particles through

$n\lambda^{n}\frac{\phi}{f_{a}}\overline{Q}\mathcal{H}q\;,\quad
n\lambda^{n}\frac{a}{f_{a}}\overline{Q}\mathcal{H}q\;,$ (11)

where $n=n_{Q}^{\rm FN}+n_{q}^{\rm FN}$ is the sum of the quark FN charges.
$\mathcal{H}$ is the SM Higgs doublet that does not carry an FN charge. We
shall see that flavons enter the SM thermal bath through these interactions.
Since we are interested in the temperature, $T\sim\mathcal{O}(m_{\phi})$, we
treat all SM particles as massless for simplicity. Focus on $\phi$ for the
moment. There are three different processes associated with this vertex, which
are

$\displaystyle\phi+\mathcal{H}$ $\displaystyle\to\bar{Q}+q\;,$
$\displaystyle\phi+Q$ $\displaystyle\to\mathcal{H}+q\;,$ $\displaystyle\phi+q$
$\displaystyle\to\mathcal{H}+Q\;.$ (12)

In the leading order, they are described as having the same amplitude, which
is given by

$|\mathcal{M}|^{2}\approx\frac{\lambda^{2}}{f_{a}^{2}}(2p_{Q}p_{q})\;,$ (13)

where $p_{j}$ is the four-momentum of $j$ particle. Here we have explicitly
chosen $n=1$ because this will be the leading contribution. The other three
particles, namely $\\{\mathcal{H},Q,q\\}$, are in thermal equilibrium. The
equilibrium number density $n_{j}^{\rm eq}$ can be approximated by integrating
the Maxwell-Boltzmann distribution $f_{j}^{\rm eq}\approx e^{-E_{j}/T}$,

$\displaystyle n_{j}^{\rm eq}$ $\displaystyle=\frac{g_{j}}{(2\pi)^{3}}\int
d^{3}p_{j}f_{j}^{\rm eq}$ (14)
$\displaystyle=\frac{g_{j}}{2\pi^{2}}m_{j}^{2}TK_{2}(m_{j}/T)\overset{m_{j}\to
0}{=}\frac{g_{j}}{\pi^{2}}T^{3}\;.$

Here $g_{j}$ is the degree of freedom of $j$ particle and $m_{j}$ is its mass.
$T$ labels the temperature. The last approximation above takes the massless
limit, which applies for $\\{\mathcal{H},Q,q\\}$. $K_{\alpha}(x)$ is the
modified Bessel function of the second kind. The number density $n_{\phi}$
evolution of $\phi$ is governed by the Boltzmann equation, which includes
three processes in Eq. (12) and their inverses,

$\displaystyle\dot{n}_{\phi}+3Hn_{\phi}$
$\displaystyle=-\tilde{\Gamma}_{\phi}\left(n_{\phi}-n_{\phi}^{\rm eq}\right)$
(15) $\displaystyle\tilde{\Gamma}_{\phi}$
$\displaystyle=\langle\sigma_{\mathcal{H}}v\rangle n_{\mathcal{H}}^{\rm
eq}+\langle\sigma_{Q}v\rangle n_{Q}^{\rm eq}+\langle\sigma_{q}v\rangle
n_{q}^{\rm eq}\;,$

where $\langle\sigma_{j}v\rangle$ labels the thermally averaged cross section
of the process whose initial states are $\phi$ and $j$ particles. The $H(T)$
is the Hubble parameter that describes the expansion of the universe. The
interaction rate $\tilde{\Gamma}_{\phi}(T)$ and the Hubble expansion rate
$3H(T)$ determines whether the $\phi$ stays in the thermal equilibrium. Since
at the relevant temperatures, the particles $\\{\mathcal{H},Q,q\\}$ are all
massless, one can simplify the expression of $\tilde{\Gamma}_{\phi}$ by
performing a relabeling of the integration variables in the second and third
terms in $\tilde{\Gamma}_{\phi}$. This allows us to write

$\displaystyle\tilde{\Gamma}_{\phi}$
$\displaystyle\approx\frac{1}{n_{\phi}^{\rm
eq}}\int\prod_{j}d\Pi_{j}\frac{\lambda^{2}}{f_{a}^{2}}\left(2p_{Q}p_{q}+2p_{\mathcal{H}}p_{q}+2p_{Q}p_{\mathcal{H}}\right)$
(16)
$\displaystyle\qquad\times(2\pi)^{4}\delta^{4}(p_{\phi}+p_{\mathcal{H}}-p_{Q}-p_{q})e^{-(E_{\phi}+E_{\mathcal{H}})/T}\;,$

where the phase space integral is defined as
$d\Pi_{j}=g_{j}d^{3}p_{j}/[(2\pi)^{3}(2E_{j})]$, and $j$ runs for all four
particles in these processes. Following the method in Ref. [30], one converts
the above formula into a single integration,

$\displaystyle\tilde{\Gamma}_{\phi}$
$\displaystyle\approx\frac{g_{\mathcal{H}}g_{Q}g_{q}\lambda^{2}}{16(2\pi)^{3}}\frac{T^{5}}{f_{a}^{2}m_{\phi}^{2}}\frac{\mathcal{I}(m_{\phi}/T)}{K_{2}(m_{\phi}/T)}\;,$
(17) $\displaystyle\mathcal{I}\left(\zeta\right)$
$\displaystyle=\int_{\zeta}^{\infty}d\xi\left(\xi^{2}-\zeta^{2}\right)\left(2\xi^{2}-\zeta^{2}\right)K_{1}(\xi)\;,$

where $\zeta$ is a real number. Taking the standard cosmology picture we will
assume that the Universe is dominated by radiation at early times. Then the
Hubble parameter can be approximated by

$H(T)\approx\frac{\pi}{3}\sqrt{\frac{g_{\star}}{10}}\frac{T^{2}}{M_{\rm
Pl}}\approx 3.4\times\frac{T^{2}}{M_{\rm Pl}}\;,$ (18)

where $M_{\rm Pl}=2.4\times 10^{18}\,{\rm GeV}$ is the reduced Planck scale
and $g_{\star}$ is the temperature-dependent effective number of energy
degrees of freedom. To get the numerical expression we take
$g_{\star}(T\gtrsim T_{\rm EWSB})\approx 107$ [31].

To estimate when $\phi$ decouples form the SM bath we compare
$\tilde{\Gamma}_{\phi}(T)$ and $3H(T)$. Similarly, this can be done for the
$a$ particle, and the interaction rate of (17) is modified by taking the
massless limit

$\tilde{\Gamma}_{a}\approx\frac{g_{\mathcal{H}}g_{Q}g_{q}\lambda^{2}}{(2\pi)^{3}}\frac{T^{3}}{f_{a}^{2}}\;.$
(19)

In Fig. 3 we show a specific realization where $\tilde{\Gamma}_{\phi}>3H$ for
temperatures below the breaking scale $f_{a}$. We see that the axion
interaction rate $\tilde{\Gamma}_{a}$ follows the simple $T^{3}$ dependence
which, decreases in a cooling Universe faster than Hubble, this leads to an
axion decoupling at around $\sim 10^{5}\,{\rm GeV}$ (left of the figure) [32].
Interestingly, for massive $\phi$ the temperature dependence changes below
$T\sim m_{\phi}$ keeping it in equilibrium and its number density can be
described by $n_{\phi}^{\rm eq}$.

Figure 3: Here we choose $f_{a}=2\times 10^{10}$ GeV and $m_{\phi}=10^{9}$
GeV. The degrees of freedom are $g_{Q}=12$, $g_{q}=6$ and $g_{\mathcal{H}}=4$.
The FN charge is chosen as $n=1$ as the leading contribution, and the FN
parameter is $\lambda=0.171$. Note the red-dashed line, the interaction rate
for the axion $\tilde{\Gamma}_{a}$, cross with the Hubble parameter at $T_{\rm
dec}^{a}\approx 5\times 10^{4}$ GeV under this choice of parameters.

###  III.2 $\chi_{3}$ freeze-in from $\phi$ decay

Given that there are parameter choices for our model that ensure that $\phi$
stays in the thermal bath, we consider the standard freeze-in for $\chi_{3}$
from $\phi$ decay [33, 34, 35, 36, 37, 38]. The Boltzmann equation for
$\chi_{3}$ is

$\dot{n}_{\chi_{3}}+3Hn_{\chi_{3}}=S(\phi\to\chi_{3}\bar{\chi}_{3})\;,$ (20)

where the source term is given by

$\displaystyle S(\phi\to\chi_{3}\bar{\chi}_{3})$
$\displaystyle=g_{\phi}g_{3}^{2}\int\frac{d^{3}p_{\phi}}{(2\pi)^{3}}\frac{m_{\phi}}{E_{\phi}}\Gamma_{\phi\to\chi_{3}\bar{\chi}_{3}}e^{-E_{\phi}/T}$
$\displaystyle=\frac{g_{\phi}g_{3}^{2}}{2\pi^{2}}\Gamma_{\phi\to\chi_{3}\bar{\chi}_{3}}m_{\phi}^{2}TK_{1}(m_{\phi}/T)\;.$
(21)

The decay width $\Gamma_{\phi\to\chi_{3}\bar{\chi}_{3}}$ is given by Eq. (8).
In the second line, we have evaluated the integration explicitly. Following
standard parametrization, we define the yield $Y_{3}=n_{\chi_{3}}/s$, where
$s$ is the entropy density. The equilibrium condition of conservation of total
entropy implies that $ds/dt=-3sH$. So the l.h.s. of the Boltzmann equation can
be written as $\dot{n}_{\chi_{3}}+3Hn_{\chi_{3}}=sdY_{3}/dt$. Furthermore,
defining the variable $x=m_{\phi}/T$ enables one to write
$dx/dt=-(x/T)(dT/dt)=xH$ for the radiation-dominated era. Then the Hubble
parameter can be written as $H(x)\approx 3.4\times m_{\phi}^{2}/M_{\rm
Pl}x^{2}$. Similarly, we extract out $x$-dependence in the entropy density,
$s\approx 27.3\times s_{0}m_{\phi}^{3}/T_{0}^{3}x^{3}$, where $T_{0}=2.3\times
10^{-13}\,{\rm GeV}$ is the CMB radiation temperature today. Therefore, the
Boltzmann equation for $\chi_{3}$ (20) becomes

$\frac{dY_{3}}{dx}\approx\frac{g_{\phi}g_{3}^{2}}{186\pi^{2}}\frac{T_{0}^{3}M_{\rm
Pl}}{s_{0}m_{\phi}^{2}}\Gamma_{\phi\to\chi_{3}\bar{\chi}_{3}}x^{3}K_{1}(x)\;,$
(22)

which can be integrated directly, either numerically or taking
$\int_{0}^{\infty}x^{3}K_{1}(x)=3\pi/2$. If we assume the initial yield of
$\chi_{3}$ is zero, we obtain the freeze-in yield of $\chi_{3}$,

$\displaystyle Y_{3}^{\rm f.i.}$ $\displaystyle\approx
Y_{3}(\infty)\approx\frac{3g_{\phi}g_{3}^{2}}{371\pi}\frac{T_{0}^{3}M_{\rm
Pl}}{s_{0}m_{\phi}^{2}}\Gamma_{\phi\to\chi_{3}\bar{\chi}_{3}}$ (23)
$\displaystyle\approx 0.32\times\frac{T_{0}^{3}M_{\rm
Pl}}{s_{0}m_{\phi}}\frac{m_{3}^{2}}{f_{a}^{2}}\left(1-\frac{4m_{3}^{2}}{m_{\phi}^{2}}\right)^{3/2}\;.$

In the last line, we use $g_{\phi}=1$ and $g_{3}=2$ and take the numerical
value for $|\mathcal{A}_{33}|^{2}$ in the decay width from Table 1.

###  III.3 Relic $\chi_{1/2}$ from $\chi_{3}$ decay

In the Gemini DM model, dark matter consists of two twins, namely $\chi_{1}$
and $\chi_{2}$. They are produced from $\chi_{3}$ decay via
$\chi_{3}\to\chi_{1/2}+a$. One can write the relic yield of $\chi_{i}$ as
$Y_{i}\approx{\rm Br}(\chi_{3}\to\chi_{i}a)Y_{3}^{\rm f.i.}$ 222The branching
ratio can be approximated by using values for $|\mathcal{B}_{ij}|^{2}$ in
Table 1. For example, ${\rm
Br}(\chi_{3}\to\chi_{2}a)\approx|\mathcal{B}_{32}|^{2}/(|\mathcal{B}_{32}|^{2}+|\mathcal{B}_{31}|^{2})$.
and we have $Y_{1}+Y_{2}\approx Y_{3}^{\rm f.i.}$. The energy density today is
given by $\rho_{i}(T_{0})=m_{1/2}Y_{i}s_{0}$, where we have approximated
$m_{1}\approx m_{2}\approx m_{1/2}$. Now the total relative relic density of
the Gemini DM is given by

$\displaystyle\Omega_{\rm DM}h^{2}$
$\displaystyle=\frac{(\rho_{1}+\rho_{2})h^{2}}{3M_{\rm
Pl}^{2}H_{0}^{2}}\approx\frac{m_{1/2}Y_{3}^{\rm f.i.}s_{0}h^{2}}{3M_{\rm
Pl}^{2}H_{0}^{2}}$ (24) $\displaystyle\approx
0.11\times\frac{h^{2}T_{0}^{3}m_{1/2}}{H_{0}^{2}M_{\rm
Pl}m_{\phi}}\frac{m_{3}^{2}}{f_{a}^{2}}\left(1-\frac{4m_{3}^{2}}{m_{\phi}^{2}}\right)^{3/2}\;.$

To resolve the $S_{8}$ tension, the decay width of
$\Gamma_{\chi_{2}\to\chi_{1}a}$ satisfying Eq. (9) puts a relation between
$m_{1/2}$ and $f_{a}$. For the benchmark model with $\epsilon=0.05$ and
$m_{a}=10^{-6}\,{\rm eV}$ (10), this gives

$\displaystyle\Omega_{\rm DM}h^{2}$ $\displaystyle\approx
4.0\times\frac{h^{2}T_{0}^{3}m_{3}^{2}}{H_{0}^{2}M_{\rm
Pl}m_{\phi}(\tau_{8}f_{a}^{4})^{1/3}}\left(1-\frac{4m_{3}^{2}}{m_{\phi}^{2}}\right)^{3/2}$
(25) $\displaystyle\approx 0.12\left(\frac{m_{3}}{1.1\times 10^{4}\,{\rm
GeV}}\right)^{2}\left(\frac{f_{a}}{2\times 10^{10}\,{\rm
GeV}}\right)^{-4/3}\;,$

where the Hubble parameter today is $H_{0}/h\approx 2.1\times 10^{-42}\,{\rm
GeV}$ (from $H_{0}\equiv 100h\,{\rm km/s/Mpc}$) [39]. For the last line, we
have chosen $\tau_{8}=10\,{\rm Gyr}$, and $m_{\phi}=10^{9}\,{\rm GeV}$. To
give the correct relic that fits the cosmology, we should have $\Omega_{\rm
DM}h^{2}\approx 0.12$. Figure 4 shows the allowed parameter space for $f_{a}$
and $m_{3}$ that resolves $S_{8}$ and indicates the correct Gemini DM relic by
the gray shaded strip. For a fixed $m_{\phi}=10^{9}$ GeV, we have a preferred
$f_{a}\sim\mathcal{O}(10^{9})$–$\mathcal{O}(10^{12})$ GeV, which indicates
that the Gemini DM mass should be roughly
$m_{1/2}\sim\mathcal{O}(1)$–$\mathcal{O}(100)$ keV, according to Fig. 2.

Figure 4: Here we choose $m_{\phi}=10^{9}$ GeV and $\Omega_{\rm DM}h^{2}=0.12$
under the benchmark model, $\epsilon=0.05$ and $m_{a}=10^{-6}\,{\rm eV}$. The
gray stripe is the parameter range that gives the correct relic and resolves
$S_{8}$. The orange region is excluded by $T_{\rm dec}^{a}>m_{\phi}/25$, and
the red region is excluded by the supernova axions [27]. The darker purple
region is excluded by overproduction of the dark radiation from $\Delta N_{\rm
eff}<0.276$ [1]. The lighter-purple region indicates $\Delta N_{\rm eff}>0.04$
that can be covered for future CMB-S4 [20], whose sensitivity could reach
$\Delta N_{\rm eff}\sim 0.02$.

In addition to potential probes of this model from structure formation, the
fact that the axion was once in thermal equilibrium means that axion relics
will potentially be observable today. References [40, 41, 42] and others
explore this possibility for the QCD axion. The greatest contribution will be
by altering the number of relativistic degrees of freedom $N_{\rm eff}$ in the
form of dark radiation. Deviations from the SM prediction $(N_{\rm eff})_{\rm
SM}\approx 3.044$ [43, 44, 45, 46, 47, 48] are parameterised by

$\Delta
N_{\text{eff}}=\left.\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{a}}{\rho_{\gamma}}\right|_{T_{\rm
rec}}\;,$ (26)

where $\rho_{\gamma}$ is the photon energy density and $T_{\rm rec}\approx
0.26\,{\rm eV}$ is the recombination temperature [49]. The energy density
ratio is evaluated at the time of the recombination because measurements of
the CMB give the most stringent constraints of $\Delta N_{\rm eff}$. For
example, the Planck collaboration gives $(N_{\rm eff})_{\rm
P18}=2.88^{+0.44}_{-0.42}$ [1], which leads to $\Delta N_{\rm eff}=(N_{\rm
eff})_{\rm P18}-(N_{\rm eff})_{\rm SM}\leq 0.276$. The energy density of the
thermal axion at decoupling is given by $\rho_{a}^{\rm th}(T_{\rm
dec}^{a})=\rho_{\gamma}(T_{\rm dec}^{a})/2=\pi^{2}(T_{\rm dec}^{a})^{4}/30$
where $T_{\rm dec}^{a}$ can be estimated by evaluating
$3H\approx\tilde{\Gamma}_{a}$. Using Eq. (19) one obtains $T_{\rm
dec}^{a}\approx 300f_{a}^{2}/M_{\rm Pl}$. The thermal axion’s contribution to
the $\Delta N_{\rm eff}$ is proportional to

$\left.\frac{\rho_{a}^{\rm th}}{\rho_{\gamma}}\right|_{T_{\rm
rec}}=\frac{1}{2}\left(\frac{g_{\star S}(T_{\rm rec})}{g_{\star S}(T_{\rm
dec}^{a})}\right)^{4/3}\;,$ (27)

where $g_{\star S}(T)$ is the relativistic degrees of freedom in entropy. This
contribution has $f_{a}$ dependence through $T_{\rm dec}^{a}$. In the Gemini
DM model, we require $f_{a}>10^{9}\,{\rm GeV}$ (see Fig. 2), this translates
to a $\Delta N_{\rm eff}\approx 0.028$. As $f_{a}$ increases, the decoupling
temperature $T_{\rm dec}^{a}$ increases, which indicates a larger $g_{\star
S}(T_{\rm dec}^{a})$ and a smaller $\Delta N_{\rm eff}$, with a minimum value
of $\approx 0.027$ assuming no extra BSM degrees of freedom in thermal
equilibrium at early times. This constitutes a falsifiable prediction for the
Gemini DM model because CMB probes of the next generation are expected to
reach a sensitivity at or below this [20, 21]. Furthermore, there are two
possible additional ways in which the axion energy density can be enhanced to
levels that are already ruled out.

  1. (i)

In addition to the decay that produces $\chi_{3}$, the $\phi$ particle also
decays into axion via $\phi\to a+a$. The decay term is generated from the
kinetic term of $\Phi$ and is therefore likely to dominate. As long as this
decay occurs when $a$ is still in thermal equilibrium, the effect will be
washed out. The production from $\phi$-decay continuously happens until
$n_{\phi}$ is exponentially suppressed when $T\lesssim m_{\phi}/25$.
Therefore, as long as $T_{\rm dec}^{a}<m_{\phi}/25$, there is no appreciable
additional production of the dark radiation from $\phi$ decay. Therefore, we
avoid the overproduction of the dark radiation from $\phi$ decay by putting an
upper bound for $f_{a}$, indicated by the orange exclusion region in Fig 4.

  2. (ii)

The decay $\chi_{3}\to\chi_{1/2}+a$ can also contribute to dark radiation. It
occurs roughly at the temperature
$3H(T_{3})\approx\Gamma_{\chi_{3}\to\chi_{2}a}$, which is,

$T_{3}\approx 9.9\times 10^{-4}\frac{\sqrt{m_{3}^{3}M_{\rm Pl}}}{f_{a}}\;.$
(28)

For the parameter space that we are interested in, we usually have
$T_{3}<T_{\rm dec}^{a}$. So this extra contribution could not be avoided. To
estimate the additional contribution we assume instantaneous decay, i.e.
$\Delta\rho_{a}(T_{3})\approx m_{3}Y_{3}^{\rm f.i.}s(T_{3})/2$, where the half
indicates the energy sharing between $\chi_{1/2}$ and $a$. This gives an
additional contribution to the $\Delta N_{\rm eff}$ at the recombination,
which is proportional to

$\displaystyle\left.\frac{\Delta\rho_{a}}{\rho_{\gamma}}\right|_{T_{\rm rec}}$
$\displaystyle=\frac{\Delta\rho_{a}(T_{3})}{\rho_{\gamma}(T_{\rm
rec})}\left(\frac{T_{\rm rec}}{T_{3}}\right)^{4}\left(\frac{g_{\star S}(T_{\rm
rec})}{g_{\star S}(T_{3})}\right)^{4/3}$ $\displaystyle\approx 6.9\times
10^{6}\frac{T_{3}}{m_{\phi}}\left(\frac{g_{\star S}(T_{\rm rec})}{g_{\star
S}(T_{3})}\right)^{4/3}\;,$ (29)

where we have neglected the phase space factor from Eq. (8) in the second line
since we are mostly interested in the limit $m_{\phi}\gg m_{3}$. This energy
density ratio depends on $m_{3}$ and $f_{a}$ through $T_{3}$ given in Eq.
(28). The constraint demands an upper bound of this ratio, which is translated
to a lower bound on $f_{a}$ in terms of $m_{3}$.

Therefore, the total radiation $\rho_{a}=\rho_{a}^{\rm th}+\Delta\rho_{a}$ at
recombination gives an $\Delta N_{\rm eff}$ that proportional to the sum of
Eq. (27) and (29). We show a region in Fig. 4 that is excluded by current CMB
results on $\Delta N_{\rm eff}>0.276$ reported by Planck [1]. Once again, the
above calculation requires $T_{\rm dec}^{a}>T_{3}$, otherwise, the decay
axions are thermalised and $\Delta N_{\rm eff}\approx 0.027$. This is the
reason the low $f_{a}$ region is not covered by $\Delta N_{\rm eff}$
constraint in Fig. 4. We have additionally shaded a region in light purple
which shows where $\Delta N_{\rm eff}>0.04$, we do this to show the level of
sensitivity future CMB probes will require in order to rule out this solution
to the $S_{8}$ tension.

Dark matter models produced from thermal processes at the with $m_{\rm
DM}\sim\mathcal{O}(10)\,{\rm keV}$ are typically highly constrained by
inferring the matter power spectrum through measurements of the Lyman-$\alpha$
forest [50, 51]. To see whether the WDM constraint applies, one has to
estimate the free-streaming scale. In the instantaneous decay approximation of
the process $\chi_{3}\to\chi_{1/2}+a$, the twins $\chi_{1/2}$ obtain the
average momentum $\langle p\rangle_{3}\approx m_{3}/2$ at temperature $T_{3}$
(28). This momentum is redshifted to today and becomes $\langle p\rangle_{0}$,
rendering the twins non-relativistic. Their average velocity today can be
expressed as,

$\displaystyle\langle v\rangle_{0}$ $\displaystyle\approx\frac{\langle
p\rangle_{0}}{m_{1/2}}\approx\frac{m_{3}}{2m_{1/2}}\frac{T_{0}}{T_{3}}\left(\frac{g_{\star
S}(T_{0})}{g_{\star S}(T_{3})}\right)^{1/3}$ (30) $\displaystyle\approx
1.4\times 10^{-6}\left(\frac{m_{\phi}}{10^{9}\,{\rm GeV}}\right)^{-1/4}\;.$

For the last approximation, we have evaluated the number under the benchmark
model applying Eq. (10) and Eq. (24) with fixed $\Omega_{\rm DM}h^{2}=0.12$
and $\tau_{8}=10\,{\rm Gyr}$, which corresponds to the parameters located in
the gray strip in Fig. 4. The dependence on $m_{3}$ and $f_{a}$ canceled out
after assuming that $m_{3}\ll m_{\phi}$. The dependence on $m_{\phi}$ is weak.
This indicates that the Gemini DM is typically very cold today. With this
result, we can estimate the free-streaming scale of the Gemini DM, which is

$\lambda_{\rm fs}\approx\frac{\langle v\rangle_{0}}{H_{0}}\approx 4.1\times
10^{-3}\,{\rm Mpc}/h\;.$ (31)

This free-streaming scale is much smaller than $\mathcal{O}(1)\,{\rm Mpc}$,
the scale of Ly-$\alpha$ matter power spectrum [52]. Hence, the usual keV WDM
constraint [53] does not apply to Gemini DM [54, 55].

##  IV Summary and discussion

In this work, we construct a decaying dark matter model, named Gemini DM, that
resolves the $S_{8}/\sigma_{8}$ tension in the large-scale structure, by
extending a consistent Froggatt-Nielsen symmetry, that explains the fermion
mass hierarchy, to the dark sector. Our model consists of three generations of
dark fermions, among which two of them have almost degenerate mass (the twins)
and the other one (the mother) is much heavier. The twins are the main dark
matter component today and the mother particle gives birth to them through its
decay. The production mechanisms of the mother particle is nonthermal freeze-
in from the Standard Model thermal bath. We have the mass of the twins roughly
$m_{1/2}\sim\mathcal{O}({\rm keV})$, which is small. However, because they
never thermalize and are produced via $\chi_{3}$ decay, they constitute a cold
dark matter model. Dark radiation is also produced during the production of
the Gemini DM, which is the prediction of our model. This radiation is
expected be probed for future CMB observations.

Unlike the usual DDM solution to the $S_{8}$ tension, only a fraction of the
total dark matter decays. This will likely alter the required
$\epsilon\in(0.01,0.1)$, but not by much, for the following reason. The cold
dark matter decays into warm dark matter and suppresses structure formation
below the free-streaming scale $k_{\rm fs}^{-1}$. This suppression follows the
relation $\delta_{\rm WDM}\approx(k_{\rm fs}^{2}/k^{2})\delta_{\rm CDM}$ [56],
where the free-streaming scale is proportional to the particle velocity,
determined by the mass separation, i.e., $k_{\rm fs}^{-1}\propto\epsilon$. If
WDM is produced from only a fraction of the total dark matter, the suppression
of structure formation is $f_{\rm W}\delta_{\rm WDM}\propto f_{\rm W}k_{\rm
fs}^{2}\propto f_{\rm W}/\epsilon^{2}$, where $f_{\rm W}$ is the fraction of
WDM relative to the total DM. This implies that a smaller fraction of DDM
requires a larger energy release to achieve the same level of structure
suppression. Therefore, the required mass-splitting factor is
$\epsilon^{\prime}\equiv\epsilon/\sqrt{f_{\rm W}}$. For Gemini DM, where
$f_{\rm W}\approx 0.5$, this requires $\sqrt{2}$ times the energy compared to
the previously favored $\epsilon$ in [12]. Hence, adopting $\epsilon=0.05$ for
our benchmark model in the main text is reasonable.

One may attempt to identify the Gemini DM as the sterile neutrino. This could
be achieved by considering the additional coupling between lepton and $\chi$
particles. After the seesaw mechanism, one ends up with a suppressed neutrino
mass matrix and a heavy dark sector. We leave this for future study.

If the $S_{8}$ tension persists and the next generation of CMB measurements
find an indication of additional relativistic degrees of freedom $\Delta
N_{\rm eff}\approx 0.027$, we will have a strong candidate solution to both
simultaneously in the Gemini DM model.

###### Acknowledgements.

We thank Yong Du, Jie Sheng, Chuan-Yang Xing, and Jiang Zhu for their valuable
discussions. The authors are supported by the National Natural Science
Foundation of China (12375101, 12090060 and 12090064) and the SJTU Double
First Class start-up fund (WF220442604).

##  Appendix A Review on Froggatt-Nielsen mechanism

The FN mechanism explains the SM fermion hierarchy via the introduction of an
extra chiral $U(1)_{\rm FN}$ symmetry under some energy scale, which could be
either global or gauged. For simplicity, we adopt the global symmetry here. We
follow the Ref. [24], and one typical $U(1)_{\rm FN}$ charge assignment of the
SM fermions is presented in Table 2. Since FN charges are assigned for SM
particles, the usual Yukawa interactions are forbidden by $U(1)_{\rm FN}$. The
leading order terms that give fermion mass are non-renormalizable. These non-
renormalizable operators generate interactions between fermions and flavons.

Table 2: FN charge assignment of quarks and leptons. Generation $i$ | $\overline{Q}_{\rm L}$ | $u_{\rm R}$ | $d_{\rm R}$ | $\overline{\ell}_{\rm L}$ | $e_{\rm R}$ | $\Phi$
---|---|---|---|---|---|---
$1$ | $3$ | $4$ | $1$ | $1$ | $4$ | -1
$2$ | $2$ | $1.5$ | $0$ | $0.5$ | $1$
$3$ | $0$ | $0$ | $0$ | $0$ | $0$

We review the FN framework by describing how it works in the lepton sector,
the quark sector can be derived in the same way Below some UV cutoff scale
$\Lambda$, the leading order operators are

$-\mathcal{L}\supset
g_{ij}\left(\frac{\Phi}{\Lambda}\right)^{n_{\ell}^{ij}}\overline{\ell}_{\rm
L}^{i}\mathcal{H}e_{\rm R}^{j}+{\rm h.c.}\;,$ (32)

where $\mathcal{H}$ is the Higgs doublet and it does not carry an FN charge.
$g_{ij}$ is some order $1$ coupling that is not responsible for the hierarchy.
$n_{\ell}^{ij}=n_{\ell}^{i}+n_{e}^{j}$ is the lepton FN charge matrix, where
$n_{\ell}^{i}$ and $n_{e}^{j}$ are FN charges of $\overline{\ell}_{\rm L}^{i}$
and $e_{\rm R}^{j}$ respectively. Explicitly,

$n_{\ell}=\begin{pmatrix}5&2&1\\\ 4.5&1.5&0.5\\\ 4&1&0\\\ \end{pmatrix}\;.$
(33)

Note that here the half charge could be absorbed by redefining the FN
parameter. After the spontaneous symmetry breaking of $U(1)_{\rm FN}$ by
$\Phi$, we can do the same field expansion as Eq. (2). After performing the
same phase rotation as described in Sec. (II), the leading term in Eq.(32)
would be the usual Yukawa interactions The Yukawa matrix of the charged
leptons is now hierarchically textured,

$y_{\ell}\sim\begin{pmatrix}\mathcal{O}(\lambda^{5})&\mathcal{O}(\lambda^{2})&\mathcal{O}(\lambda^{1})\\\
\mathcal{O}(\lambda^{4.5})&\mathcal{O}(\lambda^{1.5})&\mathcal{O}(\lambda^{0.5})\\\
\mathcal{O}(\lambda^{4})&\mathcal{O}(\lambda^{1})&\mathcal{O}(1)\end{pmatrix}\;,$
(34)

where $\lambda\approx 0.171$ could give rise to a hierarchy between three
generations and provide the best fit with the observation [24]. Therefore,
after the electroweak symmetry breaking
$\mathcal{H}\to\langle\mathcal{H}\rangle$, one obtains lepton masses and
interactions with the flavon field,

$-\mathcal{L}\supset m_{k}\bar{e}_{\rm L}^{k}e_{\rm
R}^{k}+g_{e,ij}^{\phi}\phi\bar{e}_{\rm L}^{i}e_{\rm
R}^{j}+g_{e,ij}^{a}a\bar{e}_{\rm L}^{i}e_{\rm R}^{j}+{\rm h.c.}\;.$ (35)

Here we have rotated in the lepton mass eigenstates. In this basis, one
obtains the relation

$g_{e,ij}\propto\frac{m_{i}-m_{j}+m_{i}\delta_{ij}}{f_{a}}\;,$

which is characteristic of the FN framework. The diagonal couplings between
flavon and fermions are proportional to fermion mass, and the off-diagonal
couplings are determined by their mass difference. The quarks couples to
flavons basically in the same way as the charged leptons with up to slightly
different constants. Note that neutrino mass could be generated through dim-5
operator [57, 58] in the SM effective field theory. This operator should also
be modified under the FN framework [24], which gives rise to neutrino masses
and interactions with flavons in a similar manner.

## References

  * [1] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209 [astro-ph.CO]]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [2] E. Abdalla et al., “Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies,” JHEAp 34 (2022) 49–211, [arXiv:2203.06142 [astro-ph.CO]].
  * [3] E. Di Valentino et al., “Cosmology Intertwined III: $f\sigma_{8}$ and $S_{8}$,” Astropart. Phys. 131 (2021) 102604, [arXiv:2008.11285 [astro-ph.CO]].
  * [4] M. S. Turner, “Cosmology with Decaying Particles,” Phys. Rev. D 31 (1985) 1212.
  * [5] S. Aoyama, T. Sekiguchi, K. Ichiki, and N. Sugiyama, “Evolution of perturbations and cosmological constraints in decaying dark matter models with arbitrary decay mass products,” JCAP 07 (2014) 021, [arXiv:1402.2972 [astro-ph.CO]].
  * [6] K. Enqvist, S. Nadathur, T. Sekiguchi, and T. Takahashi, “Decaying dark matter and the tension in $\sigma_{8}$,” JCAP 09 (2015) 067, [arXiv:1505.05511 [astro-ph.CO]].
  * [7] G. Franco Abellán, R. Murgia, V. Poulin, and J. Lavalle, “Implications of the $S_{8}$ tension for decaying dark matter with warm decay products,” Phys. Rev. D 105 no. 6, (2022) 063525, [arXiv:2008.09615 [astro-ph.CO]].
  * [8] G. Franco Abellán, R. Murgia, and V. Poulin, “Linear cosmological constraints on two-body decaying dark matter scenarios and the S8 tension,” Phys. Rev. D 104 no. 12, (2021) 123533, [arXiv:2102.12498 [astro-ph.CO]].
  * [9] T. Simon, G. Franco Abellán, P. Du, V. Poulin, and Y. Tsai, “Constraining decaying dark matter with BOSS data and the effective field theory of large-scale structures,” Phys. Rev. D 106 no. 2, (2022) 023516, [arXiv:2203.07440 [astro-ph.CO]].
  * [10] L. Fuß and M. Garny, “Decaying Dark Matter and Lyman-$\alpha$ forest constraints,” JCAP 10 (2023) 020, [arXiv:2210.06117 [astro-ph.CO]].
  * [11] J. Bucko, S. K. Giri, F. H. Peters, and A. Schneider, “Probing the two-body decaying dark matter scenario with weak lensing and the cosmic microwave background,” Astron. Astrophys. 683 (2024) A152, [arXiv:2307.03222 [astro-ph.CO]].
  * [12] L. Fuß, M. Garny, and A. Ibarra, “Minimal decaying dark matter: from cosmological tensions to neutrino signatures,” [arXiv:2403.15543 [hep-ph]].
  * [13] C. D. Froggatt and H. B. Nielsen, “Hierarchy of Quark Masses, Cabibbo Angles and CP Violation,” Nucl. Phys. B 147 (1979) 277–298.
  * [14] M. Leurer, Y. Nir, and N. Seiberg, “Mass matrix models: The Sequel,” Nucl. Phys. B 420 (1994) 468–504, [arXiv:hep-ph/9310320].
  * [15] M. Leurer, Y. Nir, and N. Seiberg, “Mass matrix models,” Nucl. Phys. B 398 (1993) 319–342, [arXiv:hep-ph/9212278].
  * [16] L. Calibbi, A. Crivellin, and B. Zaldívar, “Flavor portal to dark matter,” Phys. Rev. D 92 no. 1, (2015) 016004, [arXiv:1501.07268 [hep-ph]].
  * [17] A. Cheek, J. K. Osiński, L. Roszkowski, and S. Trojanowski, “Dark matter production through a non-thermal flavon portal,” JHEP 03 (2023) 149, [arXiv:2211.02057 [hep-ph]].
  * [18] K. S. Babu, S. Chakdar, N. Das, D. K. Ghosh, and P. Ghosh, “FIMP dark matter from flavon portals,” JHEP 07 (2023) 143, [arXiv:2305.03167 [hep-ph]].
  * [19] R. Mandal and T. Tong, “Exploring freeze-out and freeze-in dark matter via effective Froggatt-Nielsen theory,” JCAP 11 (2023) 074, [arXiv:2307.14972 [hep-ph]].
  * [20] CMB-S4 Collaboration, K. N. Abazajian et al., “CMB-S4 Science Book, First Edition,” [arXiv:1610.02743 [astro-ph.CO]].
  * [21] CMB-HD Collaboration, S. Aiola et al., “Snowmass2021 CMB-HD White Paper,” [arXiv:2203.05728 [astro-ph.CO]].
  * [22] T. Banks and N. Seiberg, “Symmetries and Strings in Field Theory and Gravity,” Phys. Rev. D 83 (2011) 084019, [arXiv:1011.5120 [hep-th]].
  * [23] M. Reece, “TASI Lectures: (No) Global Symmetries to Axion Physics,” PoS TASI2022 (2024) 008, [arXiv:2304.08512 [hep-ph]].
  * [24] Y.-C. Qiu, J.-W. Wang, and T. T. Yanagida, “Predictions of mee and neutrino mass from a consistent Froggatt-Nielsen model,” Phys. Rev. D 108 no. 11, (2023) 115021, [arXiv:2307.16470 [hep-ph]].
  * [25] M. Fedele, A. Mastroddi, and M. Valli, “Minimal Froggatt-Nielsen textures,” JHEP 03 (2021) 135, [arXiv:2009.05587 [hep-ph]].
  * [26] L. Hui, “Wave Dark Matter,” Ann. Rev. Astron. Astrophys. 59 (2021) 247–289, [arXiv:2101.11735 [astro-ph.CO]].
  * [27] A. Lella, P. Carenza, G. Co’, G. Lucente, M. Giannotti, A. Mirizzi, and T. Rauscher, “Getting the most on supernova axions,” Phys. Rev. D 109 no. 2, (2024) 023001, [arXiv:2306.01048 [hep-ph]].
  * [28] C. B. Adams et al., “Axion Dark Matter,” in Snowmass 2021.  3, 2022.  [arXiv:2203.14923 [hep-ex]].
  * [29] F. Calore, A. Dekker, P. D. Serpico, and T. Siegert, “Constraints on light decaying dark matter candidates from 16 yr of INTEGRAL/SPI observations,” Mon. Not. Roy. Astron. Soc. 520 no. 3, (2023) 4167–4172, [arXiv:2209.06299 [hep-ph]].
  * [30] P. Gondolo and G. Gelmini, “Cosmic abundances of stable particles: Improved analysis,” Nucl. Phys. B 360 (1991) 145–179.
  * [31] M. Drees, F. Hajkarim, and E. R. Schmitz, “The Effects of QCD Equation of State on the Relic Density of WIMP Dark Matter,” JCAP 06 (2015) 025, [arXiv:1503.03513 [hep-ph]].
  * [32] B. Lillard, M. Ratz, T. Tait, M. P., and S. Trojanowski, “The Flavor of Cosmology,” JCAP 07 (2018) 056, [arXiv:1804.03662 [hep-ph]].
  * [33] J. McDonald, “Thermally generated gauge singlet scalars as selfinteracting dark matter,” Phys. Rev. Lett. 88 (2002) 091304, [arXiv:hep-ph/0106249].
  * [34] A. Kusenko, “Sterile neutrinos, dark matter, and the pulsar velocities in models with a Higgs singlet,” Phys. Rev. Lett. 97 (2006) 241301, [arXiv:hep-ph/0609081].
  * [35] K. Petraki and A. Kusenko, “Dark-matter sterile neutrinos in models with a gauge singlet in the Higgs sector,” Phys. Rev. D 77 (2008) 065014, [arXiv:0711.4646 [hep-ph]].
  * [36] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, “Freeze-In Production of FIMP Dark Matter,” JHEP 03 (2010) 080, [arXiv:0911.1120 [hep-ph]].
  * [37] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, “The Dawn of FIMP Dark Matter: A Review of Models and Constraints,” Int. J. Mod. Phys. A 32 no. 27, (2017) 1730023, [arXiv:1706.07442 [hep-ph]].
  * [38] M. Bauer and T. Plehn, Yet Another Introduction to Dark Matter: The Particle Physics Approach, vol. 959 of Lecture Notes in Physics.  Springer, 2019.  [arXiv:1705.01987 [hep-ph]].
  * [39] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [40] M. S. Turner, “Thermal Production of Not SO Invisible Axions in the Early Universe,” Phys. Rev. Lett. 59 (1987) 2489. [Erratum: Phys.Rev.Lett. 60, 1101 (1988)].
  * [41] A. Salvio, A. Strumia, and W. Xue, “Thermal axion production,” JCAP 01 (2014) 011, [arXiv:1310.6982 [hep-ph]].
  * [42] F. D’Eramo, F. Hajkarim, and S. Yun, “Thermal QCD Axions across Thresholds,” JHEP 10 (2021) 224, [arXiv:2108.05371 [hep-ph]].
  * [43] G. Mangano, G. Miele, S. Pastor, and M. Peloso, “A Precision calculation of the effective number of cosmological neutrinos,” Phys. Lett. B 534 (2002) 8–16, [arXiv:astro-ph/0111408].
  * [44] P. F. de Salas and S. Pastor, “Relic neutrino decoupling with flavour oscillations revisited,” JCAP 07 (2016) 051, [arXiv:1606.06986 [hep-ph]].
  * [45] S. Gariazzo, P. F. de Salas, and S. Pastor, “Thermalisation of sterile neutrinos in the early Universe in the 3+1 scheme with full mixing matrix,” JCAP 07 (2019) 014, [arXiv:1905.11290 [astro-ph.CO]].
  * [46] J. Froustey, C. Pitrou, and M. C. Volpe, “Neutrino decoupling including flavour oscillations and primordial nucleosynthesis,” JCAP 12 (2020) 015, [arXiv:2008.01074 [hep-ph]].
  * [47] J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor, and Y. Y. Y. Wong, “Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED,” JCAP 04 (2021) 073, [arXiv:2012.02726 [hep-ph]].
  * [48] M. Drewes, Y. Georis, M. Klasen, L. P. Wiggering, and Y. Y. Y. Wong, “Towards a precision calculation of N eff in the Standard Model. Part III. Improved estimate of NLO contributions to the collision integral,” JCAP 06 (2024) 032, [arXiv:2402.18481 [hep-ph]].
  * [49] V. A. Rubakov and D. S. Gorbunov, Introduction to the Theory of the Early Universe: Hot big bang theory.  World Scientific, Singapore, 2017.
  * [50] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese, and A. Riotto, “Constraining warm dark matter candidates including sterile neutrinos and light gravitinos with WMAP and the Lyman-alpha forest,” Phys. Rev. D 71 (2005) 063534, [arXiv:astro-ph/0501562].
  * [51] A. Boyarsky, J. Lesgourgues, O. Ruchayskiy, and M. Viel, “Lyman-alpha constraints on warm and on warm-plus-cold dark matter models,” JCAP 05 (2009) 012, [arXiv:0812.0010 [astro-ph]].
  * [52] R. Murgia, A. Merle, M. Viel, M. Totzauer, and A. Schneider, “”Non-cold” dark matter at small scales: a general approach,” JCAP 11 (2017) 046, [arXiv:1704.07838 [astro-ph.CO]].
  * [53] M. Drewes et al., “A White Paper on keV Sterile Neutrino Dark Matter,” JCAP 01 (2017) 025, [arXiv:1602.04816 [hep-ph]].
  * [54] M. Garny and J. Heisig, “Interplay of super-WIMP and freeze-in production of dark matter,” Phys. Rev. D 98 no. 9, (2018) 095031, [arXiv:1809.10135 [hep-ph]].
  * [55] Q. Decant, J. Heisig, D. C. Hooper, and L. Lopez-Honorez, “Lyman-$\alpha$ constraints on freeze-in and superWIMPs,” JCAP 03 (2022) 041, [arXiv:2111.09321 [astro-ph.CO]].
  * [56] A. Ringwald and Y. Y. Y. Wong, “Gravitational clustering of relic neutrinos and implications for their detection,” JCAP 12 (2004) 005, [arXiv:hep-ph/0408241].
  * [57] S. Weinberg, “Baryon and Lepton Nonconserving Processes,” Phys. Rev. Lett. 43 (1979) 1566–1570.
  * [58] T. Yanagida, “Horizontal gauge symmetry and masses of neutrinos,” Conf. Proc. C 7902131 (1979) 95–99.


# Bidding Agent Design in the LinkedIn Ad Marketplace

Yuan Gao, Kaiyu Yang, Yuanlong Chen, Min Liu
LinkedIn Corporation
Mountain View, CA 94043
{yugao,kyang,yuachen<EMAIL_ADDRESS>
&Noureddine El Karoui
Berkeley, CA 94720
<EMAIL_ADDRESS>
Work done in part while at LinkedIn, Inc. and UC, Berkeley.

###### Abstract

We establish a general optimization framework for the design of automated
bidding agent in dynamic online marketplaces. It optimizes solely for the
buyer’s interest and is agnostic to the auction mechanism imposed by the
seller. As a result, the framework allows, for instance, the joint
optimization of a group of ads across multiple platforms each running its own
auction format. Bidding strategy derived from this framework automatically
guarantees the optimality of budget allocation across ad units and platforms.
Common constraints such as budget delivery schedule, return on investments and
guaranteed results, directly translates to additional parameters in the
bidding formula. We share practical learnings of the deployed bidding system
in the LinkedIn ad marketplace based on this framework.

##  1 Introduction

There is a growing trend of automation in online advertising. The manual task
of creating tailored ad for each platform and splitting budget across these
platforms are largely automated by programmatic solutions in demand-side
platforms (DSPs). The LinkedIn Marketing Solutions is a DSP that allows
advertisers to reach LinkedIn members through the LinkedIn Feed as well as its
audience network consists of trusted third party publishers. Impression
opportunities on these various placements are normally sold through different
sorts of auctions. For instance, the LinkedIn Feed placement employs
generalized second price (GSP) auctions with reserve price, while majority of
publishers in the audience network now uses first price auctions. The paper
deals with the design of an automated agent that places bid on behalf of
advertisers in these different types of auctions. The goal is to help
advertisers maximize their results under natural constraints such as budget
delivery plan and expected return on investment (ROI).

Typical advertising objectives include impression, reach, click, social
engagement, video view, lead and conversion. Advertisers can directly encode
their valuations of the objectives into the bids for each impression
opportunity. For example, advertisers who value $5 per click for their ads
could place a bid of $5 * p(click) in a second price auction, where p(click)
is the probability that the user would click after seeing the ad. However this
bidding strategy is blind to the budget. This could lead to early termination
if the budget is constrained. Throttling-based systems are introduced Agarwal
et al. (2014); Xu et al. (2015) to address this problem by probabilistically
blocking the ad from entering into auctions as a way to smooth out the spend
over the entire lifetime. Yet a better strategy in this case is to lower the
bids so that budget is spent more effectively.

To formally model this problem, in economics literature Balseiro et al.
(2015); Balseiro and Gur (2019) a profit-maximization setup is often adopted,
where profit is defined as the difference between valuation and payment. In
practice, however, valuations from advertisers are not always available.
Defining valuations itself is a difficult problem given the scale and
complexity of the targeted audience. A preferred way is to frame the bidding
problem in a result maximization setup under budget constraint. Valuations
from advertisers, if available, can be encoded in this setup via additional
constraints on the bids or ROIs.

Mathematically such a bidding problem with budget constraint can be cast in a
stochastic optimal control framework Amin et al. (2012); Gummadi et al.
(2013); Balseiro et al. (2015); Fernandez-Tapia et al. (2017). When budget is
relatively large compared to the magnitude of bids, fluid approximation
Gallego and Van Ryzin (1994); Balseiro et al. (2015); Gummadi et al. (2013);
Fernandez-Tapia et al. (2017) is usually applied which greatly simplifies the
solution. Explicit bidding formulas are given in Zhang et al. (2014) under
further assumptions on stationary market competition. In practice though,
competition in online ad marketplace is highly dynamic due to the change in
demand and supply. For second price auctions, dynamic bid update strategies
without stationary assumptions are studied in Balseiro and Gur (2019);
Karlsson (2020).

Instead of directly optimizing for individual advertiser’s interest, there is
rich literature dealing with the optimization towards objectives of the
platform or DSP, by means of online resource allocations Mehta et al. (2007);
Aggarwal et al. (2019); Chen et al. (2011); Abrams et al. (2008); Balseiro et
al. (2020); Grigas et al. (2017); Feldman et al. (2010). Bidding strategy for
each advertiser in turn arises as a by-product of this optimization. These
centralized mechanisms and their equilibrium studies are not considered in our
design.

There is limited study on the optimal budget allocation problem across
multiple placements. In the case when all placements run second price
auctions, a stochastic bandit based solution is proposed Avadhanula et al.
(2021) where the value for each impression opportunity is assumed to be
unknown. This is an overly pessimistic assumption since in reality the value
is usually given by a response prediction model. A closely related problem is
studied in Zhang et al. (2016), in which the solution relies on feedback
controllers towards pre-calculated ROI target for each placement.

In Section 2 we lay out bidding agent design under a general framework of
result maximization under budget constraint. Online methods for solving this
optimization problem are discussed in Section 3. Section 4 demonstrates that
features such as cost control, budget delivery plan and guaranteed delivery
easily translates to constraints in this optimization framework, so that they
can be incorporated holistically in one bidding formula.

It’s worth noting that the framework is inherently compatible with a
collection of ad units and a mixture of different types of selling mechanisms,
thereby allowing for the joint optimization of a group of ads across multiple
placements each running (possibly) a different auction format. In particular,
in Section 5 we show that optimal bidding solution derived from this
formulation automatically guarantees the optimality of budget allocation
across ad units and placements.

Finally, in Section 6 explicit formulas for the starting bid is derived and
various issues around the implementation of online optimization methods are
discussed.

##  2 General Formulation

In this section we establish a general optimization framework for a single
advertiser. For simplicity, the problem is formulated as maximizing total
result under a single budget constraint. As noted later in Section 4 it is
straightforward to incorporate other types of constraints as well.

Suppose there are a total of $T$ impression opportunities, triggered from
targeted users’ search or view sessions. The $t$-th impression has a value of
$v_{t}$ for the agent, which is usually produced from a user response
prediction model tied to the advertiser’s marketing objective. As an example,
for an advertiser optimizing for conversions, $v_{t}$ represents the
impression to conversion rate. Denote $G_{t}(b)$ as the probability of winning
the $t$-th impression under a bid price of $b$ and $H_{t}(b)$ as the expected
cost. With an advertising budget $B$ the optimization problem becomes

$\max_{b_{t}}\sum_{t=1}^{T}v_{t}G_{t}(b_{t}),\quad s.t.\
\sum_{t=1}^{T}H_{t}(b_{t})\leq B.$ (1)

In a complete information setting, we can model
$G_{t}=\mathbf{1}_{\left[b_{t}\geq c_{t}\right]}$ and
$H_{t}=d_{t}\mathbf{1}_{\left[b_{t}\geq c_{t}\right]}$, where $c_{t}$ is the
minimum bid to win impression $t$, and $d_{t}$ is the cost if won. This
formulation leads to a Knapsack problem Chakrabarty et al. (2008); Zhou and
Naroditskiy (2008). In practice it’s often convenient to assume distributional
foresight, where $G_{t}$ and $H_{t}$ are modeled as smooth functions. Denote
their derivatives as $g_{t}$ and $h_{t}$. We make the following assumptions:

  * •

$G_{t}(0)=0$, $G_{t}(+\infty)=1$ and $g_{t}>0$ (monotonicity).

  * •

$H_{t}(0)=0$ and $h_{t}\geq 0$.

These are mild assumptions that generally hold for all practical auction types
including first & second price auctions, generalized first price (GFP)
auctions , generalized second price (GSP) auctions, Vickrey–Clarke–Groves
(VCG) auctions and variants of them with floor price.

The Lagrangian for (1) is

$L(b_{t},\lambda)=\sum_{t=1}^{T}\left[v_{t}G_{t}(b_{t})-\lambda
H_{t}(b_{t})\right]+\lambda B.$ (2)

At fixed $\lambda>0$, we would choose $b_{t}$ that maximizes the surplus

$(v_{t}/\lambda)G_{t}(b_{t})-H_{t}(b_{t}),$ (3)

where $v_{t}/\lambda$ is the $\lambda$-adjusted value. Optimality condition in
$b_{t}$ implies

$v_{t}g_{t}(b^{*}_{t})=\lambda h_{t}(b^{*}_{t}),$ (4)

or equivalently,

$b_{t}^{*}(\lambda)=\left(\frac{h_{t}}{g_{t}}\right)^{-1}\left(\frac{v_{t}}{\lambda}\right),$
(5)

where $(\cdot)^{-1}$ denotes function inverse.

######  Remark 1 (Explicit formula for first and second price auctions).

In ad marketplace impression opportunities are usually sold via auctions, and
$H_{t}$ and $G_{t}$ are generally related in various kinds of auctions. In
second price auctions, for example,
$H_{t}(b_{t})=\int_{0}^{b_{t}}bg_{t}(b)db=b_{t}G_{t}(b_{t})-\int_{0}^{b_{t}}G_{t}(b)db$.
This implies $h_{t}(b_{t})=b_{t}g_{t}(b_{t})$. According to (4) The optimal
bidding strategy in second price auctions is simply bidding the
$\lambda$-adjusted value

$b_{t}^{*}(\lambda)=v_{t}/\lambda.$

In first price auctions $H_{t}(b_{t})=b_{t}G_{t}(b_{t})$. The surplus
maximization (3) simplifies to $\max_{b_{t}}(v_{t}/\lambda-
b_{t})G_{t}(b_{t})$. Since $h_{t}(b_{t})=G_{t}(b_{t})+b_{t}g_{t}(b_{t})$, (5)
gives the following optimal bidding strategy

$b_{t}^{*}(\lambda)=\left(I+\frac{G_{t}}{g_{t}}\right)^{-1}\left(\frac{v_{t}}{\lambda}\right),$

where $I$ stands for identity mapping. Note that $G_{t}/g_{t}\geq 0$, so we
always have $b_{t}^{*}(\lambda)\leq v_{t}/\lambda$. This means that in first
price auctions the bid has to be shaded on top of the $\lambda$-adjusted
value.

We define two useful quantities that will be used throughout.

$V_{t}(\lambda):=v_{t}G_{t}(b_{t}^{*}(\lambda)),$

the expected value obtained in the $t$-th opportunity when bidding optimally
under $\lambda$, and

$S_{t}(\lambda):=H_{t}(b_{t}^{*}(\lambda)),$

the expected spend in the $t$-th opportunity when bidding optimally under
$\lambda$. Let $V(\lambda):=\sum_{t=1}^{T}V_{t}(\lambda)$ be the total
expected value and $S(\lambda):=\sum_{t=1}^{T}S_{t}(\lambda)$ the total
expected spend.

The following proposition shows that their derivatives are linearly related,
and, under mild conditions, both of them are monotonic.

######  Proposition 1.

$V_{t}^{\prime}(\lambda)\equiv\lambda
S_{t}^{\prime}(\lambda),\forall\lambda>0$. If in addition $\forall h_{t}>0$,
we have $(\log h_{t})^{\prime}>(\log g_{t})^{\prime}$, then
$V_{t}^{\prime}(\lambda)\equiv\lambda S_{t}^{\prime}(\lambda)\leq
0,\forall\lambda>0$.

###### Proof.

The relation $V_{t}^{\prime}(\lambda)\equiv\lambda
S_{t}^{\prime}(\lambda),\forall\lambda>0$ is directly given by optimality
condition (4). Now apply implicit differentiation on (4), we get

$\frac{db_{t}^{*}}{d\lambda}=\frac{h_{t}(b_{t}^{*})}{v_{t}g_{t}^{\prime}(b_{t}^{*})-\lambda
h_{t}^{\prime}(b_{t}^{*})}.$

When $h_{t}(b_{t}^{*})=0$ it’s clear that $\frac{db_{t}^{*}}{d\lambda}=0$.
Otherwise, the condition $(\log h_{t})^{\prime}>(\log g_{t})^{\prime}$ implies

$\frac{h_{t}^{\prime}(b_{t}^{*})}{h_{t}(b_{t}^{*})}>\frac{g_{t}^{\prime}(b_{t}^{*})}{g_{t}(b_{t}^{*})}=\frac{g_{t}^{\prime}(b_{t}^{*})}{(\lambda/v)h_{t}(b_{t}^{*})}\Rightarrow\frac{db_{t}^{*}}{d\lambda}<0.$

Finally by chain rule, we have
$V_{t}^{\prime}(\lambda)=v_{t}g_{t}(b_{t}^{*})\frac{db_{t}^{*}}{d\lambda}\leq
0$. ∎

######  Remark 2.

In second price auctions, the extra condition $\forall h_{t}>0,(\log
h_{t})^{\prime}>(\log g_{t})^{\prime}$ in Proposition 1 always holds due to
the fact that $h_{t}(b_{t})=b_{t}g_{t}(b_{t})$. In first price auctions, it’s
easy to show that the condition translates to log-concavity in $G_{t}$, i.e.
$G_{t}G_{t}^{\prime\prime}\leq g_{t}^{2}$. Note that $G_{t}$ can be seen as
the cumulative density function (CDF) of competitors’ bid distribution,
therefore the condition holds as long as the distribution is log-concave.

Plugging $b_{t}^{*}(\lambda)$ in the Lagrangian (2), we arrive at the dual
problem.

$\min_{\lambda\geq 0}L(\lambda):=\sum_{t=1}^{T}\left[V_{t}(\lambda)-\lambda
S_{t}(\lambda)+\lambda\frac{B}{T}\right].$ (6)

Setting $L^{\prime}(\lambda^{*})=0$, and applying Proposition 1 we get

$\sum_{t=1}^{T}S_{t}(\lambda^{*})=B.$ (7)

Equation (7) says that an optimal solution $\lambda^{*}>0$ would result in a
match of the expected spend $S(\lambda)$ with budget. In other words, the
solution $\lambda^{*}$ is either 0 when budget is unconstrained, or one that
spends the budget, not surprising given Karush–Kuhn–Tucker (KKT) conditions
Kuhn and Tucker (2014).

In addition, Proposition 1 implies the monotonicity of spend $S(\lambda)$. Due
to this monotonicity, bisection method can be used to find $\lambda^{*}$ given
$G_{t}$, $H_{t}$ and $v_{t}$ for $t=1,\ldots,T$. For example, $\lambda^{*}$
hence the optimal bidding sequence $b_{t}$ can be obtained in hindsight given
historical auction logs. However in practice it’s difficult to generate
accurate forecast for $G_{t}$ and $H_{t}$ in dynamic marketplaces. In the next
section we consider online optimization methods in the dual that does not
require access to explicit models of these quantities.

##  3 Online Methods in the Dual

In this section we present various online optimization methods in the dual
(6). Denote the dual loss at impression opportunity $t$ as

$L_{t}(\lambda):=V_{t}(\lambda)-\lambda S_{t}(\lambda)+\lambda\frac{B}{T}.$

The dual objective (6) is then $\sum_{t}L_{t}(\lambda_{t})$. Derivative of the
online loss

$L_{t}^{\prime}(\lambda)=B/T-S_{t}(\lambda),$

is the difference between the average budget in each opportunity and the
expected spend at $t$.

###  3.1 Follow the Leader

In follow the leader (FTL) algorithm at each iteration we choose the best
solution in hindsight, i.e.,

$\lambda_{t+1}=\text{argmin}_{\lambda\geq
0}\sum_{\tau=1}^{t}L_{\tau}(\lambda).$

Similar to optimality condition (7) we have

$\frac{1}{t}\sum_{\tau=1}^{t}S_{\tau}(\lambda_{t+1})=\frac{B}{T}.$

In practice, this can be achieved by replaying the past auctions to find
$\lambda_{t+1}$ such that average cost per auction is $B/T$. As the problem
size becomes larger over time, other variants such as those using a fixed size
look back window from $t-\Delta t$ to $t$ might be preferred.

In the stationary case (i.e. $G_{t}=G,H_{t}=H$), the FTL approach converges to
the optimal solution in one iteration (i.e.
$\lambda_{2}=...=\lambda_{T}=\lambda^{*}$), but in reality $\lambda_{t}$
produced by FTL can be unstable when $G_{t}$ and $H_{t}$ are nonstationary.

###  3.2 Linearization and Follow the Regularized Leader

It’s obvious that in both first and second price auctions $L_{t}(\lambda)$ is
convex since $L_{t}^{\prime\prime}(\lambda)\geq 0$ according to Proposition 1.
Therefore the linearized loss around $\lambda_{t}$

$\tilde{L_{t}}(\lambda):=L_{t}(\lambda_{t})+L^{\prime}_{t}(\lambda_{t})(\lambda-\lambda_{t})$

is an lower bound of $L_{t}(\lambda)$. In follow the regularized leader (FTRL)
we choose

$\lambda_{t+1}=\text{argmin}_{\lambda\geq
0}\sum_{\tau=1}^{t}\tilde{L}_{\tau}(\lambda)+\Omega(\lambda),$

where $\Omega(\lambda)$ is a regularization term.

######  Theorem 1 (Dual Online Mirror Descent).

When $\Omega(\lambda)=\lambda^{2}/(2\epsilon)$, the online update is additive
$\lambda_{t+1}=\lambda_{t}-\epsilon L^{\prime}_{t}(\lambda_{t})$; when
$\Omega(\lambda)=\lambda\ln\lambda/\epsilon$, the online update is
multiplicative
$\lambda_{t+1}=\lambda_{t}e^{-\epsilon\L^{\prime}_{t}(\lambda_{t})}$.

Second order approximations to $L(\lambda)$ can be employed to obtain faster
convergence. In particular,
$L^{\prime\prime}(\lambda)=-\sum_{t=1}^{T}S^{\prime}_{t}(\lambda)=-S^{\prime}(\lambda)$
represents the sensitivity of expected spend around $\lambda$. Intuitively
when the variation of expected spend is small one can make larger updates to
$\lambda$. Online optimization methods that use second order information, such
as natural gradient descent Amari et al. (2000) or online Newton’s method
Hazan et al. (2007), can be used in the update.

##  4 Other Types of Constraints

Apart from the overall budget constraint, in reality there are sometimes
preferences on budget delivery schedule as well. Advertisers may also prefer
some sort of cost control, e.g., keeping the cost per result under a certain
threshold. Other types of constraints include the request to ensure certain
number of results are delivered within a time range. In this section we show
it’s straightforward to incorporate them in the general formulation.

###  4.1 Cost Control

A cost per result target $C$ could be specified representing a given ROI goal.
The agent attempts to control its cost per result under $C$ by solving the
following problem

$\max_{b_{t}}\sum_{t=1}^{T}v_{t}G_{t}(b_{t}),\quad s.t.\
\sum_{t=1}^{T}H_{t}(b_{t})\leq\min\left(B,C\sum_{t=1}^{T}v_{t}G_{t}(b_{t})\right).$

Equipping a multiplier for each constraint, the Lagrangian becomes

$\sum_{t=1}^{T}\left[(1+\mu
C)v_{t}G_{t}(b_{t})-(\lambda+\mu)H_{t}(b_{t})\right]+\lambda B.$

Solution in the optimal bidding formula $b_{t}$ is similar, where

$b_{t}^{*}(\lambda,\mu)=\left(\frac{h_{t}}{g_{t}}\right)^{-1}\left(\frac{1+\mu
C}{\lambda+\mu}\cdot v_{t}\right).$

Again KKT conditions imply either $\lambda^{*}=0$ (in which case the budget is
not a binding constraint) or the budget is spent, i.e.,

$\sum_{t=1}^{T}H_{t}(b_{t}^{*}(\lambda^{*},\mu^{*}))=B.$

Similarly, either $\mu^{*}=0$ (in which case the cost per result target is not
a binding constraint), or the cost per result is equal to $C$, i.e.,

$\sum_{t=1}^{T}H_{t}(b_{t}^{*}(\lambda^{*},\mu^{*}))=C\sum_{t=1}^{T}v_{t}G_{t}(b_{t}^{*}(\lambda^{*},\mu^{*})).$

###  4.2 Budget Delivery Control

We could define budget delivery constraints on some subintervals
$T_{k},k=1,\ldots,K$, where $\sum_{t\in T_{k}}\leq T$. This can reflect
advertisers’ delivery preferences, for example, to limit the spend on
weekends.

$\max_{b_{t}}\sum_{t=1}^{T}v_{t}G_{t}(b_{t}),\
s.t.\sum_{t=1}^{T}H_{t}(b_{t})\leq B,\sum_{t\in T_{k}}H_{t}(b_{t})\leq
B_{k},\forall k.$

In practice these constraints can also be dynamic, where they can be added and
adjusted in realtime based on inputs from advertisers. The Lagrangian is then

$\sum_{k}\sum_{t\in
T_{k}}\left[v_{t}G_{t}(b_{t})-(\lambda+\lambda_{k})H_{t}(b_{t})\right]+\sum_{k}\lambda_{k}B_{k}+\lambda
B,$

where $\lambda_{k}\geq 0$. The new constraints influence the bidding formula.
For $t\in T_{k}$,

$b_{t}^{*}(\lambda,\lambda_{k})=\left(\frac{h_{t}}{g_{t}}\right)^{-1}\left(\frac{v_{t}}{\lambda+\lambda_{k}}\right).$

Note that each $\lambda_{k}$ is only active during the period of $T_{k}$,
i.e., when the constraint is in place.

###  4.3 Guaranteed Delivery

Advertisers may also request that certain number of results to be delivered,
for instance, during holiday season. Instead of manually increasing the budget
during the time period, more precise control can be achieved by encoding the
requirement as constraints directly in the optimization framework. More
formally, suppose we’d like $V_{k}$ number of results to be delivered in each
subinterval $T_{k}$ (defined in the previous section), then the problem
becomes

$\max_{b_{t}}\sum_{t=1}^{T}v_{t}G_{t}(b_{t}),\
s.t.\sum_{t=1}^{T}H_{t}(b_{t})\leq B,\sum_{t\in T_{k}}v_{t}G_{t}(b_{t})\geq
V_{k},\forall k.$

The Lagrangian is now

$\sum_{k}\sum_{t\in T_{k}}\left[(1+\mu_{k})v_{t}G_{t}(b_{t})-\lambda
H_{t}(b_{t})\right]-\sum_{k}\mu_{k}V_{k}+\lambda B.$

For $t\in T_{k}$, the bidding formula is

$b_{t}^{*}(\lambda,\mu_{k})=\left(\frac{h_{t}}{g_{t}}\right)^{-1}\left(\frac{1+\mu_{k}}{\lambda}\cdot
v_{t}\right).$

Since $\mu_{k}\geq 0$, a nonzero multiplier $\mu_{k}^{*}$ would give a boost
to the bid during the period $T_{k}$ to help achieve the delivery requirement.

A combination of the constraints can be implemented at the same time. As an
illustration, if all constraints discussed in this section are involved, the
bidding formula would then become

$b_{t}^{*}(\lambda,\mu,\lambda_{k},\mu_{k})=\left(\frac{h_{t}}{g_{t}}\right)^{-1}\left(\frac{1+\mu
C+\mu_{k}}{\lambda+\lambda_{k}+\mu}\cdot v_{t}\right).$

This offers a systematic way to handle a complex set of constraints
simultaneously.

##  5 Multiple Placements, Group of Ads and Equality of Marginal ROI

In advertising, each impression opportunity is associated with a placement.
For instance, an advertiser might want to serve his/her ad on multiple
publishers at the same time, where each publisher website is a placement. Our
framework naturally allows the simultaneous optimization on multiple
placements as equation (1) does not differentiate the $t$-th impression
opportunity based on its placement. Concretely, suppose there are a total of
$K$ placements and $T_{k},k=1,\ldots,K$ represents the set of impression
opportunities associated with the $k$-th placement, with $\sum_{t\in
T_{k}}=T$. Then (1) is essentially

$\max_{b_{t}}\sum_{k=1}^{K}\sum_{t\in T_{k}}v_{t}G_{t}(b_{t}),\quad s.t.\
\sum_{k=1}^{K}\sum_{t\in T_{k}}H_{t}(b_{t})\leq B.$

Note that the impression opportunities across different placements can be
interleaved in time. According to the Lagrangian (2) the total number of
results obtained under optimal bidding is

$L^{*}:=L(b^{*}_{t},\lambda^{*})=\sum_{k=1}^{K}\left[\lambda^{*}B_{k}+\sum_{t\in
T_{k}}V_{t}(\lambda^{*})-\lambda^{*}S_{t}(\lambda^{*})\right],$

where $B_{k}$ is the budget spent on the $k$-th placement. In this case
$\partial L^{*}/\partial B_{k}$, the marginal return on investment (ROI) for
all placements $1,\ldots,K$ are equal (to $\lambda^{*}$), which is a necessary
condition for optimal budget allocation across multiple placements.

Similar arguments apply to the scenario of optimizing a group of ads with a
global budget. Under the optimal bidding strategy in our framework, budgets
are allocated to each ad in an optimal manner to generate the most results. In
fact, one major advantage of this framework is its composability. Automated
budget allocation across placements, ad units as well as features described in
Section 4 are all elegantly handled via a unified bidding strategy.

##  6 Practical Considerations

In this section we aim to bridge the gaps between theory and practical
implementation. Important topics such as starting bid, forecasting error and
various issues around the online updates are discussed. In particular, Section
6.1 gives explicit formulas for the optimal starting bid based on statistics
from the targeted audience. Section 6.2 describes several variants of the
batch online gradient descent formula and applicability of them in different
scenarios.

###  6.1 Initialization of $\lambda$

Section 3 provides incremental update rules for the multiplier $\lambda$. The
cold start problem, namely the starting value of the multiplier, is also very
important as it impacts the time it takes to converge. On the macroscopic
level it also influences the price stability of the entire marketplace. Given
a set of targeted users, we can learn from these users’ past auction logs to
initialize $\lambda$. In the following theorem we provide explicit formulas
for second price auctions.

######  Theorem 2.

Suppose the competitors’ bids follow a log-normal distribution with parameters
$\mu$ and $\sigma$. Additionally, suppose the value of the ad (independently)
follows a log-normal distribution with parameters $\mu^{\prime}$ and
$\sigma^{\prime}$. Given a budget of $B$ and a total opportunity forecast of
$T$, the $\lambda^{*}$ for second price auctions is the solution to

$e^{\mu+\frac{\sigma^{2}}{2}}\Phi\left(\frac{\mu^{\prime}-\mu-\ln\lambda^{*}-\sigma^{2}}{\sqrt{(\sigma^{\prime})^{2}+\sigma^{2}}}\right)=\frac{B}{T},$

where $\Phi$ is the CDF of a standard normal distribution $\mathcal{N}(0,1)$.

In practice $\mu$ and $\sigma$ are derived using auction logs from targeted
users. The $\mu^{\prime}$ and $\sigma^{\prime}$ come from a combination of
targeted users’ value distribution and ad-specific features. In the multiple
placement scenario, the estimation of these quantities can be done at each
placement level. The solution $\lambda^{*}$ will then be derived based on the
traffic forecast $T_{k}$ from each placement $k$ and a global budget using
Theorem 2.

###  6.2 Implementation of Online Methods

Section 3 provides various online update methods that converge to
$\lambda^{*}$. In practice we employ a batch version of the dual online mirror
descent algorithm, where the multiplier $\lambda_{t}$ is updated every time
interval $dt$. In the following we focus on the additive formula in Theorem 1
though the arguments apply to the multiplicative update as well. Let $N_{dt}$
be the number of impression opportunities in time $dt$, then the (batch)
online update becomes

$\lambda_{t+N_{dt}}=\lambda_{t}-\epsilon\left(\frac{B}{T}\cdot
N_{dt}-S_{dt}\right),$ (8)

where $S_{dt}$ is the sum of expected spend during the past update period
under $\lambda_{t}$. This can simply be set as the observed spend in $dt$.
However if the charge events are sparse, some estimator based on the observed
spend is needed to reduce the variance. (8) is equivalent to

$\lambda_{t+N_{dt}}=\lambda_{t}-\epsilon_{dt}\left(1-R_{dt}\right),$

where

$\epsilon_{dt}:=\epsilon\frac{BN_{dt}}{T}$

is the step size, and

$R_{dt}:=\frac{S_{dt}/N_{dt}}{B/T}$ (9)

is the ratio of average cost per opportunity and average budget per
opportunity.

The choice of $dt$ is important as it controls the tradeoff between frequency
of update and variance in $S_{dt}$. An alternative is to keep track of the
observed impression opportunities and to trigger an update after a fixed
amount of observations.

####  6.2.1 Normalization

The scale of $\lambda$ various a lot across advertisers, due to the diversity
in budget and targeting setups. We therefore choose a normalizing constant
$\lambda^{\prime}$ for each advertiser and define
$\tilde{\lambda}_{t}=\lambda_{t}/\lambda^{\prime}$ to make the update
dimensionless:

$\tilde{\lambda}_{t+N_{dt}}=\tilde{\lambda}_{t}-\eta_{dt}\left(1-R_{dt}\right),$

where

$\eta_{dt}:=\epsilon\frac{BN_{dt}}{\lambda^{\prime}T}$ (10)

becomes a dimensionless step size that is easier to tune in practice. However,
the choice of $\lambda^{\prime}$ for each advertiser plays an important role
in convergence if a global parameter $\xi:=\epsilon B/\lambda^{\prime}$ is
selected such that $\eta_{dt}=\xi N_{dt}/T$. To see that, note that the
inverse of $\lambda_{t}$ will be used to compute the bids in auctions, and

$\frac{1}{\lambda_{t+N_{dt}}}=\frac{1}{\lambda^{\prime}(\tilde{\lambda}_{t}-\eta_{dt}\left(1-R_{dt}\right))}=\frac{1/\lambda_{t}}{1-(\xi
N_{dt}/T)\left(1-R_{dt}\right)(\lambda^{\prime}/\lambda_{t})}.$

Therefore a $\lambda^{\prime}<\lambda^{*}$ would slow down the convergence
whereas a $\lambda^{\prime}>\lambda^{*}$ tends to create oscillatory behavior
asymptotically due to magnification of the noise in $1-R_{dt}$. In practice we
find setting the normalization factor as the initialization provided in
Section 6.1 works well.

####  6.2.2 Total Forecast vs. Relative Forecast

Sometimes it is easier to forecast the relative traffic pattern over time
rather than the absolute number of opportunities $T$. In that case one can use
$\hat{r}_{dt}$, the forecasted proportion of traffic in time $dt$, to replace
$N_{dt}$ and $T$ in the formula of $R_{dt}$ in (9) and $\eta_{dt}$ in (10).
The formula then becomes $R_{dt}=S_{dt}/(B\hat{r}_{dt})$ and step size
$\eta_{dt}=\epsilon B\hat{r}_{dt}/\lambda^{\prime}$. The estimator
$\hat{r}_{dt}$ can also be dynamically adjusted based on past observed traffic
$N_{dt}$.

####  6.2.3 Model Predictive Control

One can reset the average budget per opportunity target $B/T$ at each update.
A model predictive control (MPC) version of the update would set $R_{dt}$ as

$\frac{S_{dt}/N_{dt}}{(B-B_{t+N_{dt}})/(T-t-N_{dt})}.$

This version of the update would encourage budget exhaustion near the end of
the ad’s lifetime.

##  7 Experimental Results

Bidding agent based on this framework was implemented on the LinkedIn ad
marketplace and is compared with a prior feedback-control agent which adjusts
the bid to track the forecast traffic curve. For unbiased evaluation the two
agents are compared using the budget-split experimentation platform Liu et al.
(2020). We observe statistically significant increase of 8.25% in advertiser
ROI (with neutral platform revenue) using the new method.

##  8 Acknowledgement

We would like to thank Yi Zhang and Onkar Dalal for their support. We
appreciate Wen Pu, Qian Yao and Ricardo Salmon for helpful discussions.

## References

  * (1)
  * Abrams et al. (2008) Zoö Abrams, S Sathiya Keerthi, Ofer Mendelevitch, and John A Tomlin. 2008.  Ad delivery with budgeted advertisers: A comprehensive LP approach.  _Journal of Electronic Commerce Research_ 9, 1 (2008).
  * Agarwal et al. (2014) Deepak Agarwal, Souvik Ghosh, Kai Wei, and Siyu You. 2014\.  Budget pacing for targeted online advertisements at linkedin. In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_. 1613–1619.
  * Aggarwal et al. (2019) Gagan Aggarwal, Ashwinkumar Badanidiyuru, and Aranyak Mehta. 2019.  Autobidding with constraints. In _International Conference on Web and Internet Economics_. Springer, 17–30.
  * Amari et al. (2000) Shun-ichi Amari, Hyeyoung Park, and Kenji Fukumizu. 2000\.  Adaptive method of realizing natural gradient learning for multilayer perceptrons.  _Neural computation_ 12, 6 (2000), 1399–1409.
  * Amin et al. (2012) Kareem Amin, Michael Kearns, Peter Key, and Anton Schwaighofer. 2012.  Budget optimization for sponsored search: Censored learning in MDPs.  _arXiv preprint arXiv:1210.4847_ (2012).
  * Avadhanula et al. (2021) Vashist Avadhanula, Riccardo Colini Baldeschi, Stefano Leonardi, Karthik Abinav Sankararaman, and Okke Schrijvers. 2021.  Stochastic bandits for multi-platform budget optimization in online advertising. In _Proceedings of the Web Conference 2021_. 2805–2817.
  * Balseiro et al. (2020) Santiago Balseiro, Haihao Lu, and Vahab Mirrokni. 2020\.  The best of many worlds: Dual mirror descent for online allocation problems.  _arXiv preprint arXiv:2011.10124_ (2020).
  * Balseiro et al. (2015) Santiago R Balseiro, Omar Besbes, and Gabriel Y Weintraub. 2015\.  Repeated auctions with budgets in ad exchanges: Approximations and design.  _Management Science_ 61, 4 (2015), 864–884.
  * Balseiro and Gur (2019) Santiago R Balseiro and Yonatan Gur. 2019.  Learning in repeated auctions with budgets: Regret minimization and equilibrium.  _Management Science_ 65, 9 (2019), 3952–3968.
  * Chakrabarty et al. (2008) Deeparnab Chakrabarty, Yunhong Zhou, and Rajan Lukose. 2008\.  Online knapsack problems. In _Workshop on internet and network economics (WINE)_.
  * Chen et al. (2011) Ye Chen, Pavel Berkhin, Bo Anderson, and Nikhil R Devanur. 2011\.  Real-time bidding algorithms for performance-based display ad allocation. In _Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining_. 1307–1315.
  * Feldman et al. (2010) Jon Feldman, Monika Henzinger, Nitish Korula, Vahab S Mirrokni, and Cliff Stein. 2010\.  Online stochastic packing applied to display ad allocation. In _European Symposium on Algorithms_. Springer, 182–194.
  * Fernandez-Tapia et al. (2017) Joaquin Fernandez-Tapia, Olivier Guéant, and Jean-Michel Lasry. 2017.  Optimal real-time bidding strategies.  _Applied Mathematics Research eXpress_ 2017, 1 (2017), 142–183.
  * Gallego and Van Ryzin (1994) Guillermo Gallego and Garrett Van Ryzin. 1994.  Optimal dynamic pricing of inventories with stochastic demand over finite horizons.  _Management science_ 40, 8 (1994), 999–1020.
  * Grigas et al. (2017) Paul Grigas, Alfonso Lobos, Zheng Wen, and Kuang-chih Lee. 2017\.  Profit maximization for online advertising demand-side platforms.  In _Proceedings of the ADKDD’17_. 1–7.
  * Gummadi et al. (2013) Ramki Gummadi, Peter Key, and Alexandre Proutiere. 2013\.  Optimal bidding strategies and equilibria in dynamic auctions with budget constraints.  _Available at SSRN 2066175_ (2013).
  * Hazan (2019) Elad Hazan. 2019\.  Introduction to online convex optimization.  _arXiv preprint arXiv:1909.05207_ (2019).
  * Hazan et al. (2007) Elad Hazan, Amit Agarwal, and Satyen Kale. 2007.  Logarithmic regret algorithms for online convex optimization.  _Machine Learning_ 69, 2-3 (2007), 169–192.
  * Karlsson (2020) Niklas Karlsson. 2020\.  Feedback Control in Programmatic Advertising: The Frontier of Optimization in Real-Time Bidding.  _IEEE Control Systems Magazine_ 40, 5 (2020), 40–77.
  * Kuhn and Tucker (2014) Harold W Kuhn and Albert W Tucker. 2014.  Nonlinear programming.  In _Traces and emergence of nonlinear programming_. Springer, 247–258.
  * Liu et al. (2020) Min Liu, Jialiang Mao, and Kang Kang. 2020.  Trustworthy online marketplace experimentation with budget-split design.  _arXiv preprint arXiv:2012.08724_ (2020).
  * Mehta et al. (2007) Aranyak Mehta, Amin Saberi, Umesh Vazirani, and Vijay Vazirani. 2007.  Adwords and generalized online matching.  _Journal of the ACM (JACM)_ 54, 5 (2007), 22–es.
  * Xu et al. (2015) Jian Xu, Kuang-chih Lee, Wentong Li, Hang Qi, and Quan Lu. 2015.  Smart pacing for effective online ad campaign optimization. In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. 2217–2226.
  * Zhang et al. (2016) Weinan Zhang, Yifei Rong, Jun Wang, Tianchi Zhu, and Xiaofan Wang. 2016.  Feedback control of real-time display advertising. In _Proceedings of the Ninth ACM International Conference on Web Search and Data Mining_. 407–416.
  * Zhang et al. (2014) Weinan Zhang, Shuai Yuan, and Jun Wang. 2014.  Optimal real-time bidding for display advertising. In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_. 1077–1086.
  * Zhou and Naroditskiy (2008) Yunhong Zhou and Victor Naroditskiy. 2008.  Algorithm for stochastic multiple-choice knapsack problem and application to keywords bidding. In _Proceedings of the 17th international conference on world wide web_. 1175–1176.

##  Appendix A Proof of Theorem 1

This is standard result in online optimization Hazan (2019). Ignoring the
constant terms, we have

$\lambda_{t+1}=\text{argmin}_{\lambda\geq
0}\sum_{s=1}^{t}L^{\prime}_{s}(\lambda_{s})\cdot\lambda+R(\lambda).$

First order optimality condition gives
$R^{\prime}(\lambda_{t+1})=-\sum_{s=1}^{t}L^{\prime}_{s}(\lambda_{s})$. In
case $R(\lambda)=\lambda^{2}/(2\epsilon)\Rightarrow
R^{\prime}(\lambda)=\lambda/\epsilon$, then

$\lambda_{t+1}=-\epsilon\left(\sum_{s=1}^{t-1}L^{\prime}_{s}(\lambda_{s})+L^{\prime}_{t}(\lambda_{t})\right)=-\epsilon\left(\frac{\lambda_{t}}{-\epsilon}+L_{t}^{\prime}(\lambda_{t})\right)=\lambda_{t}-\epsilon
L^{\prime}_{t}(\lambda_{t}).$

Similarly, $R(\lambda)=\lambda\ln\lambda/\epsilon\Rightarrow
R^{\prime}(\lambda)=(\ln\lambda+1)/\epsilon$, then

$\lambda_{t+1}=e^{-1-\epsilon\sum_{s=1}^{t-1}L^{\prime}_{s}(\lambda_{s})-\epsilon
L^{\prime}_{t}(\lambda_{t})}=\lambda_{t}e^{-\epsilon
L^{\prime}_{t}(\lambda_{t})}.$

##  Appendix B Proof of Theorem 2

The following lemma is useful in the proof. Let $\phi$ and $\Phi$ be the PDF
and CDF of a standard normal distribution $\mathcal{N}(0,1)$.

######  Lemma 1.

Suppose $X\sim\mathcal{N}(0,1)$, then
$E[\Phi(aX+b)]=\Phi(\frac{b}{\sqrt{1+a^{2}}})$.

###### Proof.

$\displaystyle E[\Phi(aX+b)]$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{ax+b}\phi(y)dy\phi(x)dx$
$\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\mathbf{1}_{\left[y\leq
ax+b\right]}\cdot\phi(y)\phi(x)dydx.$

This is simply the probability of $Y\leq aX+b$ given that $X,Y$ are i.i.d
standard normally distributed variables. Since
$Y-aX\sim\mathcal{N}(0,1+a^{2})$, this completes the proof. ∎

Now we proceed with the proof of the theorem.

###### Proof.

Denote $p_{\mu,\sigma}(x)$ as the probability density function for
$Lognormal(\mu,\sigma)$. The expected cost per opportunity under a given
$\lambda$ is

$S(\lambda):=\int_{0}^{\infty}\left(\int_{0}^{v/\lambda}zp_{\mu,\sigma}(z)dz\right)p_{\mu^{\prime},\sigma^{\prime}}(v)dv,$

where at each given value $v$ we would bid $v/\lambda$ (see Remark 1). First
note that the inner integral $\int_{0}^{v/\lambda}zp_{\mu,\sigma}(z)dz$ is the
partial expectation of a log-normal random variable. By standard results

$\int_{0}^{v/\lambda}zp_{\mu,\sigma}(z)dz=e^{\mu+\frac{\sigma^{2}}{2}}\Phi\left(\frac{\ln
v-\ln\lambda-\mu-\sigma^{2}}{\sigma}\right).$

That implies
$S(\lambda)=e^{\mu+\frac{\sigma^{2}}{2}}E\left[\Phi\left(\frac{\ln
V-\ln\lambda-\mu-\sigma^{2}}{\sigma}\right)\right]$, where $V$ is r.v. with
$Lognormal(\mu^{\prime},\sigma^{\prime})$. Now with a change of variable
$X:=\frac{\ln V-\mu^{\prime}}{\sigma^{\prime}}\sim\mathcal{N}(0,1)$, we have

$S(\lambda)=e^{\mu+\frac{\sigma^{2}}{2}}E\left[\Phi\left(\frac{\sigma^{\prime}X+\mu^{\prime}-\mu-\ln\lambda-\sigma^{2}}{\sigma}\right)\right].$

Now apply the results from Lemma 1 ,

$S(\lambda)=e^{\mu+\frac{\sigma^{2}}{2}}\Phi\left(\frac{\mu^{\prime}-\mu-\ln\lambda-\sigma^{2}}{\sqrt{(\sigma^{\prime})^{2}+\sigma^{2}}}\right).$

∎