# Dynamics of Dark Matter Misalignment Through the Higgs Portal

Brian Batell<EMAIL_ADDRESS>Pittsburgh Particle Physics, Astrophysics, and
Cosmology Center, Department of Physics and Astronomy, University of
Pittsburgh, Pittsburgh, USA  Akshay Ghalsasi<EMAIL_ADDRESS>Pittsburgh
Particle Physics, Astrophysics, and Cosmology Center, Department of Physics
and Astronomy, University of Pittsburgh, Pittsburgh, USA  Mudit Rai
<EMAIL_ADDRESS>Pittsburgh Particle Physics, Astrophysics, and Cosmology Center,
Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, USA

###### Abstract

A light singlet scalar field feebly coupled through the super-renormalizable
Higgs portal provides a minimal and well-motivated realization of ultra-light
bosonic dark matter. We study the cosmological production of dark matter in
this model by elucidating the dynamics of two sources of scalar field
misalignment generated during the radiation era. For large scalar masses
(above ${\cal O}(10^{-3}\,{\rm eV}$)), dark matter is produced through thermal
misalignment, by which the scalar field is driven towards large field values
as a result of the finite-temperature effective potential. The dominance of
thermal misalignment in this mass range leads to a sharp relic abundance
prediction which, is to a significant extent, insensitive to the initial
conditions of the scalar field. On the other hand, for low mass scalars (below
${\cal O}(10^{-5}\,{\rm eV}$)), dark matter is produced via VEV misalignment,
which is caused by the induced scalar field vacuum expectation value triggered
by the electroweak phase transition. We show that the relic abundance in this
low mass range is sensitive to the scalar field initial conditions. In the
intermediate mass range, the relic abundance is a consequence of a competition
between thermal misalignment and VEV misalignment, leading to novel forced
resonance effects which cause a recurring enhancement and suppression in the
late time oscillation amplitude as a function of the scalar mass. We compare
our relic abundance predictions with constraints and projections from
equivalence principle and inverse square law tests, stellar cooling, resonant
molecular absorption, and observations of extra-galactic background light and
diffuse X-ray backgrounds. New experimental ideas are needed to probe most of
the cosmologically motivated regions of parameter space.

††preprint: PITT-PACC-2213

##  I Introduction

Understanding the nature of the cosmic dark matter (DM), which constitutes
roughly a quarter of the energy density in the universe Aghanim  _et al._
(2020), is among the most pressing problems in particle physics and cosmology
today. Despite its clear influence on a host of astrophysical and cosmological
phenomena, the most basic properties of DM remain mysterious, including it
fundamental dynamics (spin, mass, interactions, etc.) as well as its origin
during the earliest epochs of the universe.

One particularly compelling idea is that DM is an ultra-light, feebly-coupled
scalar field $\phi$ Antypas  _et al._ (2022). As is well-known, such scalar
field DM is generically created in the early universe through the misalignment
mechanism Preskill  _et al._ (1983); Abbott and Sikivie (1983); Dine and
Fischler (1983). The scalar field starts from some initial field value at the
end of inflation that is misaligned with respect to its eventual potential
minimum (i.e., its vacuum expectation value (VEV)). In the early stages of the
radiation-dominated era, the scalar field is held up by Hubble friction and,
notably, does not experience any additional dynamical misalignment during this
epoch. As the universe expands, the scalar field eventually begins oscillating
once the Hubble expansion rate drops below its mass. The oscillating scalar
field subsequently forms a pressureless, non-relativistic fluid and thus has
the desired bulk properties of cold DM. Thus, in this standard misalignment
scenario, the DM density today is controlled by the initial value $\phi_{i}$
and does not depend at all on the coupling of $\phi$ to the Standard Model
(SM).

A simple model realization of ultra-light bosonic DM consists of a real
singlet scalar field coupled to the SM through the super-renormalizable Higgs
portal, $A\phi H^{\dagger}H$, as first proposed more than a decade ago by
Piazza and Pospelov Piazza and Pospelov (2010). This model is attractive from
several perspectives. First, it is among the most minimal, UV-complete
extensions of the SM, with the cosmology (up to initial conditions) and
phenomenology determined by two parameters, namely, the mass of the scalar,
$m_{\phi}$ and its Higgs portal coupling $A$. Second, the feebleness of this
coupling allows the light scalar mass (and scalar potential) to be stable
against radiative corrections. Furthermore, the model has a distinctive
phenomenology, with a variety of existing probes from both terrestrial
experiments and astrophysical observations, and has served as a well-motivated
benchmark model for several proposed ultralight scalar DM detection concepts;
see, e.g., Refs. Graham  _et al._ (2016); Arvanitaki  _et al._ (2018); Antypas
_et al._ (2022).

On the other hand, the cosmology of this simple model remains relatively
underexplored, and it is the primary aim of this work to fill that gap. As our
study will make clear, the cosmology of the scalar field in the Higgs portal
model is not generally encapsulated by the standard misalignment scenario
discussed above. In contrast to the standard misalignment scenario, the scalar
field undergoes a nontrivial evolution during the radiation era by which new
sources of misalignment are dynamically generated, impacting the scalar relic
abundance in an essential way. In particular, we will investigate two distinct
dynamical misalignment mechanisms that are inherent in the Higgs portal model:
(i) thermal misalignment, and (ii) VEV misalignment.

The thermal misalignment mechanism was recently explored in detail in Ref.
Batell and Ghalsasi (2021) (see also Refs. Chun (2022); Brzeminski  _et al._
(2021) for related work) in the context of simple Yukawa-type theory, but the
essential aspects of the mechanism carry over to the Higgs portal model. In
the earliest stages of the radiation era, the scalar field responds to a
finite-temperature correction to its effective potential, associated with the
free energy density of the Higgs degrees of freedom in the thermal bath to
which the scalar is feebly coupled. This tends to drive the scalar field
toward its high-temperature minimum at large field values, thus generating
misalignment. As we will demonstrate below, thermal misalignment dominates for
scalar masses larger than a few meV. Furthermore, provided the initial field
value $\phi_{i}$ at the end of inflation is small (in magnitude) in comparison
to the displacement generated by thermal misalignment, the late-time scalar
oscillation amplitude and associated relic density are insensitive to initial
conditions, instead being tightly controlled by the DM mass and Higgs portal
coupling. Thus, a sharp prediction can be made for the model parameters
(scalar mass and Higgs portal coupling) yielding the observed relic abundance
through thermal misalignment, providing a cosmologically motivated target that
can be compared with experimental tests.

VEV misalignment is a second source of dynamical misalignment which is built
in to the Higgs portal model. It arises as a consequence of the electroweak
phase transition (EWPT), during which the Higgs field quickly turns on as
electroweak symmetry is broken and induces a rapid shift in the $\phi$
potential minimum towards its zero-temperature VEV. As we will see, VEV
misalignment dominates for scalar masses below about $10^{-5}$ eV. For this
region of parameter space, we will find that the scalar oscillation amplitude
and corresponding relic abundance sensitively depends on the initial field
value $\phi_{i}$, in particular, whether the scalar field is initially close
to its zero temperature VEV or far away from it.

In the intermediate scalar mass range, an intricate interplay between thermal
misalignment and VEV misalignment leads to a novel forced resonance effect,
with the scalar oscillation amplitude experiencing either an enhancement or
suppression following the EWPT depending on the $\phi$ mass. This manifests as
a striking series of peaks and valleys in the relic density target line as the
scalar mass is varied.

Our main results are presented in Fig. 1 as contours in the mass-coupling
parameter space yielding the observed DM relic abundance, providing
cosmologically motivated targets for experiments searching for ultralight
scalar DM. In view of this, we also compile the existing bounds and
sensitivity projections for a variety of experimental and astrophysical probes
of the model, including equivalence principle and inverse square law tests,
stellar cooling, resonant molecular absorption, and observations of extra-
galactic background light and diffuse X-ray backgrounds. Still, much of the
parameter space remains unexplored, and we hope our results will stimulate new
innovative experimental concepts for ultra-light DM searches.

Figure 1:  DM relic abundance predictions along with constraints and
projections displayed in the $m_{\phi}-A$ (mass-coupling) plane. Parameter
choices leading to the observed relic abundance, $\Omega_{\phi,0}=\Omega_{\rm
DM}=0.26$, computed by numerically evolving Eq. (29) and using Eq. (30), are
shown for two choices of initial conditions: $\phi_{i}=\phi_{0}$ (blue solid)
and $\phi_{i}=0$ (red solid). In Region I, thermal misalignment dominates and
the approximate relic abundance is given by Eq. (III.5.1) (purple dashed),
independently of the choice of initial conditions. In Region III, VEV
misalignment dominates over thermal misalignment, and the relic abundance is
sensitive to the initial conditions, being given by Eq. (III.5.3) (blue
dashed) for $\phi_{i}=\phi_{0}$ and by Eq. (III.5.3) (red dashed) for
$\phi_{i}=0$. In the intermediate mass range the relic density lines feature a
series of peaks and valleys, the result of a forced resonance phenomena due to
the EWPT (see inset plot). Also shown are bounds and projections from
equivalence principle tests (EP, green), inverse square law tests (ISL,
yellow), stellar cooling (orange), resonant absorption in molecules [SII (dark
red dashed) and BII (dark blue dashed)], extra-galactic background light (EBL,
light blue), and X-ray observatories (X-Ray, pink), and a viable electroweak
vacuum (4) (gray). See the text for further details.

It is important to note that the potential impact of the finite-temperature
effects and the EWPT on the scalar field evolution and relic abundance were
discussed in the original study of Ref. Piazza and Pospelov (2010). The effect
of the EWPT on scalar field misalignment in Higgs portal models was also
investigated in Ref. Arkani-Hamed  _et al._ (2021), though thermal effects
were not considered in that work. We believe our study elucidates the dynamics
of the scalar field during the radiation era by including both the thermal
misalignment and VEV misalignment effects, discerning the regions of parameter
space where each is relevant, and exploring the role of and sensitivity to
initial conditions. In particular, we have carefully derived estimates of the
scalar relic density, both by numerically solving its equation of motion
during the radiation era through the EWPT and by developing various
approximate solutions as appropriate for the region of parameter space,
cosmological epoch, and initial conditions under consideration. Both our
numerical and analytical results are shown in Fig. 1 and are in good
agreement.

This paper is organized as follows. In Sec. II we describe the super-
renormalizable Higgs portal model, including its zero temperature vacuum and
scalar spectrum. Sec. III describes the cosmology of the scalar field and
forms the bulk of the paper. It contains several subsections covering the
effective potential, EWPT, sources of misalignment and initial conditions,
scalar field equation of motion and evolution, and analytic relic density
estimates. Sec. IV compares our relic abundance predictions with existing and
proposed probes of the model. Finally, our conclusions are presented in Sec.V.

##  II Super-Renormalizable Higgs Portal Model

We consider the model of Ref. Piazza and Pospelov (2010) which contains a real
singlet scalar field $\phi$ that couples to the SM via the super-
renormalizable Higgs portal. The scalar potential to be considered is

$V=-\mu^{2}\,H^{\dagger}H+\lambda\,(H^{\dagger}H)^{2}+\frac{1}{2}m_{\phi}^{2}\,\phi^{2}+A\,\phi\,H^{\dagger}H\,,$
(1)

where the first two terms constitute the usual SM Higgs potential, $m_{\phi}$
is the scalar mass parameter, and $A$ is a dimensionful coupling of the scalar
to the Higgs. Note that we have omitted a potential linear term in $\phi$,
which can always be achieved by performing a field redefinition.

We will mainly be interested in feeble couplings such that $A\ll m_{\phi}$.
This implies that the light scalar mass is technically natural, with loops
generating $\delta m_{\phi}^{2}\sim A^{2}/(16\pi^{2})\log{\Lambda_{\rm UV}}\ll
m_{\phi}$. This is one of the attractive features of the super-renormalizable
portal. Furthermore, radiatively generated nonlinear potential terms such as
$\phi^{3}$ and $\phi^{4}$ are also small and will be neglected for the
remainder of this work.

###  II.1 Zero temperature vacuum and spectrum

We first consider the theory at zero temperature. The background Higgs field
is parameterized as $H^{T}=(0,\tfrac{1}{\sqrt{2}}h)$. The potential then reads

$\displaystyle
V_{0}(\phi,h)=-\frac{1}{2}\,\mu^{2}\,h^{2}+\frac{1}{4}\lambda\,h^{4}+\frac{1}{2}m_{\phi}^{2}\phi^{2}+\frac{1}{2}A\,\phi\,h^{2}\,.$
(2)

Minimizing the potential, we obtain the vacuum expectation values (VEVs) for
the scalar fields $\langle h\rangle=v$ and $\langle\phi\rangle=\phi_{0}$,
given by

$\displaystyle
v^{2}=\frac{\mu^{2}}{\lambda-A^{2}/2m_{\phi}^{2}}\,,~{}~{}~{}~{}~{}~{}~{}\phi_{0}=-\frac{Av^{2}}{2m_{\phi}^{2}}\,,$
(3)

with $v=246$ GeV. It can be seen from Eq. (3) that a viable electroweak vacuum
is obtained only for

$\displaystyle\frac{A^{2}}{2\,\lambda\,m_{\phi}^{2}}<1\,.$ (4)

To study the spectrum at zero temperature, we replace $h\rightarrow
v+\tilde{h}$ and $\phi\rightarrow\phi_{0}+\tilde{\phi}$ in the potential, Eq.
(2), where $\tilde{h}$ and $\tilde{\phi}$ represent the fluctuations about the
vacuum. Expanding the potential to quadratic order, we find that there is mass
mixing between the scalars,

$\displaystyle V$ $\displaystyle\supset\frac{1}{2}(2\lambda
v^{2})\tilde{h}^{2}+\frac{1}{2}m_{\phi}^{2}\,\tilde{\phi}^{2}+A\,v\,\tilde{h}\,\tilde{\phi}\,.$
(5)

We move to the physical basis by performing a rotation,

$\displaystyle\left[\begin{array}[]{c}\tilde{h}\\\
\tilde{\phi}\end{array}\right]=\left[\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{array}\right]\left[\begin{array}[]{c}h\\\
\phi\end{array}\right],$ (12)

where, in an abuse of notation, $h$ and $\phi$ represent the physical scalar
fluctuations in the mass basis. The mixing angle $\theta$ in (12) is given by

$\tan 2\theta=\frac{2Av}{2\lambda v^{2}-m_{\phi}^{2}}\,.$ (13)

The mass eigenvalues are

$M^{2}_{h,\phi}=\frac{1}{2}\left[2\lambda v^{2}+m_{\phi}^{2}\pm\sqrt{(2\lambda
v^{2}-m_{\phi}^{2})^{2}+4A^{2}v^{2}}\right]\,.$ (14)

As we will always be working in the regime $A^{2}\leq 2\lambda
m_{\phi}^{2}\ll\lambda v^{2}$, the approximate expressions for the mixing
angle and mass eigenvalues are

$\displaystyle\theta\sim\frac{A}{2\lambda
v}\simeq\frac{Av}{M_{h}^{2}}\,,~{}~{}~{}~{}~{}~{}M_{h}^{2}\simeq 2\lambda
v^{2}+\frac{A^{2}}{2\lambda}\,,~{}~{}~{}~{}~{}~{}M_{\phi}^{2}\simeq
m_{\phi}^{2}-\frac{A^{2}}{2\lambda}\,.$ (15)

Note that the correction to the Higgs mass is always negligible and thus
$M_{h}^{2}\simeq 2\lambda v^{2}\equiv m_{h}^{2}=(125~{}\rm GeV)^{2}$ as in the
SM. Furthermore, we will find that in most of the cosmologically interesting
parameter space, $A\ll m_{\phi}$, such that $M^{2}_{\phi}\simeq m^{2}_{\phi}$.
Note that very close to the boundary where the electroweak vacuum is viable
(see Eq. (4)), $A^{2}/(2\lambda m_{\phi}^{2})\lesssim 1$, and only in this
very small region of parameter space is the physical scalar mass is
substantially modified from $m_{\phi}^{2}$.

The electroweak vacuum condition, Eq. (4) (gray shaded region in Fig. 1),
could be relaxed in a scenario with additional scalar potential terms.
However, in this case the light scalar mass $M_{\phi}$ would require fine-
tuning, as can be seen from Eq. (15) above. In this work, we will restrict our
investigation of the scalar cosmology to parameters satisfying Eq. (4).

##  III Cosmology

In this section, we describe the cosmological evolution of the scalar $\phi$
before and through the EWPT until it oscillates about its late time potential
minimum and behaves as DM. We will assume that the universe reheats to a
temperature much larger than the electroweak scale $(T_{\rm RH}\gg v)$. We
begin by describing the contribution to the scalar potential from the thermal
bath. We will then describe the EWPT experienced by the Higgs, which can
significantly impact the dynamics of $\phi$ in certain regions of parameter
space. We describe the sources of misalignment and detail the initial
conditions assumed for the scalar at the end of inflation. We then describe
the final equation of motion that governs the scalar evolution and numerically
estimate the relic $\phi$ abundance. We also derive approximate analytic
estimates of the relic density and compare them with our numerical results.

###  III.1 Effective Potential

Our analysis begins in the radiation-dominated era during which the Hubble
parameter is given by $H=1/2t=\gamma T^{2}/M_{\rm pl}$, where $t$ is the
cosmic time, $T$ is the temperature of the SM thermal bath, $M_{\rm
pl}=2.43\times 10^{18}$ GeV is the reduced Planck mass, and
$\gamma(T)\equiv\sqrt{\pi^{2}g_{*}(T)/90}$ with $g_{*(S)}$ the effective
number of relativistic (entropy) degrees of freedom. During this epoch, the
effective potential of the scalar fields is given by

$\displaystyle V_{\rm eff}(\phi,h,T)$ $\displaystyle=V_{0}(\phi,h)+V_{\rm
CW}(\phi,h)+V_{T}(\phi,h,T)\,.$ (16)

The first piece is the tree-level potential given in Eq. (1). The second term
represents the zero-temperature correction to the effective potential, i.e.,
the Coleman-Weinberg potential Coleman and Weinberg (1973), which we will
comment on shortly. The final term in (16) is the finite-temperature
correction associated with the thermal free energy density of the SM particles
in the plasma Dolan and Jackiw (1974); Weinberg (1974). The expression for the
finite-temperature potential in Landau gauge is Quiros (1999)

$\displaystyle V_{T}(\phi,h,T)$
$\displaystyle\supset\frac{1}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{h}^{2}(\phi,h,T)}{T^{2}}\right]+\frac{3}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{\chi}^{2}(\phi,h,T)}{T^{2}}\right]+\frac{4}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{W_{T}}^{2}(h)}{T^{2}}\right]$
$\displaystyle+\frac{2}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{Z_{T}}^{2}(h)}{T^{2}}\right]+\frac{2}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{W_{L}}^{2}(h,T)}{T^{2}}\right]+\frac{1}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{Z_{L}}^{2}(h,T)}{T^{2}}\right]$
$\displaystyle+\frac{1}{2\pi^{2}}T^{4}J_{B}\left[\frac{m_{A_{L}}^{2}(h,T)}{T^{2}}\right]-\frac{12}{2\pi^{2}}T^{4}J_{F}\left[\frac{m_{t}^{2}(h)}{T^{2}}\right]\,,$
(17)

where the $J_{B,F}$ functions are defined as

$\displaystyle J_{B,F}(w^{2})$
$\displaystyle=\int_{0}^{\infty}\\!dx\,x^{2}\,\log\bigg{[}1\mp\exp\left({-\sqrt{x^{2}+w^{2}}}\right)\bigg{]}\,.$
(18)

The thermal squared masses, $m_{i}^{2}(\phi,h,T)$, entering in the arguments
of these functions depend in general on both the background field fields
$\phi$ and $h$ and the temperature. To obtain a reliable perturbative
expansion near the EWPT, we follow the prescription for daisy resummation of
Ref. Parwani (1992), writing the thermal squared mass as sum of a tree-level
term field-dependent term, $m^{2}_{0,i}(\phi,h)$, and a temperature-dependent
self-energy term $\Pi_{i}(T)$,

$\displaystyle m_{i}^{2}(\phi,h,T)=m^{2}_{0,i}(\phi,h)+\Pi_{i}(T)\,.$ (19)

In particular, for the self-energies we retain the leading contributions in
the high temperature expansion. The explicit expressions for the thermal
masses are collected in App. A. The subscripts
$(h,\chi,W_{T},Z_{T},W_{L},Z_{L},A_{L},t)$ in these functions refer to the
Higgs, Nambu-Goldstones, transverse and longitudinal gauge bosons, and top
quark. The prefactors in front of the $J_{B,F}$ functions account for the
degrees of freedom of the corresponding field.

It is important to emphasize that $\phi$ itself is not in thermal equilibrium
for the feeble values of Higgs portal coupling $A$ considered in this work,
and as such there is no term associated with $\phi$ in the finite temperature
potential (III.1). Rather, the finite temperature potential is due to the SM
degrees of freedom in thermal equilibrium. It is a function of the background
value for $\phi$ via the dependence of the squared mass parameters on $\phi$.

For simplicity, we will neglect the Coleman-Weinberg (CW) contribution $V_{\rm
CW}(\phi,h)$ to the effective potential (16) in our numerical analysis below.
The correction to the $\phi$ potential from $V_{\rm CW}$ is negligible in the
viable region of parameter space (Eq. (4)) due to the feeble coupling $A$, as
already alluded to at the beginning of Sec. II. Furthermore, as we will see
below, for large scalar masses the misalignment is generated at high
temperatures (thermal misalignment) where $V_{T}(\phi,h,T)$ dominates and
$V_{\rm CW}$ is negligible. For lower scalar masses, the misalignment is a
result of the induced scalar VEV triggered by the EWPT. Even in this case,
$V_{\rm CW}$ does not qualitatively alter the behavior of the scalar potential
and the nature of the EWPT, and will only lead to relatively small numerical
differences in our results.

Note that Eq. (III.1) is appropriate for temperatures above the QCD phase
transition. For lower temperatures, an alternative description in terms of
light hadronic states would be required. However, we will see the onset of
scalar oscillations in the parameter space we consider occurs at temperatures
well above the QCD scale, so that a description at lower temperatures is not
required to calculate the relic abundance. For similar reasons, the quarks
lighter than the top quark as well as leptons will not influence our results.

We will find it convenient to work with the following set of dimensionless
parameters:

$y\equiv\frac{T}{\mu}\,,\quad~{}~{}\hat{\phi}\equiv\frac{\phi}{M_{\rm
pl}}\,,\quad~{}~{}\hat{h}\equiv\frac{h}{\mu}\,,\quad~{}~{}\kappa\equiv\frac{m_{\phi}M_{\rm
pl}}{\mu^{2}}\,,\quad~{}~{}\beta\equiv\frac{AM_{\rm pl}}{\mu^{2}}\,.~{}~{}$
(20)

In particular, the parameter $y$ is a proxy for the temperature, with
$y\sim{\cal O}(1)$ for temperatures near the electroweak scale. Furthermore,
$\hat{\phi}$ and $\hat{h}$ are simply the ratios of the field variables to
their respective characteristic scales. Finally, the scalar mass $m_{\phi}$
and Higgs portal coupling $A$ are expressed in terms of the dimensionless
variables $\kappa$ and $\beta$, respectively.

Combining the $V_{0}$ and $V_{1}^{T}$ in Eqs. (1,III.1) and writing the
effective potential in terms of dimensionless quantities ($\hat{V}_{\rm
eff}=V_{\rm eff}/\mu^{4}$), we obtain

$\displaystyle\hat{V}_{\rm eff}$
$\displaystyle=-\frac{1}{2}\hat{h}^{2}(1-\beta\hat{\phi})+\frac{1}{4}\lambda\hat{h}^{4}+\frac{1}{2}\kappa^{2}\hat{\phi}^{2}+\frac{y^{4}}{2\pi^{2}}\bigg{\\{}J_{B}[\eta_{h}(\hat{\phi},\hat{h},\hat{y})]+3J_{B}[\eta_{\chi}(\hat{\phi},\hat{h},\hat{y})]+4J_{B}[\eta_{W_{T}}(\hat{h})]$
$\displaystyle+2J_{B}[\eta_{Z_{T}}(\hat{h})]+2J_{B}[\eta_{W_{L}}(\hat{h},\hat{y})]+J_{B}[\eta_{Z_{L}}(\hat{h},\hat{y})]+J_{B}[\eta_{A_{L}}(\hat{h},\hat{y})]-12J_{F}[\eta_{t}(\hat{h},\hat{y})]\bigg{\\}}\,,$
(21)

where we have defined the dimensionless arguments of the $J_{B,F}$ functions,
$\eta_{i}(\hat{\phi},\hat{h},y)\equiv m^{2}_{i}(\phi,h,T)/T^{2}$. Explicit
expressions for these quantities are provided in App. A.

###  III.2 Higgs field and Electroweak Phase Transition

We first discuss the evolution of the Higgs field $h$, which serves as the
order parameter for the EWPT. In the SM, the EWPT is a smooth crossover
Kajantie  _et al._ (1996) characterized by the critical temperature
$T_{c}\sim{\cal O}(v)$ ($y_{c}\equiv T_{c}/\mu\sim{\cal O}(1)$). The feeble
portal coupling between $\phi$ and $h$ will not alter the nature of the phase
transition, at least in the $\beta\hat{\phi}\ll 1$ regime of primary relevance
for this work. The Higgs field is assumed to track the potential minimum
throughout the phase transition, starting at $h=0$ at high temperatures $(y\gg
y_{c})$, then taking nonzero values for $y<y_{c}$, and finally settling at
$h=v$ at low temperatures $(y\ll 1)$.

Thus, for temperatures below $y_{c}$, the evolution of $\hat{h}(\hat{\phi},y)$
is determined by the minimization condition $\partial\hat{V}_{\rm
eff}/\partial\hat{h}=0$, with $\hat{V}_{\rm eff}$ given in Eq. (III.1):

$\displaystyle 0$
$\displaystyle=\lambda\hat{h}^{2}-(1-\beta\hat{\phi})+\frac{y^{2}}{2\pi^{2}}\bigg{\\{}6\lambda(J^{\prime}_{B}[\eta_{h}]+J^{\prime}_{B}[\eta_{\chi}])+g^{2}\left(2J^{\prime}_{B}[\eta_{W_{T}}]+J^{\prime}_{B}[\eta_{W_{L}}]\right)$
$\displaystyle+(g^{2}+g^{\prime
2})J^{\prime}_{B}[\eta_{Z_{T}}]+2y^{2}\frac{\partial\eta_{Z_{L}}}{\partial\hat{h}^{2}}J^{\prime}_{B}[\eta_{Z_{L}}]+2y^{2}\frac{\partial\eta_{A_{L}}}{\partial\hat{h}^{2}}J^{\prime}_{B}[\eta_{A_{L}}]-12y_{t}^{2}J^{\prime}_{F}[\eta_{t}]\bigg{\\}}\,.$
(22)

Here $g$, $g^{\prime}$, and $y_{t}$ are the $SU(2)_{L}$ gauge coupling,
$U(1)_{Y}$ gauge coupling, and top Yukawa coupling, respectively, and the
respective $\eta_{i}$ are given in App. A. Note that the $\eta_{i}$ are
functions of $(\hat{\phi},\hat{h},y)$. Our procedure is to solve Eq. (III.2)
for $\hat{h}(\hat{\phi},y)$ and then use this solution in equation of motion
of $\hat{\phi}$, to be discussed below in Sec. III.4. In principle, Eq.
(III.2) can be solved numerically. However, since $\beta\hat{\phi}\ll 1$ for
essentially all of our parameter space, in practice we will treat
$\beta\hat{\phi}$ as a perturbation. We therefore define $\hat{h}_{0}(y)$ as
the solution to Eq. (III.2) with $\beta\hat{\phi}=0$. Including the
$\beta\hat{\phi}$ dependence will change the evolution and $y_{c}$ slightly
($\propto\mathcal{O}(\beta\hat{\phi})$). The evolution of $\hat{h}_{0}(y)$ is
shown in Fig. 2.

Figure 2: Dimensionless Higgs profile $\hat{h}_{0}=\hat{h}(\hat{\phi}=0,y)$ as
a function of $y=T/\mu$. The critical temperature corresponding to the EWPT
can be seen at $y_{c}\simeq 1.6$ .

###  III.3 Sources of Misalignment and Initial Conditions

Before we study the evolution of the scalar field $\phi$ in detail, it is
worth discussing the various sources of misalignment that contribute to the
late-time $\phi$ oscillation amplitude and associated relic density. To gain
insight, it is helpful to examine the minimum of the effective potential with
respect to $\hat{\phi}$ as a function of temperature, $\hat{\phi}_{\rm
min}(y)$, which can be obtained from Eq. (III.1) by setting
$\partial\hat{V}_{\rm eff}/\partial\hat{\phi}=0$. The solution can be written
as

$\displaystyle\hat{\phi}_{\rm
min}=-\frac{\beta}{2\kappa^{2}}\left[\hat{h}^{2}+\frac{y^{2}}{\pi^{2}}\left(J^{\prime}_{B}[\eta_{h}]+3J^{\prime}_{B}[\eta_{\chi}]\right)\right]\,.$
(23)

Since $\beta\hat{\phi}\ll 1$ for the bulk of our parameter space, it is
typically a good approximation to neglect this dependence in $\eta_{h,\chi}$
and $\hat{h}$ which enter in the r.h.s. of Eq. (23). In Fig. 3 we show the
variation of $\hat{\phi}_{\rm min}(y)$ with $y=T/\mu$ for a representative
benchmark model. We see that at high temperatures, $y\gg y_{c}$, the minimum
is located at large (negative) scalar field values, reflecting the dominance
of the second term in Eq. (23). As the temperature drops, the minimum
decreases as $|\hat{\phi}_{\rm min}|\propto y^{2}$ until the EWPT at $y_{c}$.
At this point, the Higgs field turns on and the $\hat{\phi}$ potential minimum
rapidly transitions to its zero temperature value, Eq. (3), or in terms of the
dimensionless variables,

$\hat{\phi}_{0}=-\frac{\beta}{2\lambda\kappa^{2}-\beta^{2}}\simeq-\frac{\beta}{2\lambda\kappa^{2}}\,.$
(24)

The qualitative behavior of $\hat{\phi}_{\rm min}$ is similar to that in Fig.
3 for other parameter choices.

Figure 3: Evolution of the temperature-dependent scalar field minimum
$|\hat{\phi}_{\rm min}(y)|$ as a function of $y=T/\mu$ for the benchmark
$\kappa=10^{3}$, $\beta=1$. The $\phi$ minimum is controlled by the finite-
temperature contribution to the effective potential for $y>y_{c}$, while
instead for $y<y_{c}$ it transitions quickly towards its zero temperature
value induced by the the Higgs VEV, Eq. (3). We have used
$\hat{h}_{0}=\hat{h}(\hat{\phi}=0,y)$ as a function of $y=T/\mu$.

As we now discuss, there are three qualitatively distinct contributions to the
misalignment of $\hat{\phi}$:

  * •

Misalignment due to initial conditions - This contribution is given by the
value $\hat{\phi}_{i}$ the scalar field takes at the end of inflation.

  * •

Thermal misalignment - At temperatures much larger than the temperature of
EWPT, $y\gg y_{c}$ , the scalar experiences a finite-temperature contribution
to the effective potential due to its feeble coupling to the thermalized Higgs
degrees of freedom given by Eq. (III.1). This causes the scalar to roll
towards the high-temperature minimum set by Eq. (23) and illustrated in Fig.
3, dynamically generating misalignment. We refer to this contribution as
thermal misalignment, $\hat{\phi}_{T}$, and we will show below in Sec. III.5.1
that

$\hat{\phi}_{T}\propto\frac{\beta}{\kappa}\,.$ (25)

  * •

VEV misalignment - Below the temperature of the EWPT, $y<y_{c}$ , the Higgs
field switches on and induces a rapid shift in the minimum of the scalar field
$\hat{\phi}$ (see Fig. 3) to its zero temperature VEV (24). This dynamically
generates and additional source of misalignment, $\hat{\phi}_{V}$, which we
refer to as VEV misalignment.

The final misalignment is in general influenced by all of the above
contributions. Note that both thermal misalignment and VEV misalignment are
dynamical processes occurring during the radiation era. As we will see, these
sources of misalignment depend on the DM model parameters, i.e., the mass
$m_{\phi}$ (or equivalently $\kappa$) and the coupling $A$ (or equivalently
$\beta$). On the other hand, the initial value $\hat{\phi}_{i}$ after
reheating depends on the detailed dynamics during inflation and can in
principle be arbitrary. An important question then is to what extent the relic
abundance prediction depends on the assumed initial condition. If the relic
density prediction is insensitive to the initial conditions and only depends
on the DM model parameters, this opens up the exciting prospect of confronting
the cosmological production mechanism with experiments, since the signatures
and predictions for the latter are determined by the same model parameters
(mass and coupling).

As we will show in detail in Sec. III.5.1, for larger scalar masses,
$m_{\phi}\gtrsim$ few meV ($\kappa\gtrsim 10^{3}$), thermal misalignment
dominates over VEV misalignment, $\hat{\phi}_{T}\gg\hat{\phi}_{V}$. Then
provided $\hat{\phi}_{i}\ll\hat{\phi}_{T}$, which can be satisfied for a broad
range of conceivable initial field values given Eq. (25), the thermal
misalignment generated during the radiation era will overwhelm (or “erase”)
the contribution from the initial field value. If this condition holds, the
late-time oscillation amplitude and resulting scalar relic abundance is
controlled entirely by the DM model parameters. We can then define a relic
density target corresponding to a line in the $m_{\phi}-A$ plane yielding a DM
abundance in accord with the measured value, $\Omega_{\rm DM}\simeq 0.26$
Aghanim  _et al._ (2020). Modulo fine tuning of initial conditions, parameters
above the relic density target line will lead to even larger scalar
oscillation amplitudes generated by thermal misalignment and would thus
overclose the universe. These features, such as correlations between DM relic
density and its mass and coupling, as well as insensitivity to initial
conditions and UV physics, are reminiscent of similar attractive aspects found
in models of weakly-interacting-massive particle (WIMP) DM.

For lower masses, $m_{\phi}\lesssim 10^{-5}$ eV ($\kappa\lesssim 1$), VEV
misalignment dominates over thermal misalignment, $\phi_{V}\gg\phi_{T}$. As
will be clearly shown in Sec. III.5, the insensitivity of the scalar relic
density to the initial field value $\hat{\phi}_{i}$ does not carry over to
this lower mass regime.

Below we will detail two choices of initial conditions that will be used in
our numerical and analytic estimates in the coming sections. These two choices
will allow us to demonstrate the insensitivity of the relic density to initial
conditions for scalars with relatively large masses where thermal misalignment
dominates, as well as study the impact of different initial conditions for
scalars with relatively low masses where VEV misalignment is important.

  * •

$\phi_{i}=\phi_{0}$:

In this case, the scalar is initially at its zero-temperature VEV, $\phi_{0}$,
given by Eq. (3). This initial condition is naturally realized through a
sufficiently long enough period of inflation with a low enough inflationary
Hubble parameter, $H_{I}\ll v$. In this case, electroweak symmetry is broken
during inflation, the Higgs is close to its VEV $v$, and the $\phi$ potential
minimum is approximately at $\phi_{0}$. Given a long enough period of
inflation the distribution of $\phi$ will relax to the average
$\langle\phi\rangle|_{H_{I}}=\phi_{i}=\phi_{0}$. This relaxation will take
order $\mathcal{N}\sim H^{2}_{I}/m^{2}_{\phi}$ $e$-folds, and the variance of
$\phi$ will be given by $\sigma_{\phi_{i}}=3H^{4}_{I}/(8\pi m^{2}_{\phi})$.
Note that $H_{I}$ cannot be arbitrarily low if thermal misalignment is to be
important, which requires the reheat temperature to satisfy $T_{\rm RH}\gg v$,
though this condition can easily be met. We also require
$\sigma_{\phi_{i}}\ll\phi_{T}$, $(\phi_{V})$ for the thermal (VEV)
misalignment to be important. We note that a low inflationary Hubble scale
naturally suppresses scalar fluctuations (see, e.g., Refs. Graham and Scherlis
(2018); Takahashi  _et al._ (2018); Tenkanen (2019)), thus easing otherwise
stringent CMB constraints on isocurvature perturbations Akrami  _et al._
(2020).

As emphasized above, for heavier scalars the displacement generated by thermal
misalignment overwhelms the initial field value assumed here. On the other
hand, for lighter scalars thermal misalignment is negligible and VEV
misalignment operates, though we must carefully consider the impact of the
assumed initial condition. The scalar is held fixed by Hubble friction for
some time near its initial value, $\hat{\phi_{i}}\sim\hat{\phi_{0}}$. Then at
temperatures just above the EWPT, the $\hat{\phi}$ potential minimum adjusts
to field values that are much smaller than its initial value $\hat{\phi}_{0}$,
as can be seen in in Fig. 3. At this stage, the scalar field slow-rolls away
from its initial position, generating misalignment. After the EWPT, the
location of the $\hat{\phi}$ minimum quickly transitions to its zero-
temperature VEV $\hat{\phi}_{0}$, and eventually oscillations commence after
the expansion rate drops below the scalar mass. It will be shown that the
parametric dependence of the VEV misalignment for this inital condition is

$\displaystyle\hat{\phi}_{V}\propto\beta,~{}~{}~{}~{}~{}~{}~{}~{}(\hat{\phi}_{i}=\hat{\phi}_{0})\,.$
(26)

A more thorough treatment of this dynamics will be presented in Sec. III.5.3.

  * •

$\phi_{i}=0$:

This initial condition is chosen as a representative example of the general
situation where $|\phi_{i}|$ is vastly different than $|\phi_{0}|$. For larger
scalar masses, thermal misalignment again dominates the evolution and the
scalar relic abundance is insensitive to the assumed initial field value. In
contrast, for smaller scalar masses, VEV misalignment dominates over thermal
misalignment, although we must again understand the role of the initial
condition. In fact, the evolution of the scalar field in this case is quite
simple. Initially, the scalar is held at its initial location at the origin by
Hubble friction. The EWPT triggers a shift in the $\phi$ potential minimum to
its zero-temperature value $\phi_{0}$, generating VEV misalignment of
parametric size (see Sec. III.5.3):

$\displaystyle\hat{\phi}_{V}\simeq-\frac{\beta^{2}}{2\kappa^{2}\lambda}~{}~{}~{}~{}~{}~{}~{}~{}(\hat{\phi}_{i}=0)\,.$
(27)

At some later time the Hubble rate falls below the scalar mass and
oscillations begin.

It can be seen in Fig. 1 that the relic density lines in the lower scalar mass
region are strongly sensitive to the initial field value, while instead those
in the higher scalar mass range are the same for both assumed initial
conditions.

###  III.4 Scalar field dynamics

We are now ready to discuss the evolution of the scalar field $\phi$, which is
governed by the following equation of motion during the radiation-dominated
era

$\displaystyle\ddot{\phi}+3H\dot{\phi}+\frac{\partial V_{\rm
eff}}{\partial\phi}=0\,,$ (28)

where again $V_{\rm eff}=\mu^{4}\,\hat{V}_{\rm eff}$ is given in Eqs.
(16,III.1) and the Hubble parameter $H$ has been defined at the beginning of
Sec. III.1.

For simplicity, we will ignore the small time variation of $g_{*}$ during the
initial phase of the $\phi$ evolution until $\phi$ oscillations begin, fixing
it to the high temperature SM value of $g_{*}^{\rm SM}=106.75$
($\gamma=\sqrt{\pi^{2}g_{*}/90}\simeq 3.4$). For large scalar masses,
oscillations begin well before the EWPT, and as such $g_{*}=g_{*}^{\rm SM}$ is
exactly correct. Instead, for smaller masses, oscillations begin below the
EWPT and $g_{*}$ starts to decreases as $t,W,Z,h$ fall out of thermal
equilibrium, giving for example $g_{*}\sim 80$ at $y\sim 0.01$. Thus for
smaller masses our constant $g_{*}=g_{*}^{\rm SM}$ approximation will result
in a $\mathcal{O}(10\%)$ error in the relic abundance prediction.

Let us also briefly discuss a potential additional source of scalar field
damping. As the scalar field evolves, the effective mass of the Higgs
particles in the bath changes though their momenta does not, causing a small
deviation in the Higgs phase space distribution from its equilibrium value.
This can manifest as an additional effective friction term in the scalar
equation of motion; see e.g., Ref. Yokoyama and Linde (1999). However, the
Higgs quickly relaxes toward equilibrium through fast number changing
processes, such that this source of damping is much smaller than the other
terms in Eq. (28). See also Ref. Bödeker and Nienaber (2022) for discussion.

We can re-express the $\phi$ equation of motion (28) in terms of the
dimensionless variables defined in Eq. (20), with independent variable of time
$t$ traded for $y=T/\mu$. We obtain

$\displaystyle\hat{\phi}^{\prime\prime}+\frac{1}{\gamma^{2}y^{6}}\left[\kappa^{2}\hat{\phi}+\frac{\beta\hat{h}^{2}}{2}+\frac{\beta
y^{2}}{2\pi^{2}}\left(J^{\prime}_{B}[\eta_{h}]+3J^{\prime}_{B}[\eta_{\chi}]\right)\right]=0\,.$
(29)

Here, $\hat{h}=\hat{h}(\hat{\phi},y)$ is the solution to Eq. (III.2), as
discussed in detail in Sec. III.2. We then solve the $\hat{\phi}$ equation of
motion (29) subject to the initial conditions as described in the previous
section. These are defined at $y_{i}=T_{i}/\mu$, where $y_{i}\gg 1$
corresponds to some suitably early time well before the scalar starts to
oscillate. The scalar field is evolved until $y_{f}=T_{f}/\mu$ corresponding
to some suitable time after it begins oscillating and redshifts as matter.
From this time on, the $\phi$ comoving energy density, $\rho_{\phi}/s$, where
$\rho_{\phi}$ is the $\phi$ energy density and
$s(T)=(2\pi^{2}/45)g_{*S}(T)T^{3}$ is the total entropy density of the plasma,
is conserved. Using this fact, we then arrive at our estimate for the $\phi$
energy density in the present epoch,

$\displaystyle\rho_{\phi,0}=\frac{g_{*S}(y_{0})}{g_{*S}(y_{f})}\left(\frac{y_{0}}{y_{f}}\right)^{3}\rho_{\phi}(y_{f})\,,$
(30)

where $y_{0}=T_{0}/\mu$ with $T_{0}\simeq 2.7$ K. We use
$g_{*S}(y_{f})=106.75$ and $g_{*S}(y_{0})\simeq 3.91$ in our estimate.
Expressing the relic density in terms of the density parameter
$\Omega_{\phi,0}=\rho_{\phi,0}/\rho_{c,0}$, with $\rho_{c,0}=3M_{\rm
pl}^{2}H_{0}^{2}$ the critical density today, we show the parameters leading
to the observed DM relic abundance, $\Omega_{\phi,0}=\Omega_{\rm DM}\simeq
0.26$ Aghanim  _et al._ (2020) in Fig. 1.

Before discussing these results and comparing with the various experimental
probes of the model, it is valuable to gain analytical insight into the $\phi$
evolution and resulting relic density estimate, which we explore in the next
section.

###  III.5 Relic Density Estimation

As the universe expands, the Hubble parameter decreases until it eventually
falls below the effective $\phi$ mass, marking the onset of scalar
oscillations. This can be quantified by the condition $[3H(y_{\rm
osc})]^{2}=m_{\phi}^{2}(y_{\rm osc})\simeq m_{\phi}^{2}$, which gives the
oscillation temperature as

$\displaystyle y_{\rm osc}=\frac{T_{\rm
osc}}{\mu}=\sqrt{\frac{\kappa}{3\gamma}}\,.$ (31)

Here we have used the fact that the effective squared scalar mass is dominated
by the tree-level contribution $m_{\phi}^{2}$ over the entire parameter space
we consider and for all temperatures. Several epochs after the oscillations
begin, we can safely calculate the energy density stored in the scalar field
and then redshift it to present times, which will yield the $\phi$ relic
density given by Eq. (30). Typically, once the scalar starts to oscillate, it
behaves like DM with an oscillation amplitude set by

$\displaystyle\phi_{\rm osc}\equiv\phi(y_{\rm osc})\,,$ (32)

i.e., the scalar amplitude at the beginning of oscillations. To determine
$\phi_{\rm osc}$ we must study the evolution of the scalar field during the
radiation era up until the onset of oscillations. As discussed at length in
Sec. III.3, there are additional dynamical sources of misalignment present in
the Higgs portal model and contribute to the oscillation amplitude. For large
scalar masses, thermal misalignment provides the main contribution, while for
small scalar masses VEV misalignment (along with initial conditions) is
dominant. In the intermediate mass range, both effects are important. In the
following, we will therefore divide the parameter space into distinct regions
characterized by which sources of misalignment dominate and then develop a
suitable analytic estimate for the $\phi$ evolution and relic abundance. These
estimates will be compared with our relic density prediction obtained from the
exact numerical solution of the $\phi$ equation of motion (28), as described
in the previous section.

####  III.5.1 Region I $(\kappa\gtrsim 10^{3},~{}m_{\phi}\gtrsim 3\times
10^{-3}\,{\rm eV})$

Region I is defined to be the region of large scalar masses $m_{\phi}$ (large
$\kappa$). In this region, the the amplitude of oscillations at $y_{c}$ due to
thermal misalignment dominates over the “kick” imparted to the scalar
resulting from the EWPT (VEV misalignmnet). To understand this point, we first
need to estimate the thermal misalignment generated at high temperatures by
considering Eq. (29). For $y\gg y_{\rm osc}\gg y_{c}$,
$\hat{h}(\hat{\phi},y)=0$ and we can ignore the $\kappa^{2}\hat{\phi}$ term
since it is subdominant to the finite-temperature contribution to the equation
of motion. Hence, in this regime, Eq. (29) takes the approximate form

$\displaystyle\hat{\phi}^{\prime\prime}(y)+\frac{\beta}{2\pi^{2}\gamma^{2}y^{4}}\left(J^{\prime}_{B}[\eta_{h}]+3(J^{\prime}_{B}[\eta_{\chi}]\right)$
$\displaystyle=0,$
$\displaystyle\Longrightarrow~{}~{}~{}~{}\hat{\phi}^{\prime\prime}(y)+\frac{\beta}{\pi^{2}\gamma^{2}y^{4}}$
$\displaystyle=0\,,$ (33)

where in the second line we have used the fact that
$J^{\prime}_{B}[\eta_{h}]+3J^{\prime}_{B}[\eta_{\chi}]\approx 2$ for $y\gg
y_{c}$ (see App. A for definitions of $\eta_{h},\eta_{\chi}$). The solution to
Eq. (III.5.1), assuming negligible initial field velocity, is

$\displaystyle\hat{\phi}(y)=-\frac{\beta}{6\pi^{2}\gamma^{2}y^{2}}+\hat{\phi}_{i}\,,$
(34)

where the first term represents the thermal misalignment,
$\hat{\phi}_{T}(y)=-\beta/(6\pi^{2}\gamma^{2}y^{2})$, while the second is the
initial condition $\hat{\phi_{i}}$. As the temperature drops below $y_{\rm
osc}=\sqrt{\kappa/3\gamma}$, the mass term $(\kappa^{2}\hat{\phi})$ in Eq.
(29) begins to dominate over the thermal contribution and the scalar starts to
oscillate and redshift like DM. Plugging $y=y_{\rm osc}$ into Eq. (34), we
find the scalar amplitude at the onset of oscillations is given by

$\displaystyle\hat{\phi}_{\rm osc}=\hat{\phi}(y_{\rm
osc})=-\frac{\beta}{2\pi^{2}\gamma\kappa}+\hat{\phi}_{i}\,.$ (35)

For the first initial condition
$\hat{\phi}_{i}=\hat{\phi}_{0}\simeq-\beta/(2\lambda\kappa^{2})$ motivated in
Sec. III.3, we observe that thermal misalignment $\hat{\phi}_{T}$ dominates
over $\hat{\phi}_{i}=\hat{\phi}_{0}$ for
$\kappa\gtrsim\pi^{2}\gamma/\lambda\approx 300$. Thus, our definition of
Region I, $\kappa>10^{3}$, satisfies this criterion.

In Fig. 4 we display the numerical evolution of $\hat{\phi}(y)$ relative to
its temperature-dependent Higgs-induced minimum,
$\delta\phi(y)\equiv\hat{\phi}(y)-[-\beta\hat{h}_{0}^{2}(y)/(2\kappa^{2})]$,
for two benchmark models in Region I and the initial conditions
$\hat{\phi}_{i}=\hat{\phi}_{0}$ and $\hat{\phi}_{i}=0$ discussed in Sec.
III.3. This quantity provides a convenient visualization of both the thermal
misalignment generated at high temperatures as well as the negligible impact
of VEV misalignment for $y<y_{c}$. The temperature-dependent Higgs-induced
minimum, $[-\beta\hat{h}_{0}^{2}(y)/(2\kappa^{2})]$, is simply the first term
in the full temperature-dependent minimum $\hat{\phi}_{\rm min}$ in Eq. (23),
which dominates the expression at low temperatures, i.e.,

$\displaystyle\hat{\phi}_{\rm
min}(y)\simeq-\frac{\beta\hat{h}_{0}^{2}(y)}{2\kappa^{2}}\,~{}~{}~{}~{}~{}~{}~{}~{}(y\lesssim
y_{c}).$ (36)

This is also clear from Fig. 3. The dominance of thermal misalignment and the
independence of the the oscillation amplitude on the initial conditions are
clearly seen in Fig. 4.

Figure 4: Scalar field evolution $\hat{\phi}(y)$ relative to its temperature-
dependent Higgs-induced minimum,
$\delta\phi(y)\equiv\hat{\phi}(y)-[-\beta\hat{h}_{0}^{2}(y)/(2\kappa^{2})]$
for two benchmark models in Region I: $\beta=10^{-2},\kappa=10^{4}$ (blue) and
$\beta=10^{-3},\kappa=10^{4}$ (orange). The black dotted lines show the
corresponding approximate initial thermal misalignment trajectories of Eq.
(35). The $\beta^{2}$ scaling due to thermal misalignment is clearly observed.
The evolution is shown for the two choices of initial conditions discussed in
Sec. III.3: $\hat{\phi}_{i}=\hat{\phi}_{0}$ (solid) and $\hat{\phi}=0$
(dashed). It is evident that the late time $\hat{\phi}$ evolution is
independent of the initial value of the scalar field for these choices.

We will now justify the condition $\kappa\gtrsim 10^{3}$ used to define Region
I. For $\kappa\gtrsim 10^{3}$, oscillations of the scalar field around its
high temperature minimum begin at $y_{\rm osc}\gg y_{c}$. The amplitude of the
oscillations then decreases with temperature as $y^{3/2}$ (i.e., $\phi$
redshifts like matter) until $y=y_{c}$. At $y_{c}$ the Higgs field transitions
nearly instantaneously, causing a shift in the minimum of $\hat{\phi}$ from
the its minimum at $y>y_{c}$ to $\hat{\phi}_{0}$. If the change in the
$\hat{\phi}$ minimum near $y=y_{c}$ within one oscillation time period is
larger than the oscillation amplitude, then the relic density is controlled by
VEV misalignment. Conversely, if the oscillation amplitude at $y=y_{c}$ is
much larger than the shift in the $\hat{\phi}$ minimum due to the EWPT, then
thermal misalignment dominates. For large $\kappa$ (large $m_{\phi}$), the
latter scenario of thermal misalignment applies. In what follows, we will
derive a lower bound on $\kappa$ above which the thermal misalignment
dominates (defined to be Region I) and then develop an analytical
approximation for the $\phi$ relic density.

At $y_{c}$, the minimum of $\hat{\phi}$ transitions from being governed by the
thermal contribution to the potential to being controlled by the Higgs VEV
$\hat{h}$ (see Fig. 3). In this region, the $\hat{\phi}$ minimum swiftly
changes according $\hat{\phi}_{\rm min}$ given in Eq. (36). This rapid shift
in $\hat{\phi}_{\rm min}$ provides a “kick” by effectively changing the
amplitude of the oscillation. The magnitude of this kick is the variation of
$\hat{\phi}_{\rm min}$ within one half oscillation of the scalar field i.e.
within $\Delta t=\pi/m_{\phi}$. It can be shown that this time period
corresponds to $\Delta y=-y_{c}^{3}\gamma\pi/\kappa$. The change in the
$\hat{\phi}$ minimum is then given by

$\displaystyle\Delta\hat{\phi}_{\rm min}$
$\displaystyle\simeq\frac{\partial\hat{\phi}_{\rm min}}{\partial y}\Delta
y=-\frac{\beta}{2\kappa^{2}}\frac{\partial\hat{h}^{2}}{\partial y}\Delta y$
$\displaystyle\simeq-\frac{\pi\gamma y_{c}^{4}\beta}{2\lambda\kappa^{3}},$
(37)

where in the last line we have used $\partial\hat{h}^{2}/\partial y\simeq-
y_{c}/\lambda$ (this can be seen from the evolution equation for $\hat{h}^{2}$
in Eq. (III.2)). This is to be compared to the oscillation amplitude near
$y_{c}$ due to thermal misalignment, which is given by

$\displaystyle\hat{\phi}_{\rm osc}(y_{c})=\hat{\phi}(y_{\rm
osc})\left(\frac{y_{c}}{y_{\rm
osc}}\right)^{3/2}=-\frac{(\sqrt{3}y_{c})^{3/2}\beta}{2\pi^{2}\gamma^{1/4}\kappa^{7/4}}\,,$
(38)

where we have used Eq. (35) and redshifted it to $y_{c}$. For Region I, we
require thermal misalignment to dominate over the “kick”, taking as our
criterion $\hat{\phi}_{\rm osc}(y_{c})>3\Delta\hat{\phi}_{\rm min}$. Using
Eqs. (III.5.1,38) we find that this is satisfied for

$\displaystyle\kappa>\gamma
y_{c}^{2}\left(\frac{3\pi^{12}}{\lambda^{4}}\right)^{1/5}\simeq 10^{3}\,.$
(39)

In practice we find that numerically we agree well with the estimates of
Region I for $\kappa>10^{3}$. We thus take $\kappa=10^{3}$ as the boundary for
Region I.

Having established that thermal misalignment dominates in Region I,
$\kappa\geq 10^{3}$, we are now ready to provide an analytic estimate of the
$\phi$ relic abundance. We use Eq. (30) to estimate $\rho_{\phi,0}$, the
$\phi$ energy density today. As input to this equation, we take $y_{f}=y_{\rm
osc}$ as the temperature at which oscillations begins, with $y_{\rm osc}$
given by Eq. (31). We also use $\rho(y_{f})=\rho(y_{\rm
osc})\simeq\tfrac{1}{2}m_{\phi}^{2}\phi(y_{\rm osc})^{2}$, with $\phi(y_{\rm
osc})$ given by the first term in Eq. (35). Other inputs needed are described
below Eq. (30). Putting everything together, we arrive at the following
estimate for the DM density parameter,

$\displaystyle\Omega_{\phi,0}$
$\displaystyle\simeq\frac{g_{*S}^{0}}{g_{*S}^{\rm
osc}}\,\frac{y_{0}^{3}\,\mu^{4}}{\rho_{c,0}}\,\frac{3^{3/2}}{8\pi^{4}}\,\frac{\beta^{2}}{\gamma^{1/2}\kappa^{3/2}}$
$\displaystyle\simeq
0.26\left(\frac{\beta}{0.05}\right)^{2}\left(\frac{1000}{\kappa}\right)^{3/2}\,.$
(40)

Parameters leading to the observed DM abundance according to the approximate
analytic estimate, Eq. (III.5.1), are displayed as a dashed line in the
$m_{\phi}-A$ plane in Fig. 1 for $m_{\phi}\gtrsim 10^{-3}$ eV, agreeing well
with the calculation using the exact numerical solution to the $\phi$ equation
of motion (28) in Region I. Eq. (III.5.1) (and the associated numerical
estimate) represents one of the most important results of this work, providing
a prediction for the observed DM relic density in terms of only the DM model
parameters ($m_{\phi}$ and $A$) which is insensitive to initial conditions
provided $|\phi_{i}|\ll|\phi_{T}|$, as discussed in Sec. III.3. In particular,
note that for parameters above the relic density line in Fig. 1, initial
conditions with $|\phi_{i}|\gg|\phi_{T}|$ will overclose the universe. Since
thermal misalignment is unavoidable given a standard radiation dominated
cosmology and a high enough reheating temperature $T_{RH}\gg v$, this gives a
cosmologically motivated target in the parameter space for experiments to
pursue. In fact, requiring DM not to overclose the universe currently provides
the best constraint on the model (with the caveats mentioned above) over the
vast majority of the parameter space.

####  III.5.2 Region II $(10\lesssim\kappa\lesssim 10^{3},~{}3\times
10^{-5}\lesssim m_{\phi}\lesssim 3\times 10^{-3}{\rm eV})$

We now move to Region II, corresponding to the intermediate mass range
$(10\lesssim\kappa\lesssim 10^{3},~{}3\times 10^{-5}\lesssim m_{\phi}\lesssim
3\times 10^{-3}{\rm eV})$. In this region, the final misalignment amplitude is
influenced in an important way by the change in the minimum of $\phi$
triggered by the EWPT at $y=y_{c}$.

Notably, for masses near the upper end of Region II, $m_{\phi}\sim{\cal
O}(10^{-3}{\rm eV})$, a competition between thermal misalignment and VEV
misalignment leads to a novel forced resonance phenomena which impacts the
relic abundance in a rather dramatic fashion. At the EWPT, the Higgs field
rapidly moves from the origin towards $h=v$, simultaneously inducing a shift
the $\phi$ potential minimum towards its zero-temperature value $\phi_{0}$.
The acts as a step-like forcing term in the $\phi$ equation of motion (29),
causing a suppression or enhancement in the oscillation amplitude depending on
the relative phase between the oscillations and Higgs source term, which in
turn depends on the the scalar mass. For instance, if the scalar is near its
oscillation maximum as this shift happens, the effective oscillation amplitude
is reduced, thus requiring a larger coupling $\beta$ to produce the observed
DM abundance. This behavior is shown in Fig. 5 for one example benchmark
point. A striking series of recurring peaks and valleys in the coupling
$\beta$ yielding the correct DM abundance is observed as the scalar mass is
varied, as shown the inset of Fig. 1.

For masses near the lower end of Region II, VEV misalignment dominates and the
predicted DM abundance starts to be sensitive to initial conditions, which is
also evident from Fig. 1. The impact of VEV misalignment and the role of the
initial conditions for the low mass scalars will be clarified in the next
through our analysis of Region III.

Figure 5: Left panel: The evolution of $\phi$ (blue) along with the
progression of the temperature-dependent minimum of $\phi$ (black dashed)
during the EWPT. The true minimum jumps quickly after the transition, reducing
the effective oscillation amplitude. Right panel: This amplitude reduction is
clearly seen in $\Delta\phi(y)\equiv\phi(y)-\phi_{\rm min}(y)$, which shows
the oscillations relative to the temperature dependent $\phi$ minimum, Eq.
(36). In both plots we have chosen the benchmark parameters
$\kappa=10^{2.985}$, $\beta=10^{-1}$.

####  III.5.3 Region III $(\kappa\lesssim 1,m_{\phi}\lesssim 3\times
10^{-6}\,{\rm eV})$

Finally, we come to Region III, which may roughly be defined by the condition
that oscillations begin below the EWPT, $y_{\rm osc}\lesssim
1\implies\kappa\lesssim 3\gamma$. We will simply take $\kappa\lesssim 1$
($m_{\phi}\lesssim 3\times 10^{-6}\,{\rm eV}$) for concreteness. In this low
mass regime, VEV misalignment caused by the EWPT dominates over thermal
misalignment. However, the final $\phi$ misalignment and number density
depends sensitively on the initial conditions (see Fig. 1). In the following,
we will examine two distinct choices for the initial condition,
$\hat{\phi}_{i}=\hat{\phi}_{0}$ and $\hat{\phi}_{i}=0$, as discussed in Sec.
III.3.

(i) $\hat{\phi}_{i}=\hat{\phi}_{0}$: Let us first consider the initial
condition $\hat{\phi}_{i}=\hat{\phi}_{0}\simeq-\beta/(2\lambda\kappa^{2})$.
Examining the $\hat{\phi}$ equation of motion (29) for
$\hat{\phi}\approx\hat{\phi}_{0}$, we observe that the mass term
$\kappa^{2}\hat{\phi}$ dominates over the thermal term for
$\frac{y^{2}}{2\pi^{2}}\beta\left(J^{\prime}_{B}[\eta_{h}]+3J^{\prime}_{B}[\eta_{\chi}]\right)$
for $y^{2}\lesssim\pi^{2}/(2\lambda)$. Thus, the scalar field first
experiences thermal misalignment
$\hat{\phi}_{T}\sim-\lambda\beta/(3\pi^{4}\gamma^{2})$, where we have used Eq.
(35) and $y^{2}\sim\pi^{2}/(2\lambda)$. However, we note that the displacement
of the scalar from $\hat{\phi}_{0}$ generated by thermal misalignment is
minuscule, $|\hat{\phi}_{T}|\ll\hat{\phi}_{0}$, since $\kappa\lesssim 1$ in
this region. For $y_{c}<y\lesssim\sqrt{\pi^{2}/(2\lambda)}\simeq 7$ the mass
term dominates and the equation of motion becomes

$\displaystyle\hat{\phi}^{\prime\prime}(y)+\frac{1}{\gamma^{2}y^{6}}\left(\kappa^{2}\hat{\phi}\right)=0\,,$
$\displaystyle\Longrightarrow~{}~{}~{}~{}\hat{\phi}^{\prime\prime}(y)-\frac{\beta}{2\gamma^{2}\lambda}\frac{1}{y^{6}}=0\,,$
(41)

where in the second line we have assumed $\hat{\phi}\simeq\hat{\phi}_{0}$. The
solution to this equation is given by

$\displaystyle\hat{\phi}(y)=\frac{1}{y^{4}}\frac{\beta}{40\gamma^{2}\lambda}+\hat{\phi}_{0}\,.$
(42)

The trajectory of $\hat{\phi}$ follows Eq. (42) until $y\sim 1$ at which point
the contribution from $\beta\hat{h}^{2}/2$ in Eq. (29) turns on and, to a good
approximation, cancels the $\kappa^{2}\hat{\phi}$ term. At this point, all
source terms in Eq. (42) are negligible and the evolution is dictated by
$\hat{\phi}^{\prime\prime}(y)=0$. This has the solution

$\displaystyle\hat{\phi}(y)=\hat{\phi}^{\prime}(1)(y-1)+\hat{\phi}(1)\,,$ (43)

where the initial conditions for the field and its velocity are obtained by
matching the solution in Eq. (42) at $y=1$. We find
$\hat{\phi}^{\prime}(1)=-\beta/(10\gamma^{2}\lambda)$ and
$\hat{\phi}(1)=\beta/(40\gamma^{2}\lambda)+\hat{\phi}_{0}$. Then for small
$y$, before beginning of oscillations, we obtain the asymptotic value of
$\hat{\phi}(y\ll 1)=5\beta/(40\gamma^{2}\lambda)+\hat{\phi}_{0}$, where the
first term is the amplitude of oscillations. The evolution of the scalar field
for $\hat{\phi}_{i}=\hat{\phi}_{0}$ is displayed for a representative
benchmark model in the left panel of Fig. 6.

Figure 6: Scalar field evolution in Region III, where VEV misalignment is
important, for two choices of initial conditions. Left panel: The initial
condition is chosen to be $\phi_{i}=\phi_{0}$. The evolution
$\Delta\phi(y)\equiv\phi(y)-\phi_{0}$ is shown for two benchmark models:
$\kappa=10^{-2}$, $\beta=10^{-4}$ (blue curves), and $\kappa=10^{-2}$,
$\beta=10^{-5}$ (orange curves). Right panel: The initial condition is chosen
to be $\phi_{i}=0$. The evolution $\phi(y)$ is shown for two benchmark models:
$\kappa=10^{-2}$, $\beta=10^{-6}$ (blue curves), and $\kappa=10^{-2}$,
$\beta=10^{-7}$ (orange curves). The dotted lines show the corresponding zero
temperature minima $\phi_{0}$.

With an estimate for the oscillation amplitude, $\hat{\phi}(y_{\rm osc})\simeq
5\beta/(40\gamma^{2}\lambda)$, we are now in a position to compute the $\phi$
relic abundance, following similar steps to those used earlier in Region I. In
particular, we employ Eq. (30) to estimate $\rho_{\phi,0}$ , taking
$y_{f}=y_{\rm osc}=\sqrt{\kappa/(3\gamma)}$, Eq. (31), as the temperature
corresponding to the onset of oscillations. Furthermore, we use
$\rho(y_{f})=\rho(y_{\rm osc})\simeq\tfrac{1}{2}m_{\phi}^{2}\phi(y_{\rm
osc})^{2}$. Collecting the other inputs described around Eq. (30), we obtain

$\displaystyle\Omega_{\phi,0}$
$\displaystyle\simeq\frac{g_{*S}^{0}}{g_{*S}^{\rm
osc}}\,\frac{y_{0}^{3}\,\mu^{4}}{\rho_{c,0}}\,\frac{3^{3/2}}{128}\,\frac{\beta^{2}\,\kappa^{1/2}}{\gamma^{5/2}\,\lambda^{2}}$
$\displaystyle\simeq
0.26\left(\frac{\beta}{10^{-4}}\right)^{2}\left(\frac{\kappa}{4\times
10^{-2}}\right)^{1/2}\,.$ (44)

It should be noted that since the amplitude of oscillations is comparatively
small for $\phi_{i}=\phi_{0}$, a much larger value of $\beta$ is needed to
obtain the correct relic density. This initial condition is representative of
an optimistic scenario for experimental detection of low mass scalars.

(ii) $\hat{\phi}_{i}=0$: Next, we consider the initial condition
$\hat{\phi}_{i}=0$. It is far more straightforward to estimate the scalar
relic density in this case. Thermal misalignment is again negligible, and in
the earliest stages of its evolution the scalar is held near at its initial
value by Hubble friction until the EWPT at $y_{c}$. At this point, the
$\hat{\phi}$ VEV rapidly evolves to its zero temperature value
$\hat{\phi}_{0}$, triggering VEV misalignment of the scalar, until eventually
oscillations begin (see Fig. 6 right panel). The oscillation amplitude is
therefore given by $\hat{\phi}(y_{\rm
osc})\simeq\hat{\phi}_{0}\simeq-\beta/(2\lambda\kappa^{2})$. Thus the $\phi$
density parameter today is given by

$\displaystyle\Omega_{\phi,0}$
$\displaystyle\simeq\frac{g_{*S}^{0}}{g_{*S}^{\rm
osc}}\,\frac{y_{0}^{3}\,\mu^{4}}{\rho_{c,0}}\,\frac{3^{3/2}}{8}\,\frac{\gamma^{3/2}\beta^{2}}{\lambda^{2}\,\kappa^{7/2}}$
$\displaystyle\simeq 0.26\left(\frac{\beta}{3\times
10^{-10}}\right)^{2}\left(\frac{10^{-2}}{\kappa}\right)^{7/2}\,.$ (45)

For this initial condition, we note the drastically lower values of $\beta$
required to get the correct relic density. Thus the initial condition
$\phi_{i}=0$ is representative of low mass scalar scenario that is challenging
to probe experimentally.

###  III.6 Summary of Results

A summary of our results is presented in Fig. 1, which shows contours leading
to the observed relic abundance in the $m_{\phi}-A$ plane for the two choices
of initial conditions, $\phi_{i}=\phi_{0}$ and $\phi_{i}=0$. Both the
numerical results from solving Eq. (29) and using Eq. (30) (solid lines), as
well as the analytic estimates, Eqs. (III.5.1,III.5.3,III.5.3) (dashed lines)
are presented and agree well.

##  IV Constraints

We now discuss the existing experimental constraints and future prospects for
the ultra-light Higgs portal DM model. These bounds are shown in the
$m_{\phi}-A$ plane and can be compared with our numerical and analytic
predictions for the relic abundance lines ($\Omega_{\phi,0}=\Omega_{\rm
DM}=0.26$) in Fig. 1. We now provide a brief overview for each of these
bounds.

###  IV.1 Equivalence Principle and Inverse Square Law Tests

The light scalar mixes with the Higgs and thus couples to ordinary matter,
leading to an additional Yukawa-like potential between nucleons. This
additional force can be constrained by experiments searching for violations of
the gravitational inverse square law (ISL) or tests of the equivalence
principle (EP).

For the ISL bounds, we follow the approach of Ref. Piazza and Pospelov (2010).
The potential between two test bodies (labeled 1 and 2) is given by

$\displaystyle
V(r)=-\frac{Gm_{1}m_{2}}{r}\left(1+\alpha_{1}\alpha_{2}\,e^{-m_{\phi}r}\right)\,.$
(46)

We will assume that the scalar has the same coupling to protons and neutrons,
with the scalar-nucleon coupling given by

$\displaystyle g_{\phi NN}=g_{hNN}\,\frac{A\,v}{m^{2}_{h}}\simeq
10^{-3}\,\frac{A\,v}{m^{2}_{h}}\,,$ (47)

where $g_{hNN}\simeq 10^{-3}$ is known with order one uncertainty Graham  _et
al._ (2016). The scalar coupling to the test body $i$ is given by

$\displaystyle\frac{\alpha_{i}}{\sqrt{2}M_{\rm pl}}=\frac{d\ln
m_{N}(\phi)}{d\phi}=\frac{g_{\phi NN}}{m_{N}}\,.$ (48)

We can then translate constraints reported on ISL tests. The constraints are
typically reported on the product $\tilde{\alpha}=\alpha_{1}\alpha_{2}$, which
in terms of our model parameters is given by

$\displaystyle\tilde{\alpha}=\alpha_{1}\alpha_{2}=10^{-6}\,\frac{2M^{2}_{\rm
pl}}{m^{2}_{N}}\left(\frac{A\,v}{m^{2}_{h}}\right)^{2}=\left(\frac{A}{1.7\times
10^{-5}\,{\rm eV}}\right)^{2}.$ (49)

The constraints from ISL tests from Refs. Hoskins  _et al._ (1985); Yang  _et
al._ (2012) are shown in Fig. 1. As can be seen in the inset, the ISL tests
are already beginning to probe interesting regions of the parameter space due
to the resonance effects discussed in Sec. III.5.2.

Constraints from violation of EP are more involved. We follow Ref. Graham  _et
al._ (2016) in the following. EP constraints involve measurements of
differential accelerations in test bodies which have different charge-to-mass
ratios, i.e., the ratio of the scalar coupling strength to the mass of a given
atom ($\Delta_{\phi H^{2}}$). The EP violation constraints on a mediator
coupled to a $B-L$ charge obtained by measuring the differential acceleration
between Beryllium (Be) and Aluminium (Al) were derived in Ref. Wagner  _et
al._ (2012). Translating these bounds to our model, we can write

$\displaystyle A=4.4\times 10^{13}{\rm
eV}\,g_{B-L}\,\sqrt{\frac{\Delta_{B-L}}{\Delta_{\Phi h^{2}}}}\,,$ (50)

where $g_{B-L}$ is the the coupling strength of the scalar to $B-L$ charge and
$\Delta_{B-L}$ is the differential charge-to-mass ratio between Be and Al. We
use $\Delta_{B-L}=0.037$ and $\Delta_{\phi H^{2}}=4\times 10^{-4}$ Graham  _et
al._ (2016). The EP constraints on the Higgs portal model are shown in Fig. 1.

###  IV.2 Stellar Cooling Bounds

For the Higgs portal, the stellar cooling bounds primarily arise from the
effective coupling of $\phi$ to electrons Hardy and Lasenby (2017). The
dominant constraints on the model come from red giant (RG) and horizontal
branch (HB) stars. In red giant cores, the dominant energy loss mechanism is
neutrino production. Any additional stronger energy loss mechanism will delay
the onset of helium fusion in cores of RG stars. Since this is not observed,
it implies a constraint on the Higgs-scalar mixing angle $\sin\theta\simeq
Av/m^{2}_{h}\leq 3\times 10^{-10}$ for $m_{\phi}<2~{}{\rm keV}$ Hardy and
Lasenby (2017). For HB stars, any additional strong energy loss mechanism
causes the cores to contract and heat up, increasing the rate of helium fusion
and shortening the lifetime of the stars, which is not observed. The combined
constraints from RG and HB stars are shown in Fig. 1.

###  IV.3 Resonant Absorption in Molecules

Resonant absorption of $\phi$ by polyatomic molecules has been proposed as a
DM detection concept in Ref. Arvanitaki  _et al._ (2018). The basic idea is to
search for DM absorption via the subsequent transition of these molecules from
the excited to ground state through the emission of photons. Concerning the
Higgs portal model, this technique is most sensitive to the effective coupling
of $\phi$ to electrons, which in terms of the coupling $A$ is given by

$\displaystyle g_{\phi
ee}=\frac{A\,v}{m^{2}_{h}}y_{e}=\frac{A\,m_{e}}{m^{2}_{h}}.$ (51)

The approach has the potential to probe scalar DM coupled to electrons in the
mass range $0.2\,{\rm eV}\lesssim m_{\phi}\lesssim 20\,{\rm eV}$. In the
parameterization employed in Arvanitaki  _et al._ (2018), the constraint is
given by $d_{m_{e}}\lesssim 10^{5}$ where $d_{m_{e}}=\sqrt{2}g_{\phi ee}M_{\rm
pl}/m_{e}\simeq 2.2\times 10^{5}\times(A/{\rm eV})$.

We show projections corresponding to a phase II ‘bulk’ (BII) and ‘stack’ (SII)
experimental configurations in Fig. 1; for more details we refer the reader to
Ref. Arvanitaki  _et al._ (2018).

###  IV.4 Extragalactic Background Light and Reionization Constraints

For larger masses and couplings, the scalar $\phi$ can decay into photons on
cosmological time scales. The decay rate of the scalar to photons is given by
Fradette  _et al._ (2019)

$\displaystyle\Gamma_{\phi\gamma\gamma}=\frac{\theta^{2}\alpha^{2}}{256\pi^{2}}\frac{m^{3}_{\phi}}{v^{2}}|C|^{2}$
(52)

where $\theta\simeq Av/m^{2}_{h}$ and the factor $C\simeq 50/27$ accounts for
the charged particles running in the loop. Decay of scalars would contribute
to the extragalactic background light (EBL). EBL constraints have been
calculated in Cadamuro and Redondo (2012); Arias  _et al._ (2012) and recast
for a Higgs portal in Flacke  _et al._ (2017). Note that these constrains are
calculated assuming a given present day DM density
$Y_{\phi}=n_{\phi}/s_{\gamma}$, where $s_{\gamma,0}$ is the entropy density of
photons today. The actual parameter constrained is
$m_{\phi}n_{\phi}\Gamma_{\phi}$ as a function of $m_{\phi}$. We rescale our
constraints assuming $\Omega_{\phi}=\Omega_{DM}\simeq 0.26$ i.e.
$Y_{\phi}=0.26\frac{\rho_{c}}{s_{\gamma,0}}m_{\phi}$, where $\rho_{c}$ is the
critical energy density of universe today.

###  IV.5 X-ray Constraints

For $m_{\phi}\gtrsim$ few keV, the scalar decay to photons on cosmological
time scales is constrained by the HEAO Gruber  _et al._ (1999) and INTEGRAL
Bouchet  _et al._ (2008) satellites. The data from these satellites was
translated into lifetime of scalar decaying into photons in Essig  _et al._
(2013) and we use translate these constraints on the lifetime of our scalar
$\tau_{\phi}=\Gamma^{-1}_{\phi\gamma\gamma}$ given by Eq. (52).

###  IV.6 Other Probes

There are other constraints that we have not shown because they lie outside
the shown mass range, or lie in the forbidden region
$A>\sqrt{2\lambda}m_{\phi}$. These include atomic and nuclear clocks, atomic
interferometers, mechanical resonators, superradiance etc. A comprehensive
summary of all constraints can be found in Antypas  _et al._ (2022).

##  V Conclusions

In this work we have revisited the cosmology of a light scalar feebly coupled
to the SM through the super-renormalizable Higgs portal. This is among the
most minimal UV complete extensions of the SM and is described by only two
additional parameters, the scalar mass $m_{\phi}$ and its dimensionful
coupling to the Higgs $A$. The scalar field in this model can naturally be
light due to its feeble super-renormalizable coupling. Despite being minimal,
the model presents a rich cosmology, with the scalar field experiencing a non-
trivial dynamical evolution during the radiation era before and/or during
scalar oscillations. This generates two additional dynamical sources of
misalignment, which we have referred to as thermal misalignment and VEV
misalignment, beyond that from its initial field value at the end of
inflation. We have studied the evolution of the scalar field both numerically
and analytically over a broad range of scalar masses and couplings,
delineating the parameters in the $m_{\phi}-A$ plane that are consistent with
the observed DM relic density. We have also investigated two qualitatively
distinct choices of initial scalar field values in order to discern the role
played by initial conditions. Our main results are presented in Fig. 1.

For large masses, the scalar field evolution is mainly governed by thermal
misalignment. The finite-temperature contribution to the effective potential
drives the scalar towards large field values at high temperatures, dynamically
generating misalignment. Provided the initial scalar field value is smaller in
magnitude than its displacement caused by thermal misalignment, the scalar
relic density in this higher mass range is insensitive to the initial
conditions. The scalar abundance today is then tightly controlled by the model
parameters $m_{\phi}$ and $A$. Therefore, parameters in the mass-coupling
plane yielding the observed relic abundance then provide a robust target that
can be compared against experimental searches for ultra-light scalar DM.
Furthermore, the relic density line in this mass range can be viewed as an
upper bound on the scalar-Higgs coupling assuming a standard cosmology with a
high enough reheat temperature, as any larger coupling would overclose the
universe. These features are quite similar in character to those of the
popular WIMP DM scenario, in which the relic abundance is insensitive to
initial conditions and UV dynamics, instead being set by the WIMP mass and its
couplings to SM particles through the mechanism of thermal freeze-out. The
insensitivity to initial conditions in the high scalar mass range is evident
from Fig. 1.

For lower scalar masses, the impact of the EWPT on the scalar field evolution
is important. The phase transition leads to a rapid shift in the scalar field
potential minimum towards its zero-temperature VEV, generating VEV
misalignment. However, in this regime the oscillation amplitude and
corresponding relic density is also sensitive to the initial conditions and
depends on whether the scalar field is initially in close proximity to its
zero temperature minimum or not. In the former case, which can be naturally
realized for a long period of inflation with a low inflationary Hubble
parameter (smaller than the electroweak scale $v$), relatively large couplings
are needed to significantly displace the field from its minimum and produce
the observed DM abundance, thus offering more promising detection prospects.
Instead, if the field is initially distant from its zero temperature VEV, the
required coupling can be extremely small while still leading to a consistent
DM cosmology.

In the intermediate mass range we have discovered a novel forced resonance
phenomena resulting from a competition between thermal misalignment and VEV
misalignment. This causes enhancements or suppressions in the oscillation
amplitude depending on the scalar mass, thus requiring smaller or larger
couplings, respectively, to achieve the correct DM abundance. This is
reflected by a series of peaks and valleys in the DM relic density contours as
the scalar mass varies; see Fig. 1.

We have also compared our relic density predictions with constraints and
projections from a variety of terrestrial and astrophysical probes, including
the equivalence principle and inverse square law tests, stellar cooling,
resonant molecular absorption, and observations of extra-galactic background
light and diffuse X-ray backgrounds. While some parts of the parameter space
are starting to be explored, new experimental ideas are needed to probe the
cosmologically motivated regions of parameter space. We hope our results will
provide the impetus for new creative approaches to ultralight scalar DM
detection.

###### Acknowledgements.

We thank Daniel Boyanovsky, Raymond Co, and Erwin Tanin for helpful
discussions. The work of B.B., A.G., and M.R. is supported by the U.S.
Department of Energy under grant No. DE–SC0007914. M.R. thanks ICTP, Trieste
for hospitality while this work was being completed.

##  Appendix A Thermal masses

In this appendix we discuss the thermal masses $m_{i}^{2}(\phi,h,T)$ entering
in the finite temperature potential, Eq. (III.1). As discussed in Section
III.1, we employ the daisy resummation scheme of Ref. Parwani (1992), writing

$\displaystyle m_{i}^{2}(\phi,h,T)=m^{2}_{0,i}(\phi,h)+\Pi_{i}(T)\,.$ (53)

We keep the leading contributions in the high temperature expansion for the
self-energies, using the results presented in Ref. Espinosa  _et al._ (1993).

For the Higgs and Nambu-Goldstone scalars, these are

$\displaystyle m_{0,h}^{2}(\phi,h)$
$\displaystyle=-\mu^{2}+3\,\lambda\,h^{2}+A\,\phi\,,$ (54) $\displaystyle
m_{0,\chi}^{2}(\phi,h)$ $\displaystyle=-\mu^{2}+\lambda\,h^{2}+A\,\phi\,,$
(55) $\displaystyle\Pi_{h}(T)$
$\displaystyle=\Pi_{\chi}(T)=T^{2}\left(\frac{3}{16}g^{2}+\frac{1}{16}g^{\prime
2}+\frac{1}{4}y_{t}^{2}+\frac{1}{2}\lambda^{2}\right)\,,$ (56)

where $\mu^{2}$, $\lambda$, $A$ are parameters entering in the scalar
potential, Eq. (1), and $g$, $g^{\prime}$, and $y_{t}$ are the $SU(2)_{L}$
gauge coupling, $U(1)_{Y}$ gauge coupling, and top Yukawa coupling,
respectively.

In gauge boson sector, the squared masses for the transverse components are
simply

$\displaystyle m^{2}_{W_{T}}(h)=m_{0,W}^{2}(h)$
$\displaystyle=\frac{1}{4}g^{2}h^{2}\,,$ (57) $\displaystyle
m^{2}_{Z_{T}}(h)=m_{0,Z}^{2}(h)$ $\displaystyle=\frac{1}{4}(g^{2}+g^{\prime
2})h^{2}\,.$ (58)

The longitudinal gauge boson components receive a self-energy correction,
expressed in the gauge basis as $\Pi^{L}_{GB}(T)=\frac{11}{6}T^{2}{\rm
diag}(g^{2},g^{2},g^{2},g^{\prime 2})$. The thermal masses are thus given by
the eigenvalues of the squared mass matrix

$\displaystyle\left(\begin{array}[]{cccc}\tfrac{1}{4}g^{2}h^{2}&0&0&0\\\
0&\tfrac{1}{4}g^{2}h^{2}&0&0\\\
0&0&\tfrac{1}{4}g^{2}h^{2}&-\tfrac{1}{4}gg^{\prime}h^{2}\\\
0&0&-\tfrac{1}{4}gg^{\prime}h^{2}&\tfrac{1}{4}g^{\prime
2}h^{2}\end{array}\right)+\left(\begin{array}[]{cccc}\tfrac{11}{6}g^{2}T^{2}&0&0&0\\\
0&\tfrac{11}{6}g^{2}T^{2}&0&0\\\ 0&0&\tfrac{11}{6}g^{2}T^{2}&0\\\
0&0&0&\tfrac{11}{6}g^{\prime 2}T^{2}\end{array}\right).$ (67)

Finally, for the top quark we have

$\displaystyle m^{2}_{t}(h)=m_{0,t}^{2}(h)$
$\displaystyle=\frac{1}{2}y_{t}^{2}h^{2}\,.$ (68)

For completeness, we also collect here the expressions for the dimensionless
arguments $\eta_{i}(\hat{\phi},\hat{h},y)\equiv m^{2}_{i}(\phi,h,T)/T^{2}$ of
the $J_{B,F}$ functions. Using the results above and the dimensionless
parameters defined in Eq. (20), we obtain

$\displaystyle\eta_{h}$
$\displaystyle=\frac{1}{y^{2}}\left[3\lambda\hat{h}^{2}-(1-\beta\hat{\phi})\right]+\left(\frac{3}{16}g^{2}+\frac{1}{16}g^{\prime
2}+\frac{1}{4}y^{2}_{t}+\frac{1}{2}\lambda\right),$ (69)
$\displaystyle\eta_{\chi}$
$\displaystyle=\frac{1}{y^{2}}\left[\lambda\hat{h}^{2}-(1-\beta\hat{\phi})\right]+\left(\frac{3}{16}g^{2}+\frac{1}{16}g^{\prime
2}+\frac{1}{4}y^{2}_{t}+\frac{1}{2}\lambda\right),$ (70)
$\displaystyle\eta_{W_{T}}$ $\displaystyle=\frac{g^{2}\hat{h}^{2}}{4y^{2}},$
(71) $\displaystyle\eta_{Z_{T}}$ $\displaystyle=\frac{(g^{2}+g^{\prime
2})\hat{h}^{2}}{4y^{2}},$ (72) $\displaystyle\eta_{W_{L}}$
$\displaystyle=\frac{g^{2}\hat{h}^{2}}{4y^{2}}+\frac{11}{6}g^{2},$ (73)
$\displaystyle\eta_{Z_{L}}$
$\displaystyle=\frac{1}{2\,y^{2}}\left[\left(g^{2}+g^{\prime
2}\right)\left(\tfrac{1}{4}\hat{h}^{2}+\tfrac{11}{6}y^{2}\right)+\sqrt{\left(g^{2}-g^{\prime
2}\right)^{2}\left(\tfrac{1}{4}\hat{h}^{2}+\tfrac{11}{6}y^{2}\right)^{2}+\tfrac{1}{4}g^{2}g^{\prime
2}{\hat{h}}^{4}~{}}\,\right],$ (74) $\displaystyle\eta_{A_{L}}$
$\displaystyle=\frac{1}{2\,y^{2}}\left[\left(g^{2}+g^{\prime
2}\right)\left(\tfrac{1}{4}\hat{h}^{2}+\tfrac{11}{6}y^{2}\right)-\sqrt{\left(g^{2}-g^{\prime
2}\right)^{2}\left(\tfrac{1}{4}\hat{h}^{2}+\tfrac{11}{6}y^{2}\right)^{2}+\tfrac{1}{4}g^{2}g^{\prime
2}{\hat{h}}^{4}~{}}\,\right],$ (75) $\displaystyle\eta_{t}$
$\displaystyle=\frac{y_{t}^{2}\hat{h}^{2}}{2y^{2}}.$ (76)

## References

  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Antypas  _et al._ (2022) D. Antypas _et al._ , (2022), arXiv:2203.14915 [hep-ex] .
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Piazza and Pospelov (2010) F. Piazza and M. Pospelov, Phys. Rev. D 82, 043533 (2010), arXiv:1003.2313 [hep-ph] .
  * Graham  _et al._ (2016) P. W. Graham, D. E. Kaplan, J. Mardon, S. Rajendran, and W. A. Terrano, Phys. Rev. D 93, 075029 (2016), arXiv:1512.06165 [hep-ph] .
  * Arvanitaki  _et al._ (2018) A. Arvanitaki, S. Dimopoulos, and K. Van Tilburg, Phys. Rev. X 8, 041001 (2018), arXiv:1709.05354 [hep-ph] .
  * Batell and Ghalsasi (2021) B. Batell and A. Ghalsasi, (2021), arXiv:2109.04476 [hep-ph] .
  * Chun (2022) E. J. Chun, Phys. Lett. B 825, 136880 (2022), arXiv:2109.07423 [hep-ph] .
  * Brzeminski  _et al._ (2021) D. Brzeminski, Z. Chacko, A. Dev, and A. Hook, Phys. Rev. D 104, 075019 (2021), arXiv:2012.02787 [hep-ph] .
  * Arkani-Hamed  _et al._ (2021) N. Arkani-Hamed, R. T. D’Agnolo, and H. D. Kim, Phys. Rev. D 104, 095014 (2021), arXiv:2012.04652 [hep-ph] .
  * Coleman and Weinberg (1973) S. R. Coleman and E. J. Weinberg, Phys. Rev. D 7, 1888 (1973).
  * Dolan and Jackiw (1974) L. Dolan and R. Jackiw, Phys. Rev. D 9, 3320 (1974).
  * Weinberg (1974) S. Weinberg, Phys. Rev. D 9, 3357 (1974).
  * Quiros (1999) M. Quiros, in  _ICTP Summer School in High-Energy Physics and Cosmology_ (1999) arXiv:hep-ph/9901312 .
  * Parwani (1992) R. R. Parwani, Phys. Rev. D 45, 4695 (1992), [Erratum: Phys.Rev.D 48, 5965 (1993)], arXiv:hep-ph/9204216 .
  * Kajantie  _et al._ (1996) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Phys. Rev. Lett. 77, 2887 (1996), arXiv:hep-ph/9605288 .
  * Graham and Scherlis (2018) P. W. Graham and A. Scherlis, Phys. Rev. D 98, 035017 (2018), arXiv:1805.07362 [hep-ph] .
  * Takahashi  _et al._ (2018) F. Takahashi, W. Yin, and A. H. Guth, Phys. Rev. D 98, 015042 (2018), arXiv:1805.08763 [hep-ph] .
  * Tenkanen (2019) T. Tenkanen, Phys. Rev. Lett. 123, 061302 (2019), arXiv:1905.01214 [astro-ph.CO] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Yokoyama and Linde (1999) J. Yokoyama and A. D. Linde, Phys. Rev. D 60, 083509 (1999), arXiv:hep-ph/9809409 .
  * Bödeker and Nienaber (2022) D. Bödeker and J. Nienaber, Phys. Rev. D 106, 056016 (2022), arXiv:2205.14166 [hep-ph] .
  * Hoskins  _et al._ (1985) J. K. Hoskins, R. D. Newman, R. Spero, and J. Schultz, Phys. Rev. D 32, 3084 (1985).
  * Yang  _et al._ (2012) S.-Q. Yang, B.-F. Zhan, Q.-L. Wang, C.-G. Shao, L.-C. Tu, W.-H. Tan, and J. Luo, Phys. Rev. Lett. 108, 081101 (2012).
  * Wagner  _et al._ (2012) T. A. Wagner, S. Schlamminger, J. H. Gundlach, and E. G. Adelberger, Classical and Quantum Gravity 29, 184002 (2012).
  * Hardy and Lasenby (2017) E. Hardy and R. Lasenby, JHEP 02, 033 (2017), arXiv:1611.05852 [hep-ph] .
  * Fradette  _et al._ (2019) A. Fradette, M. Pospelov, J. Pradler, and A. Ritz, Phys. Rev. D 99, 075004 (2019), arXiv:1812.07585 [hep-ph] .
  * Cadamuro and Redondo (2012) D. Cadamuro and J. Redondo, JCAP 02, 032 (2012), arXiv:1110.2895 [hep-ph] .
  * Arias  _et al._ (2012) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, JCAP 06, 013 (2012), arXiv:1201.5902 [hep-ph] .
  * Flacke  _et al._ (2017) T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta, and G. Perez, JHEP 06, 050 (2017), arXiv:1610.02025 [hep-ph] .
  * Gruber  _et al._ (1999) D. E. Gruber, J. L. Matteson, L. E. Peterson, and G. V. Jung, Astrophys. J. 520, 124 (1999), arXiv:astro-ph/9903492 .
  * Bouchet  _et al._ (2008) L. Bouchet, E. Jourdain, J. P. Roques, A. Strong, R. Diehl, F. Lebrun, and R. Terrier, Astrophys. J. 679, 1315 (2008), arXiv:0801.2086 [astro-ph] .
  * Essig  _et al._ (2013) R. Essig, E. Kuflik, S. D. McDermott, T. Volansky, and K. M. Zurek, JHEP 11, 193 (2013), arXiv:1309.4091 [hep-ph] .
  * Espinosa  _et al._ (1993) J. R. Espinosa, M. Quiros, and F. Zwirner, Phys. Lett. B 314, 206 (1993), arXiv:hep-ph/9212248 .


$\overline{b}_{n}$, which thus should be $\approx 13.081$. However, forbidding
only simple windmills does not seem to give a close upper bound on the
exponential growth rate of guillotine rectangulations. Indeed, from the table
of the initial counting coefficients $g_{n}$, and applying acceleration of
convergence (see, e.g., [50, Sec.6]) to the ratio $g_{n+1}/g_{n}$, the
exponential growth rate of $g_{n}$ seems to be $\approx 10.24$.

Finally, we discuss lower bounds. By Proposition 4, the number of weak
guillotine rectangulations of size $n$ is the $(n-1)$th Schröder number.
Therefore, the exponential growth rate of Schröder numbers,
$3+2\sqrt{2}\approx 5.828$, is a “trivial” lower bound on the exponential
growth rate of strong guillotine rectangulations. In order to give a better
bound, we consider weak guillotine rectangulations where every 2-sided segment
(that is, a segment with at least one neighbor on each side) has weight $2$.
This will give a lower bound, since the neighbors of every 2-sided segment can
be shuffled in at least two ways. We adapt the decomposition from our proof of
Proposition 4 as follows. Let $G=G(x,y)$ be the generating function for weak
guillotine rectangulations, where the variable $x$ is for the size, and $y$
for the number of 2-sided segments. Further, let $G_{0}=G_{0}(x,y)$ be the
generating function for weak guillotine rectangulations that have no segment
with an endpoint on the left side of ${R}$, and let $G_{1}=G_{1}(x,y)$ be the
generating function for weak guillotine rectangulations that have at least one
segment with an endpoint on the left side of ${R}$. Finally, let $H=H(x,y)$
and $V=V(x,y)$ be the generating functions for horizontal and, respectively,
vertical weak guillotine rectangulations. Then we have $H=V$, $G=x+H+V$,
$G_{0}=xG+x$, and $G_{1}=(1-x)G-x$, and the decomposition of a vertical
guillotine rectangulation by its leftmost cut leads to the following weighted
version of equation (1):

$V=xG+H\Big{(}G_{0}+yG_{1}\Big{)}.$

The solution of this system yields

$G(x,2)=\frac{1+x-x^{2}-\sqrt{1-6x-5x^{2}+2x^{3}+x^{4}}}{2(2-x)},$

and its dominant singularity gives us the following lower bound.

######  Proposition 36.

The exponential growth rate of the number of strong guillotine rectangulations
is bounded from below by
$\frac{1}{2}\Big{(}1+\sqrt{13-8\sqrt{2}}\Big{)}\Big{(}3+2\sqrt{2}\Big{)}\approx
6.699$.

This strategy can be pushed further: for a fixed threshold value $t\geq 1$,
every segment with $i$ neighbors on one side and $j$ neighbors on the other
side is weighted by $\binom{i^{\prime}+j^{\prime}}{i^{\prime}}$, where
$i^{\prime}=\mathrm{min}(i,t)$ and $j^{\prime}=\mathrm{min}(j,t)$. The
exponential growth rate for any fixed $t$ should be computable (by the
approach of [43, VII.6.3]), and grow with $t$, giving better and better lower
bounds. The complexity of the decomposition and of computations, however, will
rapidly explode as $t$ grows.

## Acknowledgments

The work on this paper began during the Order & Geometry Workshop, Ciążeń,
Poland, 13–18 September 2022. The research of Andrei Asinowski is partially
supported by FWF — The Austrian Science Fund, grant P32731. The research of
Stefan Felsner is partially supported by grant DFG FE 340/13-1. The research
of Éric Fusy is partially supported by the projects ANR-20-CE48-0018 (3DMaps)
and ANR-19-CE48-0011 (COMBINÉ).

## References

  * [1] Michio Abe.  Covering the square by squares without overlapping.  J. Japan Math. Phys., 4:359–366, 1930.
  * [2] Michio Abe.  On the problem to cover simply and without gap the inside of a square with a finite number of squares which are all different from one another.  Proc. Phys.-Math. Soc. Japan, 14:385–387, 1932.
  * [3] Eyal Ackerman, Gill Barequet, and Ron Y. Pinter.  A bijection between permutations and floorplans, and its applications.  Discret. Appl. Math., 154(12):1674–1684, 2006.
  * [4] Eyal Ackerman, Gill Barequet, Ron Y. Pinter, and Dan Romik.  The number of guillotine partitions in _d_ dimensions.  Inf. Process. Lett., 98(4):162–167, 2006.
  * [5] Andrei Asinowski and Cyril Banderier.  From geometry to generating functions: rectangulations and permutations.  arXiv:2401.05558.
  * [6] Andrei Asinowski, Gill Barequet, Mireille Bousquet-Mélou, Toufik Mansour, and Ron Y. Pinter.  Orders induced by segments in floorplans and (2–14–3, 3–41–2)-avoiding permutations.  Electron. J. Comb., 20(2):35, 2013.
  * [7] Andrei Asinowski and Toufik Mansour.  Separable $d$-permutations and guillotine partitions.  Annals of Combinatorics, 14:17–43, 2010.
  * [8] Kirby A. Baker, Peter C. Fishburn, and Fred S. Roberts.  Partial orders of dimension $2$.  Networks, 2(1):11–28, 1972.
  * [9] Glen Baxter.  On fixed points of the composite of commuting functions.  Proc. AMS, 15(6):851–855, 1964.
  * [10] Glen Baxter and James T. Joichi.  On permutations induced by commuting functions, and an embedding question.  Math. Scandin., 13:140–150, 1963.
  * [11] Olivier Beaumont, Vincent Boudet, Fabrice Rastello, and Yves Robert.  Partitioning a square into rectangles: NP-completeness and approximation algorithms.  Algorithmica, 34(3):217–239, 2002.
  * [12] Benjamin B. Bederson, Ben Shneiderman, and Martin Wattenberg.  Ordered and quantum treemaps: Making effective use of 2D space to display hierarchies.  ACM Trans. Graph., 21(4):833–854, 2002.
  * [13] David Bevan.  Permutation patterns: basic definitions and notation, 2015.  arXiv:1506.06673.
  * [14] Anders Björner and Michelle L. Wachs.  Permutation statistics and linear extensions of posets.  J. Comb. Theory, Ser. A, 58(1):85–114, 1991.
  * [15] Nicolas Bonichon, Mireille Bousquet-Mélou, and Éric Fusy.  Baxter permutations and plane bipolar orientations.  Séminaire Lotharingien de Combinatoire, 61:B61Ah, 2010.
  * [16] Prosenjit Bose, Jonathan F. Buss, and Anna Lubiw.  Pattern matching for permutations.  Inf. Process. Lett., 65(5):277–283, 1998.
  * [17] Alin Bostan, Kilian Raschel, and Bruno Salvy.  Non-D-finite excursions in the quarter plane.  J. Comb. Theory, Ser. A, 121:45–63, 2014.
  * [18] Mireille Bousquet-Mélou.  Four classes of pattern-avoiding permutations under one roof: Generating trees with two labels.  The Electron. J. Comb., 9(2):R19, 2003.
  * [19] Mathilde Bouvel, Veronica Guerrini, Andrew Rechnitzer, and Simone Rinaldi.  Semi-Baxter and Strong-Baxter: Two Relatives of the Baxter Sequence.  SIAM J. Discret. Math., 32(4):2795–2819, 2018.
  * [20] Mathilde Bouvel, Veronica Guerrini, and Simone Rinaldi.  Avoiding Baxter-like patterns.  The 17th Conference on Permutation Patterns, June 2019, Zürich, Switzerland.
  * [21] Petter Brändén and Anders Claesson.  Mesh patterns and the expansion of permutation statistics as sums of permutation patterns.  Electron. J. Comb., 18(2), 2011.
  * [22] R. Leonard Brooks, Cedric A. B. Smith, Arthur H. Stone, and William T. Tutte.  The dissection of rectangles into squares.  Duke Math. J., 7:312–340, 1940.
  * [23] Kevin Buchin, Bettina Speckmann, and Sander Verdonschot.  Evolution strategies for optimizing rectangular cartograms.  In Ningchuan Xiao, Mei-Po Kwan, Michael F. Goodchild, and Shashi Shekhar, editors, Proc. GIScience 2012, volume 7478 of LNCS, pages 29–42. Springer, 2012.
  * [24] Sophie Burrill, Julien Courtiel, Éric Fusy, Stephen Melczer, and Marni Mishna.  Tableau sequences, open diagrams, and Baxter families.  Europ. J. Comb., 58:144–165, 2016.
  * [25] Felipe C. Calheiros, Abilio Lucena, and Cid C. de Souza.  Optimal rectangular partitions.  Networks, 41(1):51–67, 2003.
  * [26] Jean Cardinal, Vera Sacristán, and Rodrigo I. Silveira.  A note on flips in diagonal rectangulations.  Discret. Math. Theor. Comput. Sci., 20(2), 2018.
  * [27] Cesar Ceballos, Francisco Santos, and Günter M. Ziegler.  Many non-equivalent realizations of the associahedron.  Comb., 35(5):513–551, 2015.
  * [28] Tung-Chieh Chen and Yao-Wen Chang.  Chapter 10 — Floorplanning.  In Laung-Terng Wang, Yao-Wen Chang, and Kwang-Ting (Tim) Cheng, editors, Electronic Design Automation, pages 575–634. Morgan Kaufmann, 2009\.
  * [29] Fan-Rong King Chung, Ronald Lewis Graham, Verner Emil Hoggatt, and Mark Kleiman.  The number of Baxter permutations.  J. Comb. Theory, Ser. A, 24:382–394, 1978.
  * [30] Jim Conant and Tim Michaels.  On the number of tilings of a square by rectangles.  Annals of Comb., 18:21–34, 2014.
  * [31] Hubert de Fraysseix, Patrice Ossona de Mendez, and Pierre Rosenstiehl.  Bipolar orientations revisited.  Discret. Appl. Math., 56(2-3):157–179, 1995.
  * [32] Denis Denisov and Vitali Wachtel.  Random walks in cones.  Annals of Prob., 43(3):992–1044, 2015.
  * [33] Blanche Descartes.  Division of a square into rectangles.  Eureka, (34):31–35, 1971.
  * [34] Enrica Duchi and Corentin Henriet.  A bijection between rooted planar maps and generalized fighting fish.  arXiv:2210.16635, 2022.
  * [35] Serge Dulucq, S. Gire, and Olivier Guibert.  A combinatorial proof of J. West’s conjecture.  Discret. Math., 187(1-3):71–96, 1998.
  * [36] Serge Dulucq, S. Gire, and J. West.  Permutations with forbidden subsequences and nonseparable planar maps.  Discret. Math., 153(1-3):85–103, 1996.
  * [37] Andrzej Ehrenfeucht and Grzegorz Rozenberg.  T-structures, T-functions, and texts.  Theor. Comput. Sci., 116(2):227–290, 1993.
  * [38] David Eppstein, Elena Mumford, Bettina Speckmann, and Kevin Verbeek.  Area-universal and constrained rectangular layouts.  SIAM J. Comput., 41(3):537–564, 2012.
  * [39] Wenjie Fang.  Fighting Fish and Two-Stack Sortable Permutations.  Séminaire Lotharingien de Combinatoire, 80B:P7, 2018.  Proceedings of the 30th International Conference on “Formal Power Series and Algebraic Combinatorics", July 16 - 20, 2018, Dartmouth College, Hanover, USA.
  * [40] Stefan Felsner.  Lattice structures from planar graphs.  Electron. J. Comb., 11(1), 2004.
  * [41] Stefan Felsner, Éric Fusy, Marc Noy, and David Orden.  Bijections for Baxter families and related objects.  J. Comb. Theory, Ser. A, 118(3):993–1020, 2011.
  * [42] Stefan Felsner and Lorenz Wernisch.  Markov chains for linear extensions, the two-dimensional case.  In Proc. ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 239–247. ACM/SIAM, 1997.
  * [43] Philippe Flajolet and Robert Sedgewick.  Analytic Combinatorics.  Cambridge University Press, 2009.
  * [44] Jean Françon and Gérard Viennot.  Permutations selon leurs pics, creux, doubles montées et double descentes, nombres d’Euler et nombres de Genocchi.  Discret. Math., 28(1):21–35, 1979.
  * [45] Ryo Fujimaki, Youhei Inoue, and Toshihiko Takahashi.  An asymptotic estimate of the numbers of rectangular drawings or floorplans.  In International Symposium on Circuits and Systems (ISCAS 2009), 24-17 May 2009, Taipei, Taiwan, pages 856–859. IEEE, 2009.
  * [46] Ryo Fujimaki and Toshihiko Takahashi.  A surjective mapping from permutations to room-to-room floorplans.  IEICE Trans. Fundam. Electron. Commun. Comput. Sci., 90-A(4):823–828, 2007.
  * [47] Éric Fusy.  Transversal structures on triangulations: A combinatorial study and straight-line drawings.  Discret. Math., 309(7):1870–1894, 2009.
  * [48] Éric Fusy, Erkan Narmanli, and Gilles Schaeffer.  Enumeration of corner polyhedra and 3-connected Schnyder labelings.  Electron. J. Comb., 30(2), 2023.
  * [49] Éric Fusy, Erkan Narmanli, and Gilles Schaeffer.  On the enumeration of plane bipolar posets and transversal structures.  Europ. J. Combin., 116:103870, 2024.
  * [50] Emmanuel Guitter, Charlotte F. Kristjansen, and Jacob L. Nielsen.  Hamiltonian cycles on random Eulerian triangulations.  Nuclear Phys. B, 546(3):731–750, 1999.
  * [51] Youhei Inoue, Toshihiko Takahashi, and Ryo Fujimaki.  Counting rectangular drawings or floorplans in polynomial time.  IEICE Trans. Fundam. Electron. Commun. Comput. Sci., 92-A(4):1115–1120, 2009.
  * [52] Richard Kenyon, Jason Miller, Scott Sheffield, and David B. Wilson.  Bipolar orientations on planar maps and $\mathrm{SLE}_{12}$.  Annals of Prob., 47(3):1240–1269, 2019.
  * [53] Shirley Law and Nathan Reading.  The Hopf algebra of diagonal rectangulations.  J. Comb. Theory, Ser. A, 119(3):788–824, 2012.
  * [54] Jean-Louis Loday.  Realization of the Stasheff polytope.  Archiv d. Math., 83(3):267–278, 2004.
  * [55] Colin L. Mallows.  Baxter permutations rise again.  Journal of Combinatorial Theory, Series A, 27(3):394–396, 1979\.
  * [56] Emily Meehan.  Baxter posets.  Electon. J. Comb., 26(3):3, 2019.
  * [57] Emily Meehan.  The Hopf algebra of generic rectangulations.  arXiv:1903.09874, 2019.
  * [58] Arturo Merino and Torsten Mütze.  Personal communication, 2021.
  * [59] Arturo Merino and Torsten Mütze.  Combinatorial generation via permutation languages. III. Rectangulations.  Discrete Comput. Geom., 70(1):51–122, 2023.
  * [60] Folkert Müller-Hoissen, Jean Marcel Pallo, and Jim Stasheff, editors.  Associahedra, Tamari Lattices and Related Structures.  Progress in Mathematics. Birkhäuser, 2012.
  * [61] The On-Line Encyclopedia of Integer Sequences.  Published electronically at http://oeis.org/.
  * [62] Vincent Pilaud and Francisco Santos.  Quotientopes.  Bull. London Math. Soc., 51(3):406–420, 2019.
  * [63] Lionel Pournin.  The diameter of associahedra.  Adv. Math., 259:13–42, 2014.
  * [64] James Propp.  Lattice structure for orientations of graphs.  arXiv:math/0209005, 1993.
  * [65] Nathan Reading.  Lattice congruences of the weak order.  Order, 21(4):315–344, 2004.
  * [66] Nathan Reading.  Generic rectangulations.  Europ. J. Comb., 33(4):610–623, 2012.
  * [67] Sadiq M. Sait and Habib Youssef.  VLSI physical design automation: Theory and practice.  IEEE, 1995.
  * [68] Gilles Schaeffer.  Conjugaison d’arbres et cartes combinatoires aléatoires.  PhD thesis, Bordeaux 1, 1998.
  * [69] Ernst Schröder.  Vier combinatorische Probleme.  Zeitschrift für Mathematik und Physik, Band 15:361–376, 1870.
  * [70] Zion Cien Shen and Chris C. N. Chu.  Bounds on the number of slicing, mosaic, and general floorplans.  IEEE Trans. Comput. Aided Des. Integr. Circuits Syst., 22(10):1354–1361, 2003.
  * [71] Naveed A. Sherwani.  Algorithms for VLSI physical design automation.  Kluwer, 1993.
  * [72] Daniel D. Sleator, Robert E. Tarjan, and William P. Thurston.  Rotation distance, triangulations, and hyperbolic geometry.  J. Amer. Math. Soc., 1(3):647–681, 1988.
  * [73] Richard P. Stanley.  Hipparchus, Plutarch, Schröder, and Hough.  The American Mathematical Monthly, 104(4):344–350, 1997.
  * [74] Richard P. Stanley.  Catalan Numbers.  Cambridge University Press, 2015.
  * [75] James Dillon Stasheff.  Homotopy associativity of $H$-spaces. I, II.  Trans. Amer. Math. Soc., 108:293–312, 1963.  108 (1963), 275-292; ibid.
  * [76] John Philip Steadman.  Architectural Morphology.  Pion, 1983.
  * [77] Toshihiko Takahashi and Ryo Fujimaki.  Fujimaki–Takahashi squeeze: Linear time construction of constraint graphs of floorplan for a given permutation.  IEICE Trans. Fundam. Electron. Commun. Comput. Sci., 91-A(4):1071–1076, 2008.
  * [78] Toshihiko Takahashi, Ryo Fujimaki, and Youhei Inoue.  A $(4n-4)$-bit representation of a rectangular drawing or floorplan.  In Hung Q. Ngo, editor, Proc. COCOON 2009, volume 5609 of LNCS, pages 47–55. Springer, 2009.
  * [79] Dov Tamari.  Monoïdes préordonnés et chaînes de Malcev.  Université de Paris, Paris, 1951.  Thèse.
  * [80] William Thomas Tutte.  A census of planar triangulations.  Canadian J. Math., 14:21–38, 1962.
  * [81] William Thomas Tutte.  A census of planar maps.  Canadian J. Math., 15:249–271, 1963.
  * [82] Marc J. van Kreveld and Bettina Speckmann.  On rectangular cartograms.  Comput. Geom., 37(3):175–187, 2007.
  * [83] Gérard Viennot.  A bijective proof for the number of Baxter permutations.  Séminaire Lotharingien de Combinatoire, 1-3:28–29, 1981\.
  * [84] Julian West.  Generating trees and the Catalan and Schröder numbers.  Discret. Math., 146(1-3):247–262, 1995.
  * [85] Bo Yao, Hongyu Chen, Chung-Kuan Cheng, and Ronald L. Graham.  Floorplan representations: Complexity and connections.  ACM Trans. Design Autom. Electron. Syst., 8(1):55–80, 2003.
  * [86] Doron Zeilberger.  A proof of Julian West’s conjecture that the number of two-stack-sortable permutations of length $n$ is $2(3n)!/((n+1)!(2n+1)!)$.  Discret. Math., 102(1):85–93, 1992.